with Ada.Text_IO;                       use Ada.Text_IO;

with Ada.Numerics.Generic_Real_Arrays;
with Generic_Real_Linear_Equations;

procedure Test_QR is

  subtype Real is Long_Float;
  type Integer_Vector is array (Integer range <>) of Integer;

  package RIO is new Float_IO(Real);

  package GRA is new Ada.Numerics.Generic_Real_Arrays(Real);
  package GRLE is new Generic_Real_Linear_Equations(Real, GRA, Integer_Vector);

  use GRA, GRLE, RIO;

  procedure Put(A: Real_Matrix) is
  begin
    for i in A'Range(1) loop
      for j in A'Range(2) loop
        Put(A(i,j), 4,2,0);
      end loop;
      New_Line;
    end loop;
  end Put;

  procedure Put(x: Real_Vector) is
  begin
    for i in x'Range loop
      Put(x(i), 4,2,0);
    end loop;
    New_Line;
  end Put;

  procedure Test(A: Real_Matrix; x, y: Real_Vector) is
    Q: Real_Matrix(A'Range(1), A'Range(1));
    R: Real_Matrix(A'Range(2), A'Range(2));
    Z: Real_Matrix:= A;
    res: Real:= 0.0;
    x_solved: Real_Vector(x'Range);
  begin
    QR_Decomposition ( A , Q , R ) ;
    Put_Line("**** QR test A -> QR. A is:");
    Put(A);
    Put_Line("Q is:");
    Put(Q);
    Put_Line("R is:");
    Put(R);
    Put_Line("Check Q orthogonal: QtQ should be close to unit matrix");
    Put(Transpose(Q) * Q);
    Put_Line("Check QR - A, should be close to 0 matrix");
    Z:= Q*R-A;
    Put(Z);
    for i in Z'Range(1) loop
      for j in Z'Range(2) loop
        res:= res + abs Z(i,j);
      end loop;
    end loop;
    Put("This sum should be close to 0 :");
    Put(res);
    New_Line;
    Put_Line("**** QR solver");
    x_solved:= QR_Solve(Q,R,y);
    Put_Line("  solution x  :");
    Put(x_solved);
    Put_Line("  x should be :");
    Put(x);
    Put("  Ax=");
    Put(A*x_solved);
    Put("   y=");
    Put(y);
    New_Line;
    New_Line;
  end Test;

  A33: constant Real_Matrix(1..3, 1..3):=
        (( 12.0 , -51.0,   4.0),
         (  6.0 , 167.0, -68.0),
         ( -4.0,   24.0, -41.0));

  x_a: constant Real_Vector:= (1.0, 2.0, 3.0);
  y_a: constant Real_Vector:= A33 * x_a;

  --  A32: constant Real_Matrix(1..3, 1..2):=
  --        ((  2.0, 1.0),
  --         (  4.0, 5.0),
  --         (  3.0, 7.0));

  --  x_b: constant Real_Vector:= (123.0, 456.0);
  --  y_b: constant Real_Vector:= A32 * x_b;

begin
  Test( A33, x_a, y_a );
  -- Test( A32, x_b, y_b ); -- !! QR not yet set for non-square matrices
  Put("Finished - press return");
  Skip_Line;
end Test_QR;