with Ada.Integer_Text_IO, Globals, Strings, Protected_Counters;
use Globals, Strings, Protected_Counters;

pragma Elaborate_All (Ada.Integer_Text_IO,Globals,Strings,Protected_Counters);

package body Matrices is

  package Txt_IO renames Ada.Text_IO;
  package Int_IO renames Ada.Integer_Text_IO;

  function BothFirst(A: Matrix; First: Integer) return Boolean is
  begin
    return (A'First(1)=First) and (A'First(2)=First);
  end BothFirst;

  function BothLast(A: Matrix; Last: Integer) return Boolean is
  begin
    return (A'Last(1)=Last) and (A'Last(2)=Last);
  end BothLast;

  function IsSquare(A: Matrix) return Boolean is
  begin
    return (A'First(1)=A'First(2)) and (A'Last(1)=A'Last(2));
  end IsSquare;

  function SameDim(A1,A2: Matrix) return Boolean is
  begin
    return (A1'First(1)=A2'First(1)) and (A1'Last(1)=A2'Last(1))
      and (A1'First(2)=A2'First(2)) and (A1'Last(2)=A2'Last(2));
  end SameDim;

  function TransDim(A1,A2: Matrix) return Boolean is
  begin
    return (A1'First(1)=A2'First(2)) and (A1'Last(1)=A2'Last(2))
      and (A1'First(2)=A2'First(1)) and (A1'Last(2)=A2'Last(1));
  end TransDim;

  function EffLast1(A: Matrix) return Integer is
  begin
    for I in reverse A'First(1)+1 .. A'Last(1) loop
      for J in A'Range(2) loop
        if not IsZero(A(I,J)) then return I; end if;
      end loop;
    end loop;
    return A'First(1);
  end EffLast1;

  function EffLast2(A: Matrix) return Integer is
  begin
    for J in reverse A'First(2)+1 .. A'Last(2) loop
      for I in A'Range(1) loop
        if not IsZero(A(I,J)) then return J; end if;
      end loop;
    end loop;
    return A'First(2);
  end EffLast2;

  --------------------------------

  function "="(A1,A2: Matrix) return Boolean is
  begin
    if not SameDim(A1,A2) then return False; end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        if not (A1(I,J)=A2(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end "=";

  function ZeroMatrix(F,L: Integer) return Matrix is
    A: Matrix(F..L,F..L);
  begin
    SetZero(A);
    return A;
  end ZeroMatrix;

  procedure ZeroRow(I: in Integer; A: in out Matrix) is
  begin
    for J in A'Range(2) loop
      SetZero(A(I,J));
    end loop;
  end ZeroRow;

  procedure ZeroColumn(J: in Integer; A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      SetZero(A(I,J));
    end loop;
  end ZeroColumn;

  procedure Assign(K: in Integer; A: in out Matrix) is
  begin
    if CheckRange and then not IsSquare(A) then
      raise Dimension_Error with "Matrices.Assign error";
    end if;
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        SetZero(A(I,J));
      end loop;
      Assign(K,A(I,I));
    end loop;
  end Assign;

  procedure Identity(A: in out Matrix) is
  begin
    Assign(1,A);
  end Identity;

  function Identity(L: Positive) return Matrix is
    A: Matrix(1..L,1..L);
  begin
    Assign(1,A);
    return A;
  end Identity;

  procedure SwapRows(I1,I2: in Integer; A: in out Matrix) is
  begin
    if I1 /= I2 then
      for J in A'Range(2) loop
        Swap(A(I1,J),A(I2,J));
      end loop;
    end if;
  end SwapRows;

  procedure SwapColumns(J1,J2: in Integer; A: in out Matrix) is
  begin
    if J1 /= J2 then
      for I in A'Range(1) loop
        Swap(A(I,J1),A(I,J2));
      end loop;
    end if;
  end SwapColumns;

  procedure Transpose(A: in out Matrix) is
  begin
    if CheckRange and then not IsSquare(A) then
      raise Dimension_Error with "Matrices.Transpose error";
    end if;
    for I in A'Range(1) loop
      for J in A'First .. I-1 loop
        Swap(A(I,J),A(J,I));
      end loop;
    end loop;
  end Transpose;

  procedure Transpose(A1: in Matrix; A2: in out Matrix) is
  begin
    if CheckRange and then not TransDim(A1,A2) then
      raise Dimension_Error with "Matrices.Transpose error";
    end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        Copy(A1(I,J),A2(J,I));
      end loop;
    end loop;
  end Transpose;

  function Transpose(A: Matrix) return Matrix is
    B: Matrix(A'Range(2),A'Range(1));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Copy(A(I,J),B(J,I));
      end loop;
    end loop;
    return B;
  end Transpose;

  procedure Adjoint(A1: in Matrix; A2: in out Matrix) is
  begin
    if CheckRange and then not TransDim(A1,A2) then
      raise Dimension_Error with "Matrices.Adjoint error";
    end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        Adjoint(A1(I,J),A2(J,I));
      end loop;
    end loop;
  end Adjoint;

  procedure Adjoint(A: in out Matrix) is
  begin
    Transpose(A);
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Adjoint(A(I,J));
      end loop;
    end loop;
  end Adjoint;

  function Adjoint(A: Matrix) return Matrix is
    B: Matrix(A'Range(2),A'Range(1));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Adjoint(A(I,J),B(J,I));
      end loop;
    end loop;
    return B;
  end Adjoint;

  function Trace(A: Matrix) return Scalar is
    F: constant Integer := A'First(1);
    S: Scalar;
  begin
    if CheckRange and then not IsSquare(A) then
      raise Dimension_Error with "Matrices.Trace error";
    end if;
    Copy(A(F,F),S);
    for I in F+1 .. A'Last(1) loop
      Add(A(I,I),S);
    end loop;
    return S;
  end Trace;

  procedure AddConst(I: in Integer; A: in out Matrix) is
    S: Scalar;
  begin
    if CheckRange and then not IsSquare(A) then
      raise Dimension_Error with "Matrices.AddConst error";
    end if;
    Assign(I,S);
    for K in A'Range(1) loop
      Add(S,A(K,K));
    end loop;
  end AddConst;

  procedure AddConst(S: in Scalar; A: in out Matrix) is
  begin
    if CheckRange and then not IsSquare(A) then
      raise Dimension_Error with "Matrices.AddConst error";
    end if;
    for K in A'Range(1) loop
      Add(S,A(K,K));
    end loop;
  end AddConst;

  procedure Mult(S: in Scalar; A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Mult(S,A(I,J));
      end loop;
    end loop;
  end Mult;

  procedure Prod(S: in Scalar; A1: in Matrix; A2: in out Matrix) is
  begin
    if CheckRange and then not SameDim(A1,A2) then
      raise Dimension_Error with "Matrices.Prod error";
    end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        Prod(S,A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Prod;

  function "*"(S: Scalar; A: Matrix) return Matrix is
    B: Matrix(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Prod(S,A(I,J),B(I,J));
      end loop;
    end loop;
    return B;
  end "*";

  procedure AddProd(S: in Scalar; A1: in Matrix; A2: in out Matrix) is
    Tmp: Scalar;
  begin
    if CheckRange and then not SameDim(A1,A2) then
      raise Dimension_Error with "Matrices.AddProd error";
    end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        AddProd(S,A1(I,J),A2(I,J),Tmp);
      end loop;
    end loop;
  end AddProd;

  ----------------------

  procedure Mult(A1: in Matrix; A2,Tmp: in out Matrix) is
  begin
    Copy(A2,Tmp);
    SetZero(A2);
    AddProd(A1,Tmp,A2);
  end Mult;

  procedure Mult(A1: in Matrix; A2: in out Matrix) is
    Tmp: Matrix(A2'Range(1),A2'Range(2));
  begin
    Copy(A2,Tmp);
    SetZero(A2);
    AddProd(A1,Tmp,A2);
  end Mult;

  procedure Prod(A1,A2: in Matrix; A3: in out Matrix) is
  begin
    SetZero(A3);
    AddProd(A1,A2,A3);
  end Prod;

  function "*"(A1,A2: Matrix) return Matrix is
    A3: Matrix(A1'Range(1),A2'Range(2));
  begin
    SetZero(A3);
    AddProd(A1,A2,A3);
    return A3;
  end "*";

  procedure AddProd(A1,A2: in Matrix; A3: in out Matrix) is
    F:  constant Integer := A3'First(1);
    L:  constant Integer := A3'Last(1);
    Task_Error: Boolean := False;
    PC: Protected_Counter(F);

    procedure Prod_Loop is
      I: Integer;
      Tmp: Scalar;
    begin
      Proper_Rounding;
      loop
        PC.Next(I);
        exit when I>L;
        for J in A3'Range(2) loop
          declare
            A3IJ: Scalar renames A3(I,J);
          begin
            for K in A1'Range(2) loop
              AddProd(A1(I,K),A2(K,J),A3IJ,Tmp);
            end loop;
          end;
        end loop;
      end loop;
    exception
      when others => Task_Error := True; raise;
    end Prod_Loop;

    NT: constant Natural := Reserve_Tasks(L-F,MProd_Parallel);
  begin
    if CheckRange and then ((A1'First(1) /= F) or (A1'Last(1) /= L)
                              or (A1'First(2) /= A2'First(1)) or (A1'Last(2) /= A2'Last(1))
                              or (A2'First(2) /= A3'First(2)) or (A2'Last(2) /= A3'Last(2))) then
      raise Dimension_Error with "Matrices.AddProd error";
    end if;
    if NT=0 then
      Prod_Loop;
    else
      declare
        task type Prod_Task_Type is end Prod_Task_Type;
        task body Prod_Task_Type is begin Prod_Loop; end Prod_Task_Type;
        Prod_Task: array(1 .. NT) of Prod_Task_Type;
        pragma Warnings (Off,Prod_Task);
      begin
        Prod_Loop;
      end;
      Free_Tasks(NT);
      if Task_Error then raise Sorry with "Task_Error"; end if;
    end if;
  end AddProd;

  function Commute4Sure(A1,A2: Matrix) return Boolean is
    F: constant Integer := A1'First(1);
    L: constant Integer := A1'Last(1);
    S,T: Scalar;
  begin
    if CheckRange and then not (IsSquare(A1) and SameDim(A1,A2)) then
      raise Dimension_Error with "Matrices.Commute4Sure error";
    end if;
    SetZero(S);
    for I in F .. L loop
      for J in F .. L loop
        for K in F .. L loop
          Prod(A1(I,K),A2(K,J),T);
          Add(T,S);
          Prod(A2(I,K),A1(K,J),T);
          Sub(T,S);
        end loop;
        if not IsZero(S) then return False; end if;
      end loop;
    end loop;
    return True;
  end Commute4Sure;

  ---------------------

  procedure Show1(N: in String; A: in Matrix; Hide0: in Boolean := True) is
    W: constant Positive := 4;
    AllZero: Boolean := True;
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        if not IsZero(A(I,J)) then
          Show1(N & Strng(I,W) & Strng(J,W) & " ",A(I,J));
          AllZero := False;
        end if;
      end loop;
    end loop;
    if not Hide0 then
      if AllZero then Show0(N & "0"); end if;
      Show0;
    end if;
  end Show1;

  procedure Show2(N: in String; A1,A2: in Matrix; Hide0: in Boolean := True) is
    W: constant Positive := 4;
    AllZero: Boolean := True;
  begin
    if CheckRange and then not SameDim(A1,A2) then
      raise Dimension_Error with "Matrices.Show2 error";
    end if;
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        if not (IsZero(A1(I,J)) and then IsZero(A2(I,J))) then
          Show2(N & Strng(I,W) & Strng(J,W) & " ",A1(I,J),A2(I,J));
          AllZero := False;
        end if;
      end loop;
    end loop;
    if not Hide0 then
      if AllZero then Show0(N & "0 0"); end if;
      Show0;
    end if;
  end Show2;

  procedure PutDim(F: in File_Type; L1,L2: in Positive) is
    W1: constant Positive := Width(L1)+1;
    W2: constant Positive := Width(L2)+1;
  begin
    Int_IO.Put(F,L1,W1);
    New_Line(F);
    Int_IO.Put(F,L2,W2);
    New_Line(F);
  end PutDim;

  procedure PutCoeffs(F: in File_Type; A: in Matrix; Decimal: in Boolean := False) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Put(F,A(I,J),Decimal);
      end loop;
    end loop;
  end PutCoeffs;

  procedure Put(F: in File_Type; A: in Matrix; Decimal: in Boolean := False) is
  begin
    PutDim(F,A'Length(1),A'Length(2));
    PutCoeffs(F,A,Decimal);
  end Put;

  function Get(F: File_Type; AFirst1,AFirst2: Integer; Decimal: Boolean := False) return Matrix is
    L1,L2: Integer;
  begin
    Int_IO.Get(F,L1);
    Int_IO.Get(F,L2);
    declare
      subtype ARange1 is Integer range AFirst1 .. AFirst1+L1-1;
      subtype ARange2 is Integer range AFirst2 .. AFirst2+L2-1;
      A: Matrix(ARange1,ARange2);
    begin
      for M in ARange1 loop
        for N in ARange2 loop
          Get(F,A(M,N),Decimal);
        end loop;
      end loop;
      return A;
    end;
  end Get;

  procedure Get(F: in File_Type; A: in out Matrix; Decimal: in Boolean := False) is
    Trunc: Boolean := False;
    L1,L2: Integer;
    S: Scalar;
  begin
    SetZero(A);
    Int_IO.Get(F,L1);
    Int_IO.Get(F,L2);
    declare
      subtype ARange1 is Integer range A'First(1) .. A'First(1)+L1-1;
      subtype ARange2 is Integer range A'First(2) .. A'First(2)+L2-1;
    begin
      for M in ARange1 loop
        for N in ARange2 loop
          if M <= A'Last(1) and then N <= A'Last(2) then
            Get(F,A(M,N),Decimal);
          else
            Get(F,S,Decimal);
            Trunc := Trunc or else not IsZero(S);
          end if;
        end loop;
      end loop;
    end;
    if Trunc then
      Show0("Matrices.Get Warning: input matrix truncated");
    end if;
  end Get;

  procedure Write(FileName: in String; A: in Matrix; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Writing " & FileName); end if;
    Create(F,Out_File,FileName);
    Put(F,A,Decimal);
    Close(F);
  end Write;

  function Read(FileName: String; AFirst1,AFirst2: Integer; Decimal: Boolean := False) return Matrix is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Reading " & FileName); end if;
    Open(F,In_File,FileName);
    declare
      A: constant Matrix := Get(F,AFirst1,AFirst2,Decimal);
    begin
      Close(F);
      return A;
    end;
  end Read;

  procedure Read(FileName: in String; A: in out Matrix; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Reading " & FileName); end if;
    Open(F,In_File,FileName);
    Get(F,A,Decimal);
    Close(F);
  end Read;

  procedure Append(FileName: in String; A: in Matrix; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Appending " & FileName); end if;
    begin
      Open(F,Append_File,FileName);
    exception
      when Name_Error => Create(F,Out_File,FileName);
    end;
    Put(F,A,Decimal);
    Close(F);
  end Append;

  procedure WriteProd(A1,A2: in Matrix; FileName: in String; Decimal: Boolean := False) is
    S,Tmp: Scalar;
    F: File_Type;
  begin
    if Verbosity>0 then Txt_IO.Put_Line("Writing " & FileName); end if;
    Create(F,Out_File,FileName);
    PutDim(F,A1'Length(1),A2'Length(2));
    for I in A1'Range(1) loop
      for J in A2'Range(2) loop
        SetZero(S);
        for K in A1'Range(2) loop
          AddProd(A1(I,K),A2(K,J),S,Tmp);
        end loop;
        Put(F,S,Decimal);
      end loop;
    end loop;
    Close(F);
  end WriteProd;

end Matrices;
