-- LU_DECOMP.ADA  Decompose a matrix A into a Lower triangular matrix
--                times an Upper triangluar matrix  A = L * U
--                where the diagonal of L is all ones.
--                Or more generally  P * A = L * U where P is a permutation
--                matrix (a permutation of the rows of the Identity matrix)
--                represented as an integer vector of subscripts
--
--  Methods:  1) simple LU that will not work for some matricies
--            2) from newmat, a translation from C
--            3) from Introduction to Algorithms P 759
--            4) from Numerical Recipies in Fortran P 38
--
--            Array_Index_Error is raised if matrix not square and also if
--            A, L and U  row and column ranges are not the same.
--
--            There is a test procedure at the end of this file, Test_LU_Decomp
--
-- basis of algorithm:
--
-- write out L U = A  and perform matrix multiply
--
-- | 1.0    0.0    0.0    0.0    |   | U(1,1) U(1,2) U(1,3) U(1,4) | 
-- | L(2,1) 1.0    0.0    0.0    | * | 0.0    U(2,2) U(2,3) U(2,4) | =
-- | L(3,1) L(3,2) 1.0    0.0    |   | 0.0    0.0    U(3,3) U(3,4) |
-- | L(4,1) L(4,2) L(4,3) 1.0    |   | 0.0    0.0    0.0    U(4,4) |
--
-- | A(1,1) A(1,2) A(1,3) A(1,4) |   | U(1,1)         U(1,2)
-- | A(2,1) A(2,2) A(2,3) A(2,4) | = | L(2,1)*U(1,1)  L(2,1)*U(1,2)+U(2,2)
-- | A(3,1) A(3,2) A(3,3) A(3,4) |   | L(3,1)*U(1,1)  L(3,1)*U(1,2)+L(3,2)*U(2,2)
-- | A(4,1) A(4,2) A(4,3) A(4,4) |   | L(4,1)*U(1,1)  L(4,1)*U(1,2)+L(4,2)*U(2,2)
--
-- U(1,3)
-- L(2,1)*U(1,3)+U(2,3)
-- L(3,1)*U(1,3)+L(3,2)*U(2,3)+U(3,3)
-- L(4,1)*U(1,3)+L(4,2)*U(2,3)+L(4,3)*U(3,3)
--
-- U(1,4)                                           |
-- L(2,1)*U(1,4)+U(2,4)                             |
-- L(3,1)*U(1,4)+L(3,2)*U(2,4)+U(3,4)               |
-- L(4,1)*U(1,4)+L(4,2)*U(2,4)+L(4,3)*U(3,4)+U(4,4) |
--
-- order of computation: (just solving each equation above)
--
-- U(1,1) := A(1,1)-Sum   Sum = 0.0
-- U(1,2) := A(1,2)-Sum   Sum = 0.0
-- U(1,3) := A(1,3)-Sum   Sum = 0.0
-- U(1,4) := A(1,4)-Sum   Sum = 0.0
--
-- L(2,1) := (A(2,1)-Sum)/U(1,1) Sum = 0.0
-- L(3,1) := (A(3,1)-Sum)/U(1,1) Sum = 0.0
-- L(4,1) := (A(4,1)-Sum)/U(1,1) Sum = 0.0
--
-- U(2,2) := A(2,2)-Sum          Sum = L(2,1)*U(1,2)
-- U(2,3) := A(2,3)-Sum          Sum = L(3,1)*U(1,2)
-- U(2,4) := A(2,4)-Sum          Sum = L(4,1)*U(1,2)
--
-- L(3,2) := (A(3,2)-Sum)/U(2,2) Sum = L(3,1)*U(1,2)
-- L(4,2) := (A(4,2)-Sum)/U(2,2) Sum = L(4,1)*U(1,2)
--
-- U(3,3) := A(3,3)-Sum          Sum = L(3,1)*U(1,3)+L(3,2)*U(2,3)
-- U(3,4) := A(3,4)-Sum          Sum = L(3,1)*U(1,4)+L(3,2)*U(2,4)
--
-- L(4,3) := (A(4,3)-Sum)/U(3,3) Sum = L(4,1)*U(1,3)+L(4,2)*U(2,3)
--
-- U(4,4) := A(4,4)-Sum     Sum = L(4,1)*U(1,4)+L(4,2)*U(2,4)+L(4,3)+U(3,4)
--
-- observations: 
--      since L is lower diagonal, the diagonal and upper L need never be used.
--      since U is upper diagonal and diagonal, lower U need never be used.
--
-- programmed 3/5/94 by Jon S. Squire

