with text_io,integer_io;                 use text_io,integer_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Natural_Vectors;
with Standard_Natural_Vectors_io;        use Standard_Natural_Vectors_io;
with Standard_Floating_Vectors;
with Standard_Floating_Vectors_io;       use Standard_Floating_Vectors_io;
with Standard_Floating_Matrices;
with Standard_Floating_Matrices_io;      use Standard_Floating_Matrices_io;
with Standard_Random_Vectors;            use Standard_Random_Vectors;
with Standard_Random_Matrices;           use Standard_Random_Matrices;
with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;
with Standard_Floating_Least_Squares;    use Standard_Floating_Least_Squares;
with Standard_Complex_Vectors;
with Standard_Complex_Vectors_io;        use Standard_Complex_Vectors_io;
with Standard_Complex_Matrices;
with Standard_Complex_Matrices_io;       use Standard_Complex_Matrices_io;
with Standard_Complex_QR_Decomposition;  use Standard_Complex_QR_Decomposition;
with Standard_Complex_Least_Squares;     use Standard_Complex_Least_Squares;

procedure ts_qrd is

-- DESCRIPTION :
--   This program tests the implementation of the QR decomposition
--   and least squares approximation.

-- GENERAL TESTS ON QR DECOMPOSITION :

  function Extract_Upper_Triangular
                ( a : Standard_Floating_Matrices.Matrix )
                return Standard_Floating_Matrices.Matrix is

  -- DESCRIPTION :
  --   Returns the upper triangular part of the matrix a.

    res : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));

  begin
    for i in a'range(1) loop
      for j in a'first(2)..(i-1) loop
        res(i,j) := 0.0;
      end loop;
      for j in i..a'last(2) loop
        res(i,j) := a(i,j);
      end loop;
    end loop;
    return res;
  end Extract_Upper_Triangular;

  function Extract_Upper_Triangular
                ( a : Standard_Complex_Matrices.Matrix )
                return Standard_Complex_Matrices.Matrix is

  -- DESCRIPTION :
  --   Returns the upper triangular part of the matrix a.

    res : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));

  begin
    for i in a'range(1) loop
      for j in a'first(2)..(i-1) loop
        res(i,j) := Create(0.0);
      end loop;
      for j in i..a'last(2) loop
        res(i,j) := a(i,j);
      end loop;
    end loop;
    return res;
  end Extract_Upper_Triangular;

  function Differences ( a,b : in Standard_Floating_Matrices.Matrix )
                       return double_float is

  -- DESCRIPTION :
  --   Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.

    sum : double_float := 0.0;

  begin
    for i in a'range(1) loop
      for j in a'range(2) loop
        sum := sum + abs(a(i,j)-b(i,j));
      end loop;
    end loop;
    return sum;
  end Differences;

  function Differences ( a,b : in Standard_Complex_Matrices.Matrix )
                       return double_float is

  -- DESCRIPTION :
  --   Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.

    sum : double_float := 0.0;

  begin
    for i in a'range(1) loop
      for j in a'range(2) loop
        sum := sum + AbsVal(a(i,j)-b(i,j));
      end loop;
    end loop;
    return sum;
  end Differences;

  procedure Orthogonality ( q : in Standard_Floating_Matrices.Matrix ) is

  -- DESCRIPTION :
  --   Tests whether the columns are orthogonal w.r.t. each other.

    sum,ip : double_float;

  begin
    sum := 0.0;
    for j in q'range(2) loop
      for k in j+1..q'last(2) loop
        ip := 0.0;
        for i in q'range(1) loop
          ip := ip + q(i,j)*q(i,k);
        end loop;
        sum := sum + abs(ip);
      end loop;
    end loop;
    put("Orthogonality check : "); put(sum,3,3,3); new_line;
  end Orthogonality;

  procedure Orthogonality ( q : in Standard_Complex_Matrices.Matrix ) is

  -- DESCRIPTION :
  --   Tests whether the columns are orthogonal w.r.t. each other.

    sum : double_float := 0.0;
    ip : Complex_Number;

  begin
    for j in q'range(2) loop
      for k in j+1..q'last(2) loop
        ip := Create(0.0);
        for i in q'range(1) loop
          ip := ip + Conjugate(q(i,j))*q(i,k);
        end loop;
        sum := sum + AbsVal(ip);
      end loop;
    end loop;
    put("Orthogonality check : "); put(sum,3,3,3); new_line;
  end Orthogonality;

  procedure Test_QRD ( a,q,r : in Standard_Floating_Matrices.Matrix ) is

    wrk : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));
    use Standard_Floating_Matrices;

  begin
    put_line("The upper triangular part R :"); put(r,3);
    wrk := q*r;
    put_line("q*r :"); put(wrk,3); 
    put("Difference in 1-norm between the matrix and q*r: ");
    put(Differences(a,wrk),3,3,3); new_line;
    Orthogonality(q);
  end Test_QRD;

  procedure Test_QRD ( a,q,r : in Standard_Complex_Matrices.Matrix ) is

    wrk : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));
    use Standard_Complex_Matrices;

  begin
    put_line("The upper triangular part R :"); put(r,3);
    wrk := q*r;
    put_line("q*r :"); put(wrk,3); 
    put("Difference in 1-norm between the matrix and q*r : ");
    put(Differences(a,wrk),3,3,3); new_line;
    Orthogonality(q);
  end Test_QRD;

