with Ada.Text_IO, Scalars, Logicals, Ints, Flts, Flt_Vectors;
use Ada.Text_IO, Scalars, Logicals, Ints, Flts;
with Matrices, Vectors, Vectors.Sorting, Vectors.Ops;

pragma Elaborate_All (Ada.Text_IO,Scalars,Logicals,Ints,Flts,Flt_Vectors);
pragma Elaborate_All (Matrices,Vectors.Sorting,Vectors.Ops);

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 AddProd(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure SubProd(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 <>;

  with package SM is new Matrices (Scalar => Scalar);
  with package SV is new Vectors (Component => Scalar);

  with function Info(Dummy: Scalar) return Scalar_Info is <>;
  with procedure Center(S: in out Scalar) is <>;
  with procedure BallAt0(R: in Flt; S: in out Scalar) is <>;
  with procedure Assign(R: in Flt; S: in out Scalar) is <>;
  with function "-"(S: Scalar) return Scalar is <>;
  with procedure Prod(R: in Flt; S1: in Scalar; S2: in out Scalar) is <>;
  with function "*"(R: Flt; S: Scalar) return Scalar is <>;
  with function "*"(S1,S2: Scalar) return Scalar is <>;
  with procedure Inv(S1: in Scalar; S2: in out Scalar) is <>;
  with function Inv(S: Scalar) return Scalar is <>;
  with procedure Quot(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with function Exp(S: Scalar) return Scalar is <>;
  with function Log(S: Scalar) return Scalar is <>;
  with procedure Max(S1: in Scalar; S2: in out Scalar) is <>;
  with function Inf(S: Scalar) return Flt is <>;
  with function Sup(S: Scalar) return Flt is <>;
  with procedure Norm(S1: in Scalar; S2: in out Scalar) is <>;
  with function MaxNorm(S: Scalar) return Radius is <>;
  with procedure Sqrt(S1: in Scalar; S2: in out Scalar) is <>;
  with function Sqrt(S: Scalar) return Scalar is <>;
  with function Scal(R: Flt) return Scalar is <>;
  with function Approx(S: Scalar) return Flt is <>;
  with procedure Mult(R: in Flt; S: in out Scalar) is <>;

package ScalVectors is

  package SV_Ops is new SV.Ops (Scalar => Scalar, SM => SM);
  package SV_Sorting is new SV.Sorting;

  subtype Vector is SV.Vector;
  subtype Matrix is SM.Matrix;
  subtype Pivot is Int_Vec;
  subtype Permutation is Int_Vec;
  subtype FltVec is Flt_Vectors.Vector;
  type Diag is new Vector;

  --- diagonal matrices
  function GetDiag(A: in Matrix) return Diag;
  procedure Add(D: in Diag; A: in out Matrix);                     -- A := A+D
  procedure Sub(D: in Diag; A: in out Matrix);                     -- A := A-D
  procedure Mult(D: in Diag; W: in out Vector);                    -- W := D*W
  procedure Mult(D1: in Diag; D2: in out Diag);                    -- D2 := D1*D2
  procedure Mult(D: in Diag; A: in out Matrix);                    -- A := D*A
  procedure Mult(A: in out Matrix; D: in Diag);                    -- A := A*D
  function "*"(D: Diag; A: Matrix) return Matrix;                  -- return D*A
  function "*"(A: Matrix; D: Diag) return Matrix;                  -- return A*D
  procedure Invert(D: in out Diag);                                -- invert D
  function Inverse(D: in Diag) return Diag;                        -- return D^{-1}

  --- norms
  procedure Norm1(V: in Vector; S: in out Scalar);                --- sum
  function Norm1(V: Vector) return Scalar;                        --- sum
  function MaxNorm1(V: Vector) return Radius;
  procedure WNorm1(W,V: in Vector; S: in out Scalar);             --- weighted sum
  procedure WNorm1(W: in RadVec; V: in Vector; S: in out Scalar); --- weighted sum
  function WNorm1(W: RadVec; V: Vector) return Scalar;            --- weighted sum
  procedure WBallAt0(W: in RadVec; R: in Flt; V: in out Vector);
  procedure Norm1(A: in Matrix; S: in out Scalar);                --- sum
  function Norm1(A: Matrix) return Scalar;                        --- sum
  function MaxNorm1(A: Matrix) return Radius;                     --- upper bound on Norm1(A)
  function CheckNorm1(A: Matrix; R: Radius) return Boolean;       --- check if MaxNorm1(A)<=R
  procedure WNorm1(W: in Vector; A: in Matrix; S: in out Scalar); --- weighted sum
  procedure WNorm1(W: in RadVec; A: in Matrix; S: in out Scalar); --- weighted sum
  function WNorm1(W: RadVec; A: Matrix) return Scalar;            --- weighted sum
  procedure WBallAt0(W: in RadVec; R: in Flt; A: in out Matrix);
  function NormI(V: Vector) return Scalar;                        --- max
  function Norm2(V: Vector) return Scalar;
  function MaxNorm2(A: Matrix; Pow4: Positive := 4) return Radius;
  function CheckNorm2(A: Matrix; R: Radius; Pow4: Positive := 4) return Boolean; --- check if MaxNorm2(A)<=R

  --- misc
  procedure Center(V: in out Vector);
  procedure Center(A: in out Matrix);
  procedure Dot(V1,V2: in Vector; S: in out Scalar);
  function Dot(V1,V2: Vector) return Scalar;         --- linear in both args
  function "*"(V1,V2: Vector) return Scalar;         --- antilinear in V1
  function "*"(S: Scalar; U: FltVec) return Vector;
  procedure Truncate(Eps: in Radius; V: in out Vector);
  procedure Truncate(Eps: in Radius; A: in out Matrix);
  procedure SetRow(I: in Integer; V: in Vector; A: in out Matrix);
  procedure SetColumn(J: in Integer; V: in Vector; A: in out Matrix);
  procedure GetRow(I: in Integer; A: in Matrix; V: in out Vector);
  function GetRow(I: Integer; A: Matrix) return Vector;
  procedure GetColumn(J: in Integer; A: in Matrix; V: in out Vector);
  function GetColumn(J: Integer; A: Matrix) return Vector;
  procedure Permute(P: in Permutation; V,Tmp: in out Vector);
  procedure Row_Permute(P: in Permutation; A: in out Matrix; Tmp: in out Vector);
  procedure Column_Permute(P: in Permutation; A: in out Matrix; Tmp: in out Vector);
  procedure Symmetrize(A: in out Matrix);
  function AntiSymm(A: in Matrix) return Matrix;

  procedure Orthogonalize(A: in out Matrix);                           -- Gram-Schmidt
  function Determinant(A: Matrix) return Scalar;
  procedure Solve(B: in Vector; A: in Matrix; X: in out Vector);       -- Solve A*X=B
  procedure Solve(A: in out Matrix; X: in out Vector);                 -- Solve A*X_out=X_in (changes A)
  procedure Old_Invert(A: in out Matrix);                                  -- Invert A
  procedure InvNewton(A: in Matrix; Ai: in out Matrix; Iter: in Positive);
  procedure InvNewton(A: in out Matrix; Iter: in Positive);            -- invert and then apply N Newton steps
  procedure Serial_Invert(A: in out Matrix) renames Old_Invert;
  procedure Parallel_Invert(NT,Last1: in Positive; A: in out Matrix);
  procedure Invert(A: in out Matrix);

  MInvert_Parallel: Boolean := True;

  --- QR stuff
  procedure QR(A,Q: in out Matrix);                                       -- A -> QA by Householder reflections
  procedure Francis_Upper(A,Q: in out Matrix; Symm: in Boolean := False); -- A -> Q^tAQ by Householder reflections
  procedure Francis_Lower(A,Q: in out Matrix; Symm: in Boolean := False); -- A -> Q^tAQ by Householder reflections
  procedure Francis(A,Q: in out Matrix; Symm: in Boolean := False) renames Francis_Lower;
  procedure DeterminantTri(A: in Matrix; D: in out Scalar);         -- determinant of a tridiagonal matrix
  function DeterminantTri(A: Matrix) return Scalar;                 -- determinant of a tridiagonal matrix

  --- some spectral stuff
  procedure LDLdecomp(A: in Matrix; L: in out Matrix; D: in out Diag); -- LDL decomp of Hermitian A
  procedure LDLTransposed(D: in Diag; L: in Matrix ; A: in out Matrix);
  function IsPositive(A: Matrix) return Logical;                       -- for symmetric A
  procedure MaxEigen(A,Tmp: in out Matrix; S: in out Scalar; Pow4: in Positive := 4);
  function MaxEigen(A: Matrix; Pow4: Positive := 4) return Scalar;         -- upper bound on largest |eigenvalue|
  function MinEigen(A: Matrix; Pow4: Positive := 4) return Scalar;         -- MaxEigen for 1/A
  procedure CheckSpecGap(A,Tmp: in out Matrix; Need: in Radius; Gap: out Radius);

  ----------------------------
  ----- purely numerical -----
  ----------------------------

  --- testing and info
  procedure RandomVec(V: in out Vector; Fac: in Flt := One);       -- random vector
  procedure RandomMat(A: in out Matrix; Fac: in Flt := One);       -- random matrix
  procedure RandomOrth(A: in out Matrix);                          -- random orthogonal matrix
  procedure RandomSymm(W: in out Vector; A: in out Matrix);        -- random symmetric A with eigenvalues W
  procedure RandomCircle(H: in out Vector; A: in out Matrix);      -- random cicle A with eigenvalues Wr+-i*Wi
  function SymmErr(A: in Matrix) return Radius;                    -- small if A is symetric
  procedure Show_Approx(N: in String; A: in Matrix);
  procedure CheckSymm(A: in Matrix; ShowMat: Boolean := False);    -- show norm1 of A-A^*
  procedure CheckInv(A,Ai: in Matrix; ShowMat: Boolean := False);  -- show norm1 of A*Ai-Id
  procedure CheckOrth(A: in Matrix; ShowMat: Boolean := False);    -- show norm1 of A*A^*-Id

  --- misc
  procedure Col_Normalize1(Vr,Vi: in out Matrix);        -- approx normalization
  procedure Normalize1(V: in out Vector; N: out Radius); -- approx normalization
  procedure NumSpecGap(A: in out Matrix; Eps: in Radius; R,Err: out Radius);
  function NumSpecGap(A: Matrix; Eps: Radius) return Radius;

  --- sorting
  procedure Default_Sort(W: in out Vector) renames SV_Sorting.Default_Sort;
  procedure Default_Sort(Eps: in Radius; Wr,Wi: in out Vector) renames SV_Sorting.Default_Sort;

private

  procedure CheckDim(N: in Positive; V: in Vector);
  procedure CheckDim(N: in Positive; A: in Matrix);

  procedure Dgefa(A: in out Matrix; Pvt: out Pivot);
  procedure Dgesl(Pvt: in Pivot; A: in Matrix; B: in out Vector);
  procedure Dgedi(Pvt: in Pivot; A: in out Matrix; B: in out Vector);
  procedure HHLeftMult(U,V: in Vector; A: in out Matrix; Tmp: in out Vector);  -- A := Householder(U)*A, use V=U^*
  procedure HHRightMult(U,V: in Vector; A: in out Matrix; Tmp: in out Vector); -- A := A*Householder(U), use V=U^*
  procedure InvNewtonStep(A: in Matrix; Ai,B1,B2: in out Matrix; Err: out Radius);

  Zero_Scalar:   constant Scalar      := Scal(Zero);
  One_Scalar:    constant Scalar      := Scal(One);
  Half_Scalar:   constant Scalar      := Scal(Half);
  SInfo:         constant Scalar_Info := Info(One_Scalar);
  STrunc:        constant Boolean     := SInfo.IsNumeric;
  Not_STrunc:    constant Boolean     := not STrunc;

  pragma Inline (IsZero,SetZero,Scal,Neg,Add,Sub,Mult,AddProd,Adjoint);
  pragma Inline_Always (Swap,Neg,Add,Sub,Mult,AddProd,Scal,CheckDim);

  use SV,SM,SV_Ops;

end ScalVectors;
