with Ada.Text_IO, MGeneric, Integer_Subtypes;
use Ada.Text_IO, MGeneric, Integer_Subtypes;

pragma Elaborate_All (Ada.Text_IO,MGeneric,Integer_Subtypes);

generic

  type Scalar is private;

  with function "="(S1,S2: Scalar) return Boolean is <>;
  with function IsZero(S: Scalar) return Boolean is  <>;
  with function IsReal(S: Scalar) return Boolean is  <>;
  with procedure SetZero(S: in out Scalar) is <>;
  with procedure Assign(I: in Integer; S: in out Scalar) is <>;
  with procedure Copy(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Swap(S1,S2: in out Scalar) is <>;
  with procedure Neg(S: in out Scalar) is <>;
  with procedure Neg(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Add(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Sum(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure Sub(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Diff(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure Mult(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Prod(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure AddProd(S1,S2: in Scalar; S3,Tmp: in out Scalar) is <>;
  with procedure Real_Part(S: in out Scalar) is <>;
  with procedure Imag_Part(S: in out Scalar) is <>;
  with procedure Adjoint(S: in out Scalar) is <>;
  with procedure Adjoint(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Show1(N: in String; S: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Show2(N: in String; S1,S2: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Put(F: in File_Type; S: in Scalar; Decimal: in Boolean := False) is <>;
  with procedure Get(F: in File_Type; S: in out Scalar; Decimal: in Boolean := False) is <>;
  with procedure Proper_Rounding is <>;

package Matrices is

  type Matrix is array(Integer range <>, Integer range <>) of Scalar;

  --- basic
  function "="(A1,A2: Matrix) return Boolean;
  function IsZero is new MF_BProperty (Scalar,Matrix,IsZero);
  function IsReal is new MF_BProperty (Scalar,Matrix,IsReal);
  procedure SetZero is new MP_Modify1 (Scalar,Matrix,SetZero);
  function ZeroMatrix(F,L: Integer) return Matrix;
  procedure ZeroRow(I: in Integer; A: in out Matrix);
  procedure ZeroColumn(J: in Integer; A: in out Matrix);
  procedure Assign(K: in Integer; A: in out Matrix);
  procedure Identity(A: in out Matrix);
  function Identity(L: Positive) return Matrix;
  procedure Copy is new MP_Map1 (Scalar,Matrix,Copy);
  procedure Swap is new MP_Modify2 (Scalar,Matrix,Swap);
  procedure SwapRows(I1,I2: in Integer; A: in out Matrix);
  procedure SwapColumns(J1,J2: in Integer; A: in out Matrix);
  procedure Transpose(A: in out Matrix);
  procedure Transpose(A1: in Matrix; A2: in out Matrix);
  function Transpose(A: Matrix) return Matrix;
  procedure Real_Part is new MP_Modify1 (Scalar,Matrix,Real_Part);
  procedure Imag_Part is new MP_Modify1 (Scalar,Matrix,Imag_Part);
  procedure Adjoint(A: in out Matrix);
  procedure Adjoint(A1: in Matrix; A2: in out Matrix);
  function Adjoint(A: Matrix) return Matrix;
  function Trace(A: Matrix) return Scalar;
  procedure AddConst(I: in Integer; A: in out Matrix);
  procedure AddConst(S: in Scalar; A: in out Matrix);
  procedure Neg is new MP_Modify1 (Scalar,Matrix,Neg);
  procedure Neg is new MP_Map1 (Scalar,Matrix,Neg);
  function "-" is new MF_Map1 (Scalar,Matrix,Neg);
  procedure Add is new MP_Map1 (Scalar,Matrix,Add);
  procedure Sum is new MP_Map2 (Scalar,Matrix,Sum);
  function "+" is new MF_Map2 (Scalar,Matrix,Sum);
  procedure Sub is new MP_Map1 (Scalar,Matrix,Sub);
  procedure Diff is new MP_Map2 (Scalar,Matrix,Diff);
  function "-" is new MF_Map2 (Scalar,Matrix,Diff);
  procedure Mult(S: in Scalar; A: in out Matrix);
  procedure Prod(S: in Scalar; A1: in Matrix; A2: in out Matrix);
  function "*"(S: Scalar; A: Matrix) return Matrix;
  procedure AddProd(S: in Scalar; A1: in Matrix; A2: in out Matrix);
  procedure Mult(A1: in Matrix; A2,Tmp: in out Matrix);
  procedure Mult(A1: in Matrix; A2: in out Matrix);
  procedure Prod(A1,A2: in Matrix; A3: in out Matrix);
  function "*"(A1,A2: Matrix) return Matrix;
  procedure AddProd(A1,A2: in Matrix; A3: in out Matrix);
  function Commute4Sure(A1,A2: Matrix) return Boolean;

  --- i/o and conversion
  procedure Show1(N: in String; A: in Matrix; Hide0: in Boolean := True);
  procedure Show2(N: in String; A1,A2: in Matrix; Hide0: in Boolean := True);
  procedure PutDim(F: in File_Type; L1,L2: in Positive);
  procedure PutCoeffs(F: in File_Type; A: in Matrix; Decimal: in Boolean := False);
  procedure Put(F: in File_Type; A: in Matrix; Decimal: in Boolean := False);
  function Get(F: File_Type; AFirst1,AFirst2: Integer; Decimal: Boolean := False) return Matrix;
  procedure Get(F: in File_Type; A: in out Matrix; Decimal: in Boolean := False);
  procedure Write(FileName: in String; A: in Matrix; Decimal: in Boolean := False);
  function Read(FileName: String; AFirst1,AFirst2: Integer; Decimal: Boolean := False) return Matrix;
  procedure Read(FileName: in String; A: in out Matrix; Decimal: in Boolean := False);
  procedure Append(FileName: in String; A: in Matrix; Decimal: in Boolean := False);
  procedure WriteProd(A1,A2: in Matrix; FileName: in String; Decimal: Boolean := False);

  --- auxiliary
  function BothFirst(A: Matrix; First: Integer) return Boolean;
  function BothLast(A: Matrix; Last: Integer) return Boolean;
  function IsSquare(A: Matrix) return Boolean;
  function SameDim(A1,A2: Matrix) return Boolean;     --- A1,A2 same dimensions
  function TransDim(A1,A2: Matrix) return Boolean;    --- A1,Transpose(A2) same dimensions
  function EffLast1(A: Matrix) return Integer;
  function EffLast2(A: Matrix) return Integer;

  --- baby cases

  subtype Matrix2 is Matrix(OneTwo,OneTwo);           --- 2x2
  subtype Matrix3 is Matrix(OneTwoThree,OneTwoThree); --- 3x3

  pragma Inline_Always(IsSquare,SameDim,TransDim);

  MProd_Parallel: Boolean := False;

  type Triangular is array(Natural range <>, Natural range <>) of Scalar;

end Matrices;