-- REAL TEST DRIVERS :

  procedure Real_LS_Test ( n,m : in natural; piv : in boolean;
                           a : in Standard_Floating_Matrices.Matrix;
                           b : in Standard_Floating_Vectors.Vector ) is

    wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
    qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
    jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
    sol : Standard_Floating_Vectors.Vector(1..m);
    rsd,dum : Standard_Floating_Vectors.Vector(1..n);
    info : integer;
    use Standard_Floating_Matrices;
    use Standard_Floating_Vectors;

  begin
    put_line("The matrix : "); put(a,3);
    QRD(wrk,qraux,jpvt,piv);
    put_line("The matrix after QR : "); put(wrk,3);
    put_line("The vector qraux : "); put(qraux,3); new_line;
    if piv
     then put("The vector jpvt : "); put(jpvt); new_line;
          Permute_Columns(wrk,jpvt);
    end if;
    QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
    if piv
     then Permute(sol,jpvt);
    end if;
    put_line("The solution : "); put(sol,3); new_line;
    put_line("The residual : "); put(rsd,3); new_line;
    dum := b - a*sol;
    put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
  end Real_LS_Test;          

  procedure Real_QR_Test ( n,m : in natural; piv : in boolean;
                           a : in Standard_Floating_Matrices.Matrix ) is

    wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
    bas : Standard_Floating_Matrices.Matrix(1..n,1..n);
    qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
    jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);

  begin
    put_line("The matrix : "); put(a,3);
    QRD(wrk,qraux,jpvt,piv);
    put_line("The matrix after QR : "); put(wrk,3);
    put_line("The vector qraux : "); put(qraux,3); new_line;
    if piv
     then put("The vector jpvt : "); put(jpvt); new_line;
          Permute_Columns(wrk,jpvt);
    end if;
    for i in wrk'range(1) loop
      for j in wrk'range(2) loop
        bas(i,j) := wrk(i,j);
      end loop;
      for j in n+1..m loop
        bas(i,j) := 0.0;
      end loop;
    end loop;
    Basis(bas,a);
    put_line("The orthogonal part Q of QR  :"); put(bas,3);
    Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
  end Real_QR_Test;

  procedure Interactive_Real_QR_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    put("Give a "); put(n,1); put("x"); put(m,1);   
    put_line(" matrix : "); get(a);
    Real_QR_Test(n,m,piv,a);
  end Interactive_Real_QR_Test;

  procedure Interactive_Real_LS_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Floating_Matrices.Matrix(1..n,1..m);
    b : Standard_Floating_Vectors.Vector(1..n);

  begin
    put("Give a "); put(n,1); put("x"); put(m,1);   
    put_line(" matrix : "); get(a);
    put("Give right-hand size "); put(n,1);
    put_line("-vector : "); get(b);
    Real_LS_Test(n,m,piv,a,b);
  end Interactive_Real_LS_Test;

  procedure Random_Real_QR_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);

  begin
    Real_QR_Test(n,m,piv,a);
  end Random_Real_QR_Test;

  procedure Random_Real_LS_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
    b : Standard_Floating_Vectors.Vector(1..n) := Random_Vector(n);

  begin
    Real_LS_Test(n,m,piv,a,b);
  end Random_Real_LS_Test;

