-- rmatin.adb  from 535 general eigenvalue, eigenvector fortran
with Ada.Numerics.Long_Complex_Types;
use  Ada.Numerics.Long_Complex_Types;
with complex_arrays; -- mine
use  complex_arrays; -- mine
with real_arrays;
use  real_arrays;
with ada.text_io;
use  ada.text_io;


procedure rmatin(Input:File_Type; -- open(input, in_file, file_name);
                 n:in out integer; a:in out complex_matrix) is

--  this input subroutine reads two real matrices  a  and  ai
--  from sysin of order n. n is read from file preceeding each matrix.
--  a=(ar,ai) contains a complex general matrix,
  seed : integer := 13; -- should be in a package
  val : real;
  M : Integer;
  Line : String(1..4096);
  Pos, Len : Integer;

  package Real_Io is new Float_Io(Real);
  use Real_Io;
  package Int_Io is new Integer_Io(Integer);
  use Int_Io;

  function drand return real is
    xhi, xalo, leftlo, fhi, k : integer;
    a:integer := 16807;
    b15:integer := 32768;
    b16:integer := 65536;
    p:integer := 2147483647;
  begin
      xhi := seed/b16;
      xalo := (seed-xhi*b16)*a;
      leftlo := xalo/b16;
      fhi := xhi*a+leftlo;
      k := fhi/b15;
      seed := (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k;
      if seed < 0 then
         seed := seed + p;
      end if;
      return real(seed)*4.656612875e-10;
  end drand;


  procedure mkrand(n:integer; iq:integer; a:in out complex_matrix) is
  begin
    put_line("mkrand n="&integer'image(n)&", seed="&integer'image(seed));
    for i in 1..n loop
      for j in 1..n loop
         if iq=1 then
            set_re(a(i,j),drand);
         else
            set_im(a(i,j),drand);
         end if;
      end loop;
    end loop;
  end mkrand;

  procedure  mkident(n:integer; iq:integer; a:out complex_matrix) is
  begin
    put_line("mkident n="&integer'image(n));
    for i in 1..n loop
      for j in 1..n loop
        if iq=1 then
          set_re(a(i,j),0.0);
          if i=j then
             set_re(a(i,j),1.0);
          end if;
        else
          set_im(a(i,j),0.0);
        end if;
      end loop;
    end loop;
  end mkident;

begin
  -- real part
  Get_Line(Input, Line, Len);
  Put_Line(Line(1..Len-1));
  get(Line(1..Len), N, Pos);
  -- Put_Line("n="&Integer'Image(N));

  if  n=0 then
    put_line("end of data for subroutine  rmatin");
    return; -- caller must test
  end if;
  get(Line(Pos+1..Len), M, pos);
  -- Put_Line("m="&Integer'Image(m));
  if m=-1 then
    mkrand(n,1,a);
    return;
  end if;
  if m=-2 then
    mkident(n,1,a);
    return;
  end if;
  -- read data
  for i in 1..n loop
    Get_Line(Input, Line, Len);
    -- Put_Line(Line(1..Len-1));
    Pos := 0;
    for j in 1..n loop
      get(Line(Pos+1..len), val, Pos);
      set_re(a(i,j),val);
    end loop;
  end loop;
  -- imaginary part
  Get_Line(Input, Line, Len);
  -- Put_Line(Line(1..Len-1));
  get(Line(1..Len), N, Pos);
  -- Put_Line("n="&Integer'Image(N));
  if  n=0 then
    put_line("end of data for subroutine  rmatin");
    return; -- caller must test
  end if;
  get(Line(Pos+1..Len), M, pos);
  -- Put_Line("m="&Integer'Image(m));
  if m=-1 then
    mkrand(n,2,a);
    return;
  end if;
  if m=-2 then
    mkident(n,2,a);
    return;
  end if;
  -- read data
  for i in 1..n loop
     Get_Line(Input, Line, Len);
     -- Put_Line(Line(1..Len-1));
     Pos := 0;
    for j in 1..n loop
      get(Line(Pos+1..len), val, Pos);
      set_im(a(i,j),val);
    end loop;
  end loop;
end rmatin;