-- test_535_eigen.c  based on 535 algorithm
--               check A*v = (lambda)B*v
--               Given n by n matrices A and B
--               use eigen.c to find vector of n eigenvalues
--               and matrix Z of corresponding eigenvalues
--               Check for residuals near zero.

with Ada.Numerics.Long_Complex_Types;
use  Ada.Numerics.Long_Complex_Types;
with Ada.Numerics.Generic_Elementary_Functions;
with Complex_Arrays; -- mine
use  Complex_Arrays; -- mine
with Real_Arrays;
use  Real_Arrays;
with Ada.Text_IO.Complex_IO; -- generic
-- with Complex_Arrays_IO; -- generic
with Ada.Text_IO;
use  Ada.Text_IO;
with Rmatin; -- from 535
with cxhes;  -- from 535
with cxval;  -- from 535
with cxvec;  -- from 535
with eigdet; -- mine
with cxinverse; -- mine
with evalvec; -- mine
with cxcheck; -- from 535
with Ada.Command_Line;

procedure test_535_eigen is
  n : Integer := 100; -- actual read in

  A, B :Complex_Matrix(1..N,1..N);  -- read file
  Z : Complex_matrix(1..n,1..n);    -- computed
  Alfa : Complex_Vector(1..n);
  Beta : Real_Vector(1..N);
  Norm : Real_Vector(1..N);
  Error : Integer := 0;
  Ahold, Bhold : Complex_Matrix(1..n,1..n);
  Resdul : Long_Float;
  Argc : Integer := Ada.Command_Line.Argument_Count;
  Fname : access String;
  Input : File_Type;

  package CX_IO is new Ada.Text_IO.Complex_IO(Ada.Numerics.Long_Complex_Types);
  use CX_IO;

  procedure Cxmprint(N:Integer; A : Complex_Matrix) is
  begin
    for I in 1..n loop
      for J in 1..n loop
        Put("("&Integer'Image(I)&","&Integer'Image(J)&")=");
        Put(A(I,J)); New_Line;
      end loop ;
    end loop ;
  end Cxmprint;

  procedure Cxmcopy(N:Integer; A : Complex_Matrix; B : out Complex_matrix) is
  begin
    for I in 1..n loop
      for J in 1..n loop
        B(I,J) := A(I,J);
      end loop ;
    end loop ;
  end Cxmcopy;

begin

  Put_Line("test_535_eigen.c running, check 535_double_f.out");
  for I in 1..Argc loop
    Ada.Text_IO.Put_Line(Ada.Command_Line.Argument(I));
  end loop;
  if Argc>0 then
    Fname := new String'(Ada.Command_Line.Argument(1));
  else
    Fname := new String'("535.dat");
  end if;
  Open(Input, In_File, Fname.all);
  Rmatin(Input, N, A); -- N read from file
  Put_Line("A matrix");
  cxmprint(N, A);
  cxmcopy(N, A,Ahold); -- restore A matrix
  Rmatin(Input, N, B);
  Put_Line("B matrix");
  cxmprint(N, B);
  cxmcopy(N, B,Bhold); -- restore B matrix

  -- eigen(n, A, B, Z, alfa, beta, error); all

  cxhes(n, A, B, Z);
  Put_Line("cxhes results A matrix");
  cxmprint(N, A);
  Put_Line("cxhes results B matrix");
  cxmprint(N, B);
  Put_Line("cxhes results Z matrix");
  cxmprint(N, Z);
  New_Line;
  New_Line;

  cxval(n, A, B, Z, Alfa, Beta, error);
  Put_Line("cxval results A matrix");
  cxmprint(N, A);
  Put_Line("cxval results B matrix");
  cxmprint(N, B);
  Put_Line("cxval results Z matrix");
  cxmprint(N, Z);
  New_Line;
  New_Line;
  Put_Line("cxval error="&Integer'Image(Error));
  Put_Line("              ALFR(I)              ALFI(I)               BETA(I)");
  for I in 1..N loop
    Put_Line(Integer'Image(I)&"  "&Long_Float'Image(Re(Alfa(I)))&
             "  "&Long_Float'Image(Im(Alfa(I)))&
             "  "&Long_Float'Image(Beta(I)));
  end loop;
  New_Line;
  cxvec(n, A, B, Z, Alfa, Beta);
  Put_Line("cxvec final eigenvctors");
  cxmprint(N, Z);
  New_Line;
  cxcheck(n, Ahold, Bhold, Z, Alfa, Beta, Norm, Resdul);
  Put_Line("cxcheck resdul="&Long_Float'Image(Resdul));
  New_Line;
  Put_Line("eigdet check eigenvalues");
  eigdet(N, Ahold, Alfa);
  New_Line;
  Put_Line("evalvec check eigenvectors");
  evalvec(N, Ahold, Bhold, Z, Alfa);
  New_Line;
  New_Line;
  Put_Line("test_535_eigen.adb finished");
end test_535_eigen;