-- COMPLEX TEST DRIVERS :

  procedure Complex_QR_Test ( n,m : in natural; piv : in boolean;
                              a : Standard_Complex_Matrices.Matrix ) is

    wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
    bas : Standard_Complex_Matrices.Matrix(1..n,1..n);
    qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
    jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);

  begin
    put_line("The matrix : "); put(a,3);
    QRD(wrk,qraux,jpvt,piv);
    put_line("The matrix after QR : "); put(wrk,3);
    put_line("The vector qraux : "); put(qraux,3); new_line;
   -- put("The vector jpvt : "); put(jpvt); new_line;
    if not piv
     then for i in wrk'range(1) loop
            for j in wrk'range(2) loop
              bas(i,j) := wrk(i,j);
            end loop;
            for j in n+1..m loop
              bas(i,j) := Create(0.0);
            end loop;
          end loop;
          Basis(bas,a);
          put_line("The orthogonal part Q of QR  :"); put(bas,3);
          Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
    end if;
  end Complex_QR_Test;

  procedure Complex_LS_Test ( n,m : in natural; piv : in boolean;
                              a : Standard_Complex_Matrices.Matrix;
                              b : Standard_Complex_Vectors.Vector ) is

    wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
    qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
    jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
    sol : Standard_Complex_Vectors.Vector(1..m);
    rsd,dum : Standard_Complex_Vectors.Vector(1..n);
    info : integer;
    use Standard_Complex_Matrices;
    use Standard_Complex_Vectors; 

  begin
    put_line("The matrix : "); put(a,3);
    QRD(wrk,qraux,jpvt,piv);
    put_line("The matrix after QR : "); put(wrk,3);
    put_line("The vector qraux : "); put(qraux,3); new_line;
   -- put("The vector jpvt : "); put(jpvt); new_line;
    QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
    put_line("The solution : "); put(sol,3); new_line;
    put_line("The residual : "); put(rsd,3); new_line;
    dum := b - a*sol;
    put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
  end Complex_LS_Test;

  procedure Interactive_Complex_QR_Test
                 ( n,m : in natural; piv : in boolean ) is

    a : Standard_Complex_Matrices.Matrix(1..n,1..m);

  begin
    put("Give a "); put(n,1); put("x"); put(m,1);
    put_line(" matrix : "); get(a);
    Complex_QR_Test(n,m,piv,a);
  end Interactive_Complex_QR_Test;

  procedure Interactive_Complex_LS_Test
                 ( n,m : in natural; piv : in boolean ) is

    a : Standard_Complex_Matrices.Matrix(1..n,1..m);
    b : Standard_Complex_Vectors.Vector(1..n);

  begin
    put("Give a "); put(n,1); put("x"); put(m,1);
    put_line(" matrix : "); get(a);
    put("Give right-hand size "); put(n,1);
    put_line("-vector : "); get(b); 
    Complex_LS_Test(n,m,piv,a,b);
  end Interactive_Complex_LS_Test;

  procedure Random_Complex_QR_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);

  begin
    Complex_QR_Test(n,m,piv,a);
  end Random_Complex_QR_Test;

  procedure Random_Complex_LS_Test ( n,m : in natural; piv : in boolean ) is

    a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
    b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(n);

  begin
    Complex_LS_Test(n,m,piv,a,b);
  end Random_Complex_LS_Test;

-- MAIN PROGRAM :

  procedure Main is

    n,m : natural;
    ans,choice : character;
    piv : boolean := false;

  begin
    new_line;
    put_line("Test on the QR decomposition and Least Squares.");
    loop
      new_line;
      put_line("Choose one of the following : ");
      put_line("  0. Exit this program.");
      put_line("  1. QR-decomposition on user-given floating-point matrix.");
      put_line("  2.                                complex matrix.");
      put_line("  3.                  on random floating-point matrix.");
      put_line("  4.                            complex matrix.");
      put_line("  5. Least Squares on user-given floating-point matrix.");
      put_line("  6.                             complex matrix.");
      put_line("  7.               on random floating-point matrix.");
      put_line("  8.                         complex matrix.");
	  put("Make your choice (0,1,2,3,4,5,6,7 or 8) : "); get(choice);
      exit when (choice = '0');
      put("Give the number of rows : "); get(n);
      put("Give the number of columns : "); get(m);
      put("Do you want pivoting ? (y/n) "); get(ans);
      piv := (ans = 'y');
      case choice is
        when '1' => Interactive_Real_QR_Test(n,m,piv);
        when '2' => Interactive_Complex_QR_Test(n,m,piv);
        when '3' => Random_Real_QR_Test(n,m,piv);
        when '4' => Random_Complex_QR_Test(n,m,piv);
        when '5' => Interactive_Real_LS_Test(n,m,piv);
        when '6' => Interactive_Complex_LS_Test(n,m,piv);
        when '7' => Random_Real_LS_Test(n,m,piv);
        when '8' => Random_Complex_LS_Test(n,m,piv);
        when others => null;
      end case;
    end loop;
  end Main;

begin
  Main;
end ts_qrd;