Return to generic_floating_linear_solvers.ads CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Standard_Natural_Vectors;           use Standard_Natural_Vectors;
with Abstract_Ring,Abstract_Ring.Field;
with Generic_Vectors,Generic_Matrices;

generic

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

package Generic_Floating_Linear_Solvers is

-- DESCRIPTION :
--   This package offers a few routines to solve linear systems of equations.
--   The code for lufac, lufco and lusolve is a literal translation from the
--   f77-linpack code.  We distinguish between static triangulators, which
--   take the matrix as a whole at once, and dynamic triangulators, which
--   allow to triangulate row after row.
 
  use Ring,Matrices;

-- STATIC TRIANGULATORS :

  procedure lufac ( a : in out Matrix; n : in natural;
                    ipvt : out Standard_Natural_Vectors.Vector;
                    info : out natural );

  -- DESCRIPTION :
  --   lufac factors a float matrix by gaussian elimination

  --   lufac is usually called by lufco, but it can be called
  --   directly with a saving of time if rcond is not needed.
  --   (time for lufco) = (1 + 9/n)*(time for lufac).

  -- ON ENTRY :
  --   a            a floating point matrix(1..n,1..n) to be factored
  --   n            the dimension of the matrix a

  -- ON RETURN :
  --   a            an upper triangular matrix and the multipliers
  --                which were used to obtain it.
  --                The factorization can be written a = l*u where
  --                l is a product of permutation and unit lower
  --                triangular matrices and u is upper triangular.
  --   ipvt         a natural vector of pivot indices
  --   info         = 0  normal value
  --                = k  if u(k,k) = 0.0.
  --                     This is not an error for this routine,
  --                     but it does indicate that lusolve will
  --                     divide by zero if called.  Use rcond in
  --                     lufco for a reliable indication of singularity.

  procedure lufco ( a : in out Matrix; n : in natural;
                    ipvt : out Standard_Natural_Vectors.Vector;
                    rcond : out number );

  -- DESCRIPTION :
  --   lufco factors a floating point matrix by gaussian elimination
  --   and estimates the condition of the matrix.

  -- If rcond is not needed, lufac is slightly faster.
  -- To solve a*x = b, follow lufco by lusolve.

  -- ON ENTRY :
  --   a           a floating point matrix(1..n,1..n) to be factored;
  --   n           the dimension of the matrix a.

  -- ON RETURN :
  --   a           an upper triangular matrix and the multipliers 
  --               which are used to obtain it.
  --               The factorization can be written a = l*u, where
  --               l is a product of permutation and unit lower triangular
  --               matrices and u is upper triangular.
  --   ipvt        a natural vector of pivot indices
  --   rcond       an estimate of the reciprocal condition of a.
  --               For the system a*x = b, relative perturbations
  --               in a and b of size epsilon may cause relative
  --               perturbations in x of size epsilon/rcond.
  --               If rcond is so small that the logical expression
  --                      1.0 + rcond = 1.0 
  --               is true, than a may be singular to working precision.
  --               In particular, rcond is zero if exact singularity is
  --               detected or the estimate underflows.

  procedure lusolve ( a : in Matrix; n : in natural;
                      ipvt : in Standard_Natural_Vectors.Vector;
                      b : in out Vectors.Vector );

  -- DESCRIPTION :
  --   lusolve solves the system a*x = b using the factors
  --   computed by lufac or lufco

  -- ON ENTRY :
  --   a           a floating-point matrix(1..n,1..n), the output from
  --               lufac or lufco;
  --   n           the dimension of the matrix a;
  --   ipvt        the pivot vector from lufac or lufco;
  --   b           the right hand side vector.

  -- ON RETURN :
  --   b           the solution vector x.

  procedure Triangulate ( a : in out Matrix; n,m : in integer );

  -- DESCRIPTION :
  --   triangulate makes the n*m matrix a triangular using
  --   Gaussian elimination.

  -- ON ENTRY :
  --   a           a floating-point matrix(1..n,1..m);
  --   n           the number of rows of a;
  --   m           the number of columns of a.

  -- ON RETURN :
  --   a           the triangulated matrix.

  procedure Diagonalize ( a : in out Matrix; n,m : in integer );

  -- DESCRIPTION :
  --   Diagonalize makes the n*m floating-point matrix a diagonal.

  -- ON ENTRY :
  --   a           a floating-point matrix(1..n,1..m);
  --   n           the number of rows of a;
  --   m           the number of columns of a.

  -- ON RETURN :
  --   a           the diagonalized matrix.

-- DYNAMIC TRIANGULATORS :

  procedure Upper_Triangulate
                 ( row : in natural; mat : in out Matrix; tol : in number;
                   ipvt : in out Standard_Natural_Vectors.Vector;
                   pivot : out integer );

  -- DESCRIPTION :
  --   Makes the matrix upper triangular by updating the current row of
  --   the matrix.  If pivot = 0 on return, then the matrix is singular.

  -- REQUIRED :
  --   The matrix is upper triangular up to current row, which means that
  --   abs(max(i,i)) > tol, for i in mat'first(1)..row-1.
  --   The matrix might have more columns than rows, but ipvt'range should
  --   match the range of the first columns in the matrix.

  procedure Upper_Triangulate
                 ( roweli : in natural; elim : in Matrix; tol : in number;
                   rowmat : in natural; mat : in out Matrix );
  procedure Upper_Triangulate
                 ( roweli : in natural; elim : in Matrix; tol : in number;
                   firstrow,lastrow : in natural; mat : in out Matrix );

  -- DESCRIPTION :
  --   Using the matrix elim, the unknown roweli is eliminated in mat.

  -- ON ENTRY :
  --   roweli      current unknown to be eliminated;
  --   elim        elim(1..roweli,m'range(2)) is upper triangular;
  --   firstrow    indicates start in range of mat to be updated;
  --   lastrow     indicates end in range of mat to be updated;
  --   roweli      indicates the current unknown to be updated;
  --   mat         the unknows before roweli are already eliminated.

  -- ON RETURN :
  --   mat       the updated matrix after elimination of roweli.

  procedure Switch ( ipvt : in Standard_Natural_Vectors.Vector;
                     row : in integer; mat : in out Matrix );

  -- DESCRIPTION :
  --   Applies the pivoting information ipvt to the row of the matrix.

  procedure Switch ( k,pivot,first,last : in integer; mat : in out Matrix );

  -- DESCRIPTION :
  --   Swaps the column k with the pivot column for row first..last in mat.

  function Solve ( mat : Matrix; tol : number;
                   ipvt : Standard_Natural_Vectors.Vector )
                 return Vectors.Vector;

  -- DESCRIPTION :
  --   Returns a solution to the homogenous linear system with inequalities
  --   in the rows of the matrix.

  -- REQUIRED : the matrix is upper triangular.

  function Solve ( n,col : natural; mat : Matrix )
                 return Vectors.Vector;

  -- DESCRIPTION :
  --   Solves the system with as right-hand side the column indicated
  --   by the parameter col.

  -- REQUIRED :
  --   The matrix must be upper triangular.

end Generic_Floating_Linear_Solvers;