with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
Procedure LU_DECOMP1( A: Complex_Matrix; L, U : in out Complex_Matrix) is
  Sum : Complex;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if A'Length(1) /= A'Length(2) or -- remove later by offsetting
     A'Length(1) /= L'Length(1) or
     A'Length(1) /= L'Length(2) or
     A'Length(1) /= U'Length(1) or
     A'Length(1) /= U'Length(2) or
     A'First(1)  /= A'First(2)  or
     A'First(1)  /= L'First(1)  or
     A'First(1)  /= L'First(2)  or
     A'First(1)  /= U'First(1)  or
     A'First(1)  /= U'First(2)  then
    raise Array_Index_Error;
  end if;

  L := (others=>(others=>(0.0,0.0)));
  U := (others=>(others=>(0.0,0.0)));

  for I in A'RANGE(1) loop -- initialize diagonal
    L(I,I) := (1.0,0.0);
  end loop;                -- only I <= J elements of U(I,J) to be computed
                           -- only I >  J elements of L(I,J) to be computed 

  for I in A'RANGE(1) loop
    for J in I..A'LAST(1) loop -- do a row of Upper
      SUM := A(I,J);
      for K in 1 .. I-1 loop 
        SUM := SUM - L(I,K) * U(K,J);
      end loop;
      U(I,J) := SUM;  -- no divide because L(I,I) is always one
    end loop;
    for J in I+1 .. A'LAST(2) loop  -- do a column of Lower
      SUM := A(J,I);
      for K in 1..I-1 loop
        SUM := SUM - L(J,K) * U(K,I);
      end loop;
      L(J,I) := SUM/U(I,I);
    end loop;
  end loop;
end LU_DECOMP1;

