
generic

  with function "-"(S: Scalar) return Scalar is <>;
  with function "+"(S1,S2: Scalar) return Scalar is <>;
  with function "-"(S1,S2: Scalar) return Scalar is <>;
  with function "*"(S1,S2: Scalar) return Scalar is <>;
  with function Sqrt(S: Scalar) return Scalar is <>;

  with function "<"(S1,S2: Scalar) return Boolean is <>;
  with function ">"(S1,S2: Scalar) return Boolean is <>;
  with function Min(S1,S2: Scalar) return Scalar is <>;
  with function Max(S1,S2: Scalar) return Scalar is <>;
  with function Sup(S: Scalar) return Flt is <>;
  with function Sqr(S: Scalar) return Scalar is <>;
  with function Scal(K: Integer) return Scalar is <>;
  with function "abs"(S: Scalar) return Scalar is <>;
  with function "/"(S1,S2: Scalar) return Scalar is <>;
  with function Exp(S: Scalar) return Scalar is <>;
  with function Log(S: Scalar) return Scalar is <>;
  with function Get_Precision(R: Scalar) return Positive is <>;

package ScalVectors.RNum is  --- purely numerical, for REAL matrices

  function CountLarge(B: Vector; Cut: Scalar) return Integer;                 -- Number of elements larger than Cut
  procedure Eigen(A: in Matrix; W: in out Vector);                            -- eigenvalues of symmetric A
  procedure Eigen(A: in Matrix; W: in out Vector; V: in out Matrix);          -- eigenvalues/eigenvectors of symmetric A
  procedure SetEigen(V: in Matrix; W: in Vector; A: in out Matrix);           -- A := V*diag(W)*V^{-1}
  procedure Eigen(Ar,Ai: in Matrix; Wr,Wi: in out Vector; Vr,Vi: in out Matrix); -- eigenvals/eigenvecs of general Ar+i*Ai
  procedure Eigen(A: in Matrix; Wr,Wi: in out Vector; Vr,Vi: in out Matrix);     -- eigenvalues/eigenvectors of general A
  procedure QuasiEigen(A: in Matrix; Wr,Wi: in out Vector; Vr: in out Matrix);   -- replace cc eigen by their real part (Wi is just for info)
  procedure ClosestEigen(A: in Matrix; Er,Ei: in out Scalar; Ur,Ui: in out Vector);
  function SpecRadius(A: Matrix; Eps: Radius) return Scalar;                     -- spectral radius (precision Eps)
  procedure SpecProjPos(A: in out Matrix);                                       -- spectral projection for Re(z)>0
  procedure SpecProjInterval(L,R: in Scalar; P: in out Matrix);    -- spectral projection for L<Re(z)<R

private

  function Epsilon(R: Flt) return Flt;
  procedure SortEigen(W: in out Vector; V: in out Matrix);
  procedure SortEigen(Eps: in Flt; Wr,Wi: in out Vector; Vr,Vi: in out Matrix);
  procedure FixEval(Wr,Wi: in out Vector; E: in out Scalar);
  procedure FixEvec(Wi: in Vector; Vr,Vi: in out Matrix);
  procedure NormalizeEvec(Vr,Vi: in out Matrix);

  function NormSmallerOne(A: in Matrix) return Boolean;
  procedure SpecHalf(A,B,C,D: in out Matrix);
  procedure SpecProjPos(A,B,C,D: in out Matrix);

end ScalVectors.RNum;
