with Standard_Integer_Vectors;
with Abstract_Ring,Abstract_Ring.Domain;
with Generic_Vectors,Generic_Matrices;
with Greatest_Common_Divisors;

generic

  with package Ring is new Abstract_Ring(<>);
  with package Domain is new Ring.Domain(<>);
  with package Vectors is new Generic_Vectors(Ring);
  with package Matrices is new Generic_Matrices(Ring,Vectors);
  with package Common_Divisors is new Greatest_Common_Divisors(Ring,Domain);

package Generic_Integer_Linear_Solvers is

-- DESCRIPTION :
--   This package offers a routines for solving linear systems over a domain.

  use Ring,Domain,Matrices;

-- SCALERS :

  function  Scale ( a : Matrix ) return Matrix;
  procedure Scale ( a : in out Matrix; v : out Vectors.Vector );
  procedure Scale ( a : in out Matrix );

  -- DESCRIPTION :
  --   Every row of the matrix will be divided by the greatest commom divisor
  --   of the elements on that row.  The divisors are returned in the vector v
  --   if the appropiate call has been made.

  procedure Scale ( a : in out Matrix; row,col : in integer ); 

  -- DESCRIPTION :
  --   Scales a(row,col..a'last(2)).

-- STATIC TRIANGULATORS :

  procedure Upper_Triangulate ( l : out Matrix; a : in out Matrix );
  procedure Upper_Triangulate ( a : in out Matrix );

  -- DESCRIPTION :
  --   a := l*a, l makes a upper triangular,
  --   l*a becomes then the Hermite normal form of a.

  -- ON ENTRY :
  --   a        an integer matrix.

  -- ON RETURN :
  --   l        the transformation matrix;
  --   a        the triangulated matrix.

  procedure Lower_Triangulate ( a : in out Matrix; u : out Matrix );
  procedure Lower_Triangulate ( a : in out Matrix );

  -- DESCRIPTION :
  --   a := a*u, u makes a lower triangular.

  -- ON ENTRY :
  --   a        an integer matrix.

  -- ON RETURN :
  --   a        the triangulated matrix;
  --   u        the transformation matrix.

-- DYNAMIC TRIANGULATOR :

  procedure Triangulate ( l : in integer; m : in out matrix;
                          ipvt : in out Standard_Integer_Vectors.Vector;
                          piv : out integer );

  -- DESCRIPTION :
  --   Updates lth row of m such that m remains upper triangular.
  --   Note that the last column of m will never be involved in
  --   pivoting operations.  In this way, one can regard the last
  --   column of m as the right hand side of the system a*x = b,
  --   with m = [a|-b].

  -- REQUIRED :  l in m'range(1) and m'range(2) = ipvt'range.

  -- ON ENTRY :
  --   l         indicates which row in m that has to be updated;
  --   m         first l-1 rows are upper triangular;
  --   ipvt      contains the pivoting information;

  -- ON RETURN :
  --   m         m(1..l,m'range(2)) is upper triangular;
  --   ipvt      updated pivoting information;
  --   piv       first entry on lth row of m that was nonzero,
  --             if piv /= l then two columns have been interchanged.
  --             if piv = 0, then m(l,i) = 0, for i < m'last(2).

  procedure Switch ( ipvt : in Standard_Integer_Vectors.Vector;
                     first,last : in integer; m : in out matrix);
  procedure Switch ( ipvt : in Standard_Integer_Vectors.Vector;
                     index : in integer; m : in out matrix );

  procedure Switch ( l,pivot,index : in integer; m : in out matrix );
  procedure Switch ( l,pivot : in integer;
                     first,last : in integer; m : in out matrix );

  -- DESCRIPTION :
  --   Performs column interchangements on m(first..last) or m(index),
  --   with the pivoting information generated by Triangulate,
  --   w.r.t. the lth unknown and pivot, or w.r.t. the pivoting vector ipvt.

-- SELECTORS :

  function Det ( a : Matrix ) return number;

  -- DESCRIPTION :
  --   First the matrix a will be triangulated and then
  --   the determinant of the matrix a will be returned.

  function Per ( a : Matrix; k : Standard_Integer_Vectors.Vector ) 
               return number;
  function Per ( a : Matrix; k : Standard_Integer_Vectors.Vector;
                 max : number ) return number;

  -- DESCRIPTION : Returns the permanent of a matrix.

  -- ON ENTRY
  --   a       (n*m)-matrix, with m <= n;
  --   k       vector(1..m), k(i) indicates the multiplicity
  --           of the ith column;
  --   max     the computation stops when the result >= max.

  function Rank ( a : Matrix ) return natural;

  -- DESCRIPTION :
  --   returns the rank of the matrix a.

-- SOLVERS :

  procedure Solve0 ( a : in Matrix; x : in out Vectors.Vector );

  -- DESCRIPTION :
  --   computes a solution of a*x = 0; where a is upper triangular.
  --   If the number of linear independent rows in a is less than
  --   the number of colums in a, then x will have a nonzero component.

  -- REQUIRED :
  --   a'range(2) = x'range;
  --   Scale(a) should be applied before calling this procedure.

  procedure Solve1 ( a : in Matrix; x : in out Vectors.Vector;
                     b : in Vectors.Vector; fail : out boolean );

  -- DESCRIPTION :
  --   computes the solution of a*x = b, where a is upper triangular.
  --   If Det(a)*Det(a) /= 1, then an integer solutions cannot be
  --   guaranteed, so fail can become true.

  -- REQUIRED :
  --   The matrix a must be square and all ranges must match.

  procedure Solve1 ( a : in Matrix; b : in out Vectors.Vector;
                     fail : out boolean );

  -- DESCRIPTION :
  --   computes the solution of a*x = b, where a is upper triangular,
  --   but after computation, the vector x is stored in b.
  --   If Det(a)*Det(a) /= 1, then an integer solution cannot be
  --   guaranteed, so fail can become true.

  -- REQUIRED :
  --   The matrix a must be square and all ranges must be the same.

  function Solve ( m : matrix; ipvt : Standard_Integer_Vectors.Vector )
                 return Vectors.Vector;

  -- DESCRIPTION :
  --   Solves the system defined by m*x = 0, with m an upper triangular
  --   matrix.  The vector ipvt contains the pivoting information:
  --   ipvt(k) = position of kth column.

  -- REQUIRED : m'range(2) = m'range(1).

  -- ON RETURN : solution vector with nonnegative last component.

end Generic_Integer_Linear_Solvers;