with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
function LU_SOLVE1 ( L : COMPLEX_MATRIX ;
                     U : COMPLEX_MATRIX ;
                     Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
  X : COMPLEX_VECTOR ( Y'RANGE ) ;
  B : COMPLEX_VECTOR ( Y'RANGE ) ;
  SUM : COMPLEX ;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if L'Length(1) /= L'Length(2) or -- remove later by offsetting
     L'Length(1) /= U'Length(1) or
     L'Length(1) /= U'Length(2) or
     L'Length(1) /= Y'Length    or
     L'First(1)  /= L'First(2)  or
     L'First(1)  /= U'First(1)  or
     L'First(1)  /= U'First(2)  or
     L'First(1)  /= Y'First     then
    raise Array_Index_Error;
  end if;

  for I in Y'Range loop  --   do i=1,n
    SUM := (0.0,0.0) ;
    for J in Y'First .. I-1 loop --     do j=1,i-1
      SUM := SUM + L(I,J) * B(J) ;
    end loop ;
    B(I) := Y(I) - SUM ;
  end loop ;

  for I in reverse Y'Range loop --   do i=n,1,-1
    SUM := (0.0,0.0) ;
    for J in I+1 .. Y'Last loop --     do j=i+1,n
      SUM := SUM + U(I,J) * X(J) ;
    end loop ;
    X(I) := ( B(I) - SUM ) / U(I,I) ;
  end loop ;
  return X;
end LU_SOLVE1 ;

with TEXT_IO; use TEXT_IO;
with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with COMPLEX_ARRAYS_IO; use COMPLEX_ARRAYS_IO;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
procedure LU_DECOMP2 ( A : COMPLEX_MATRIX ;
                       L : in out COMPLEX_MATRIX ;
                       U : in out COMPLEX_MATRIX ;
                       P : in out INTEGER_VECTOR) is
  BIG : FLOAT;
  BIG_INDEX : INTEGER;
  TEMP : FLOAT;
  SUM : COMPLEX;
  T : COMPLEX;
  SINGULAR_ERROR : exception;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if A'Length(1) /= A'Length(2) or -- remove later by offsetting
     A'Length(1) /= L'Length(1) or
     A'Length(1) /= L'Length(2) or
     A'Length(1) /= U'Length(1) or
     A'Length(1) /= U'Length(2) or
     A'Length(1) /= P'Length    or
     A'First(1)  /= A'First(2)  or
     A'First(1)  /= L'First(1)  or
     A'First(1)  /= L'First(2)  or
     A'First(1)  /= U'First(1)  or
     A'First(1)  /= U'First(2)  or
     A'First(1)  /= P'First     then
    raise Array_Index_Error;
  end if;

  L := (others=>(others=>(0.0,0.0)));
  U := (others=>(others=>(0.0,0.0)));
  for I in A'RANGE(1) loop -- initialize diagonal
    L(I,I) := (1.0,0.0);
    P(I) := I;
  end loop;                -- only I <= J elements of U(I,J) to be computed
                           -- only I >  J elements of L(I,J) to be computed 

  for I in A'RANGE(1) loop

    BIG := abs A(P(I),I);               -- find biggest pivot in rest of column
    for J in I+1..A'LAST(2) loop
      TEMP := abs A(P(J),I);
      if TEMP > BIG then
        BIG := TEMP;
        BIG_INDEX := P(J);
        P(J) := P(I);
        P(I) := BIG_INDEX;
      end if;
    end loop;
    if BIG < FLOAT'EPSILON ** 2 then -- singular
      PUT_LINE("A matrix sigular, detected at row " & INTEGER'IMAGE(I));
      raise SINGULAR_ERROR;
    end if;

    for J in I..A'LAST(1) loop -- do a row of Upper
      SUM := A(P(I),J);
      for K in 1 .. I-1 loop 
        SUM := SUM - L(P(I),K) * U(K,J);
      end loop;
      U(P(I),J) := SUM;  -- no divide because L(I,I) is always one
    end loop;

    for J in I+1 .. A'LAST(2) loop  -- do a Column of Lower
      SUM := A(P(J),I);
      for K in 1..I-1 loop
        SUM := SUM - L(P(J),K) * U(K,I);
      end loop;
      L(J,I) := SUM/U(P(I),I);
    end loop;

    PUT_LINE("Status on iteration I = " & INTEGER'IMAGE(I));

    PUT("row positions = ");
    for K in A'RANGE(1) loop
      PUT(INTEGER'IMAGE(P(K)));
    end loop;
    NEW_LINE;
    PUT_LINE("L");
    PUT(L);
    PUT_LINE("U");
    PUT(U);
    NEW_LINE;

  end loop;
end LU_DECOMP2;

with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
function LU_SOLVE2 ( L : COMPLEX_MATRIX ;
                     U : COMPLEX_MATRIX ;
                     P : INTEGER_VECTOR ;
                     Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
  X : COMPLEX_VECTOR ( Y'RANGE ) ;
  B : COMPLEX_VECTOR ( Y'RANGE ) ;
  SUM : COMPLEX ;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if L'Length(1) /= L'Length(2) or -- remove later by offsetting
     L'Length(1) /= U'Length(1) or
     L'Length(1) /= U'Length(2) or
     L'Length(1) /= P'Length    or
     L'Length(1) /= Y'Length    or
     L'First(1)  /= L'First(2)  or
     L'First(1)  /= U'First(1)  or
     L'First(1)  /= U'First(2)  or
     L'First(1)  /= P'First     or
     L'First(1)  /= Y'First     then
    raise Array_Index_Error;
  end if;

  for I in Y'Range loop  --   do i=1,n
    SUM := (0.0,0.0) ;
    for J in Y'First .. I-1 loop --     do j=1,i-1
      SUM := SUM + L(I,J) * B(J) ;
    end loop ;
    B(I) := Y(P(I)) - SUM ;
  end loop ;

  for I in reverse Y'Range loop --   do i=n,1,-1
    SUM := (0.0,0.0) ;
    for J in I+1 .. Y'Last loop --     do j=i+1,n
      SUM := SUM + U(I,J) * X(J) ;
    end loop ;
    X(I) := ( B(I) - SUM ) / U(I,I) ;
  end loop ;
  return X;
end LU_SOLVE2 ;

with TEXT_IO; use TEXT_IO;
with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
procedure LU_DECOMP3 ( A : COMPLEX_MATRIX ;
                       L : in out COMPLEX_MATRIX ;
                       U : in out COMPLEX_MATRIX ;
                       P : in out INTEGER_VECTOR) is
  AA : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) := A;
  ITEMP : INTEGER;
  KK    : INTEGER;
  TEMP  : COMPLEX;
  BIG   : FLOAT;
  SINGULAR_ERROR : exception;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if A'Length(1) /= A'Length(2) or -- remove later by offsetting
     A'Length(1) /= L'Length(1) or
     A'Length(1) /= L'Length(2) or
     A'Length(1) /= U'Length(1) or
     A'Length(1) /= U'Length(2) or
     A'Length(1) /= P'Length    or
     A'First(1)  /= A'First(2)  or
     A'First(1)  /= L'First(1)  or
     A'First(1)  /= L'First(2)  or
     A'First(1)  /= U'First(1)  or
     A'First(1)  /= U'First(2)  or
     A'First(1)  /= P'First     then
    raise Array_Index_Error;
  end if;

  L := (others=>(others=>(0.0,0.0)));
  U := (others=>(others=>(0.0,0.0)));
  for I in P'RANGE loop
    P(I) := I;
    L(I,I) := (1.0,0.0);
  end loop;
  for K in A'FIRST(1)..A'LAST(1)-1 loop
    BIG := 0.0;
    for I in K..A'LAST(1) loop
      if abs AA(I,K) > BIG then
        BIG := abs AA(I,K);
        KK := I;
      end if;
    end loop;
    if BIG = 0.0 then
      raise SINGULAR_ERROR;
    end if;
    ITEMP := P(K);
    P(K) := P(KK);
    P(KK) := ITEMP;
    for I in A'RANGE(1) loop
      TEMP := AA(K,I);
      AA(K,I) := AA(KK,I);
      AA(KK,I) := TEMP;
    end loop;
    for I in K+1..A'LAST(1) loop
      AA(I,K) := AA(I,K) / AA(K,K);
      for J in K+1..A'LAST(1) loop
        AA(I,J) := AA(I,J) - AA(I,K) * AA(K,J);
      end loop;
    end loop;
  end loop;
  for I in A'RANGE(1) loop  -- copy from A to U
    for J in I..A'LAST(1) loop
      U(I,J) := AA(I,J);
    end loop;
  end loop;
  for I in A'FIRST(1)+1..A'LAST(1) loop  -- copy from A to L
    for J in 1..I-1 loop
      L(I,J) := AA(I,J);
    end loop;
  end loop;
end LU_DECOMP3;

with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
function LU_SOLVE3 ( L : COMPLEX_MATRIX ;
                     U : COMPLEX_MATRIX ;
                     P : INTEGER_VECTOR ;
                     Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
  X : COMPLEX_VECTOR ( Y'RANGE ) ;
  B : COMPLEX_VECTOR ( Y'RANGE ) ;
  SUM : COMPLEX ;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if L'Length(1) /= L'Length(2) or -- remove later by offsetting
     L'Length(1) /= U'Length(1) or
     L'Length(1) /= U'Length(2) or
     L'Length(1) /= P'Length    or
     L'Length(1) /= Y'Length    or
     L'First(1)  /= L'First(2)  or
     L'First(1)  /= U'First(1)  or
     L'First(1)  /= U'First(2)  or
     L'First(1)  /= P'First     or
     L'First(1)  /= Y'First     then
    raise Array_Index_Error;
  end if;

  for I in Y'Range loop  --   do i=1,n
    SUM := (0.0,0.0) ;
    for J in Y'First .. I-1 loop --     do j=1,i-1
      SUM := SUM + L(I,J) * B(J) ;
    end loop ;
    B(I) := Y(P(I)) - SUM ;
  end loop ;

  for I in reverse Y'Range loop --   do i=n,1,-1
    SUM := (0.0,0.0) ;
    for J in I+1 .. Y'Last loop --     do j=i+1,n
      SUM := SUM + U(I,J) * X(J) ;
    end loop ;
    X(I) := ( B(I) - SUM ) / U(I,I) ;
  end loop ;
  return X;
end LU_SOLVE3 ;


with TEXT_IO; use TEXT_IO;
with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with REAL_ARRAYS; use REAL_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
procedure LU_DECOMP4 ( A : COMPLEX_MATRIX ;
                       L : in out COMPLEX_MATRIX ;
                       U : in out COMPLEX_MATRIX ;
                       P : in out INTEGER_VECTOR) is
  D     : COMPLEX := (1.0,0.0); -- determinant
  AAMAX : FLOAT;
  SUM   : COMPLEX;
  TEMP  : COMPLEX;
  DUM   : FLOAT;
  IMAX  : INTEGER;
  AA    : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) := A;
  VV    : REAL_VECTOR(A'RANGE(1));
  SINGULAR_ERROR : exception;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if A'Length(1) /= A'Length(2) or -- remove later by offsetting
     A'Length(1) /= L'Length(1) or
     A'Length(1) /= L'Length(2) or
     A'Length(1) /= U'Length(1) or
     A'Length(1) /= U'Length(2) or
     A'Length(1) /= P'Length    or
     A'First(1)  /= A'First(2)  or
     A'First(1)  /= L'First(1)  or
     A'First(1)  /= L'First(2)  or
     A'First(1)  /= U'First(1)  or
     A'First(1)  /= U'First(2)  or
     A'First(1)  /= P'First     then
    raise Array_Index_Error;
  end if;

  L := (others=>(others=>(0.0,0.0)));
  U := (others=>(others=>(0.0,0.0)));
  for I in A'RANGE(1) loop -- 12
    L(I,I) := (1.0,0.0);
    AAMAX := 0.0;
    for J in A'RANGE(2) loop
      if abs AA(I,J) > AAMAX then
        AAMAX := abs AA(I,J);
      end if;
    end loop;
    if AAMAX = 0.0 then
      raise SINGULAR_ERROR;
    end if;
    VV(I) := 1.0 / AAMAX;
  end loop;  -- 12
  for J in A'RANGE(2) loop -- 19
    for I in 1..J-1 loop -- 14
      SUM := AA(I,J);
      for K in 1..I-1 loop -- 13
        SUM := SUM - AA(I,K) * AA(K,J);
      end loop; -- 13
      AA(I,J) := SUM;
    end loop; -- 14
    AAMAX := 0.0;
    for I in J..A'LAST(1) loop -- 16
      SUM := AA(I,J);
      for K in 1..J-1 loop -- 15
        SUM := SUM - AA(I,K) * AA(K,J);
      end loop; -- 15
      AA(I,J) := SUM;
      DUM := VV(I) * abs SUM;
      if DUM >= AAMAX then
        IMAX := I;
        AAMAX := DUM;
      end if;
    end loop; -- 16
    if J /= IMAX then
      for K in A'RANGE(1) loop -- 17
        TEMP := AA(IMAX,K);
        AA(IMAX,K) := AA(J,K);
        AA(J,K) := TEMP;
      end loop; -- 17
      D := -D;
      VV(IMAX) := VV(J);
    end if;
    P(J) := IMAX;
    if abs AA(J,J) = 0.0 then
      raise SINGULAR_ERROR;
    end if;
    if J /= A'LAST(1) then
      TEMP := 1.0 / AA(J,J);
      for I in J+1..A'LAST(1) loop -- 18
        AA(I,J) := AA(I,J) * TEMP;
      end loop; -- 18
    end if;
  end loop; -- 19
  for I in A'RANGE(1) loop  -- copy from A to U
    for J in I..A'LAST(1) loop
      U(I,J) := AA(I,J);
    end loop;
  end loop;
  for I in A'FIRST(1)+1..A'LAST(1) loop  -- copy from A to L
    for J in 1..I-1 loop
      L(I,J) := AA(I,J);
    end loop;
  end loop;
end LU_DECOMP4;

with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
function LU_SOLVE4 ( L : COMPLEX_MATRIX ;
                     U : COMPLEX_MATRIX ;
                     P : INTEGER_VECTOR ;
                     Y : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
  X : COMPLEX_VECTOR ( Y'RANGE ) ;
  B : COMPLEX_VECTOR ( Y'RANGE ) ;
  SUM : COMPLEX ;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  if L'Length(1) /= L'Length(2) or -- remove later by offsetting
     L'Length(1) /= U'Length(1) or
     L'Length(1) /= U'Length(2) or
     L'Length(1) /= P'Length    or
     L'Length(1) /= Y'Length    or
     L'First(1)  /= L'First(2)  or
     L'First(1)  /= U'First(1)  or
     L'First(1)  /= U'First(2)  or
     L'First(1)  /= P'First     or
     L'First(1)  /= Y'First     then
    raise Array_Index_Error;
  end if;

  for I in Y'Range loop  --   do i=1,n
    SUM := (0.0,0.0) ;
    for J in Y'First .. I-1 loop --     do j=1,i-1
      SUM := SUM + L(I,J) * B(J) ;
    end loop ;
    B(I) := Y(P(I)) - SUM ;
  end loop ;

  for I in reverse Y'Range loop --   do i=n,1,-1
    SUM := (0.0,0.0) ;
    for J in I+1 .. Y'Last loop --     do j=i+1,n
      SUM := SUM + U(I,J) * X(J) ;
    end loop ;
    X(I) := ( B(I) - SUM ) / U(I,I) ;
  end loop ;
  return X;
end LU_SOLVE4 ;


with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
function MAKE_PERM( P : INTEGER_VECTOR ) return COMPLEX_MATRIX is
  PP : COMPLEX_MATRIX(P'RANGE,P'RANGE) := (others=>(others=>(0.0,0.0)));
begin
  for I in P'RANGE loop
    PP(I,P(I)) := (1.0,0.0);
  end loop;
  return PP;
end MAKE_PERM;

with LU_DECOMP1;
with LU_DECOMP2;
with LU_DECOMP3;
with LU_DECOMP4;
with LU_SOLVE1;
with LU_SOLVE2;
with LU_SOLVE3;
with LU_SOLVE4;
with MAKE_PERM;
with COMPLEX_TYPES; use COMPLEX_TYPES;
with COMPLEX_ARRAYS; use COMPLEX_ARRAYS;
with COMPLEX_ARRAYS_IO; use COMPLEX_ARRAYS_IO;
with INTEGER_ARRAYS; use INTEGER_ARRAYS;
with INTEGER_ARRAYS_IO; use INTEGER_ARRAYS_IO;
with TEXT_IO; use TEXT_IO;
with FLOAT_TEXT_IO; use FLOAT_TEXT_IO;
procedure TEST_LU_DECOMP is
  A : COMPLEX_MATRIX(1..4,1..4);
  L : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((1.0,0.0),(0.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((3.0,0.0),(1.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((1.0,1.0),(4.0,0.0),(1.0,0.0),(0.0,0.0)),
           ((1.0,2.0),(3.0,3.0),(2.0,0.0),(1.0,0.0)));
  L2 : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((1.0,0.0), (0.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((0.4,0.0), (1.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((-0.2,0.0),(0.5,0.0),(1.0,0.0),(0.0,0.0)),
           ((0.6,0.0), (0.0,0.0),(0.4,0.0),(1.0,0.0)));
  L3 : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((1.0,0.0),  (0.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((0.25,0.0), (1.0,0.0),(0.0,0.0),(0.0,0.0)),
           ((0.5,0.0),  (0.0,0.0),(1.0,0.0),(0.0,0.0)),
           ((0.75,0.0), (0.0,0.0),(0.0,0.0),(1.0,0.0)));

  U : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((2.0,1.0),(2.0,0.0),(0.0,3.0),(4.0,4.0)),
           ((0.0,0.0),(1.0,2.0),(3.0,0.0),(2.0,2.0)),
           ((0.0,0.0),(0.0,0.0),(1.0,3.0),(1.0,2.0)),
           ((0.0,0.0),(0.0,0.0),(0.0,0.0),(2.0,2.0)));
  U2 : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((5.0,0.0),( 5.0,0.0),(4.0,0.0),( 2.0,0.0)),
           ((0.0,0.0),(-2.0,0.0),(0.4,0.0),(-0.2,0.0)),
           ((0.0,0.0),( 0.0,0.0),(4.0,0.0),(-0.5,0.0)),
           ((0.0,0.0),( 0.0,0.0),(0.0,0.0),(-3.0,0.0)));
  U3 : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)) :=
          (((4.0,0.0),( 4.0,0.0),(4.0,0.0),(4.0,0.0)),
           ((0.0,0.0),( 1.0,0.0),(2.0,0.0),(3.0,0.0)),
           ((0.0,0.0),( 0.0,0.0),(1.0,0.0),(2.0,0.0)),
           ((0.0,0.0),( 0.0,0.0),(0.0,0.0),(1.0,0.0)));

  AA : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2));
  AS : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2));
  A2 : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2)):=
          ((( 2.0,0.0),( 0.0,0.0),(2.0,0.0),( 0.6,0.0)),
           (( 3.0,0.0),( 3.0,0.0),(4.0,0.0),(-2.0,0.0)),
           (( 5.0,0.0),( 5.0,0.0),(4.0,0.0),( 2.0,0.0)),
           ((-1.0,0.0),(-2.0,0.0),(3.4,0.0),(-1.0,0.0)));
  V : COMPLEX_VECTOR(A'RANGE(1));
  LS : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2));
  US : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2));
  P  : INTEGER_VECTOR(A'RANGE(1));
  PP : COMPLEX_MATRIX(A'RANGE(1),A'RANGE(2));
  X  : COMPLEX_VECTOR(A'RANGE(1));
  Y  : COMPLEX_VECTOR(A'RANGE(1)) := ((4.0,0.0),(1.0,0.0),(3.0,0.0),(2.0,0.0));
  YY : COMPLEX_VECTOR(A'RANGE(1));

  SINGULAR_ERROR : exception;
  ARRAY_INDEX_ERROR : exception renames COMPLEX_ARRAYS.ARRAY_INDEX_ERROR;
begin
  PUT_LINE("begin TEST_LU_DECOMP");
  PUT_LINE("Case 1");
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  A := L * U;
  PUT_LINE("A");
  PUT(A);
  NEW_LINE;
  PUT_LINE("LU_DECOMP1 has no interchange");

  LU_DECOMP1(A, L, U);

  PUT_LINE("check LU_DECOMP1");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  NEW_LINE;
  AA := A - L * U;
  PUT_LINE("AA := A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP2  interchange on largest pivot");

  LU_DECOMP2(A, L, U, P);

  PUT_LINE("check LU_DECOMP2");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP3  interchange on largest pivot");

  LU_DECOMP3(A, L, U, P);

  PUT_LINE("check LU_DECOMP3");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  NEW_LINE;
  PUT_LINE("Y");
  PUT(Y);
  X := LU_SOLVE3(L,U,P,Y);
  PUT_LINE("X = LU_SOLVE3(L,U,Y)");
  PUT(X);
  YY := Y - A * X ;
  PUT_LINE("YY := Y - A * X should be zero");
  PUT(YY);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP4  interchange on largest pivot");

  LU_DECOMP4(A, L, U, P);

  PUT_LINE("check LU_DECOMP4");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);

  PUT_LINE("Case 2");
  PUT_LINE("L2");
  PUT(L2);
  PUT_LINE("U2");
  PUT(U2);
  A := L2 * U2;
  PUT_LINE("A");
  PUT(A);
  NEW_LINE;
  PUT_LINE("LU_DECOMP1 has no interchange");

  LU_DECOMP1(A, L, U);

  PUT_LINE("check LU_DECOMP1");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  NEW_LINE;
  AA := A - L * U;
  PUT_LINE("AA := A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP2  interchange on largest pivot");
  begin
  LU_DECOMP2(A, L, U, P);

  PUT_LINE("check LU_DECOMP2");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  exception
    when others =>
      PUT_LINE("Died with exception");
  end;
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP3  interchange on largest pivot");

  LU_DECOMP3(A, L, U, P);

  PUT_LINE("check LU_DECOMP3");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  PUT_LINE("Y");
  PUT(Y);
  X := LU_SOLVE3(L,U,P,Y);
  PUT_LINE("X = LU_SOLVE3(L,U,Y)");
  PUT(X);
  YY := Y - A * X ;
  PUT_LINE("YY := Y - A * X should be zero");
  PUT(YY);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP4  interchange on largest pivot");

  LU_DECOMP4(A, L, U, P);

  PUT_LINE("check LU_DECOMP4");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);


  PUT_LINE("Case 3");
  PUT_LINE("L3");
  PUT(L3);
  PUT_LINE("U3");
  PUT(U3);
  A := L3 * U3;
  PUT_LINE("A");
  PUT(A);
  NEW_LINE;
  PUT_LINE("LU_DECOMP1 has no interchange");

  LU_DECOMP1(A, L, U);

  PUT_LINE("check LU_DECOMP1");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  NEW_LINE;
  AA := A - L * U;
  PUT_LINE("AA := A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP2  interchange on largest pivot");
  begin
  LU_DECOMP2(A, L, U, P);

  PUT_LINE("check LU_DECOMP2");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  exception
    when others =>
      PUT_LINE("Died with exception");
  end;
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP3  interchange on largest pivot");

  LU_DECOMP3(A, L, U, P);

  PUT_LINE("check LU_DECOMP3");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  PUT_LINE("Y");
  PUT(Y);
  X := LU_SOLVE3(L,U,P,Y);
  PUT_LINE("X = LU_SOLVE3(L,U,Y)");
  PUT(X);
  YY := Y - A * X ;
  PUT_LINE("YY := Y - A * X should be zero");
  PUT(YY);
  NEW_LINE(2);

  PUT_LINE("LU_DECOMP4  interchange on largest pivot");

  LU_DECOMP4(A, L, U, P);

  PUT_LINE("check LU_DECOMP4");
  NEW_LINE;
  PUT_LINE("L");
  PUT(L);
  PUT_LINE("U");
  PUT(U);
  PUT_LINE("P");
  PUT(P);
  NEW_LINE;
  PP := MAKE_PERM(P);
  PUT_LINE("PP");
  PUT(PP);
  AA := PP * A - L * U;
  PUT_LINE("AA := P * A - L * U should be zero");
  PUT(AA);
  NEW_LINE(2);


end TEST_LU_DECOMP;

--   int i,j; // C++ code from newmat8
--   real* vv=new real [nrows]; MatrixErrorNoSpace(vv);
--   real* a;
--   a=store;
--   for (i=0;i<nrows;i++)
--   {
--      real big=0.0;
--      j=nrows; while (j--) { real temp=fabs(*a++); if (temp > big) big=temp; }
--      if (big == 0.0) MatrixError("Singular matrix");
--     vv[i]=1.0/big;
--   }
--
--   real* aj=store;
--   for (j=0;j<nrows;j++)
--   {
--      real* ai=store;
--      for (i=0;i<j;i++)
--      {
--         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
--         int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
--         *ajx = sum; ai += nrows;
--      }
--
--      real big=0.0; int imax;
--      for (i=j;i<nrows;i++)
--      {
--         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
--         int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
--         *aix = sum; ai += nrows;
--         real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
--      }
--
--      if (j != imax)
--      {
--         real* amax=store+imax*nrows; real* ajx=store+j*nrows;
--         int k=nrows; while(k--) { real dum=*amax; *amax++=*ajx; *ajx++=dum; }
--         d=!d; vv[imax]=vv[j];
--      }
--
--      indx[j]=imax; ai=aj+j*nrows;
--      if (*ai == 0.0) MatrixError("Singular matrix");
--      real dum=1.0/(*ai);
--      i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
--
--      aj++;
--   }
--   delete [] vv;
--}