{ **********************************************************************
  *                          Unit MATRICES.PAS                         *
  *                            Version 2.4d                            *
  *                   (c) J. Debord, September 2003                    *
  **********************************************************************
  This unit implements vector and matrix operations with dynamic arrays.
  There are 10 types available:

    TVector,     TMatrix      for floating point arrays
    TIntVector,  TIntMatrix   for integer arrays
    TCompVector, TCompMatrix  for complex arrays
    TBoolVector, TBoolMatrix  for boolean arrays
    TStrVector,  TStrMatrix   for string arrays

  To use these arrays in your programs, you must:

  (1) Declare variables of the appropriate type, e.g.

         var
           V : TVector;
           A : TMatrix;

  (2) Allocate each array BEFORE using it:

         DimVector(V, N);         creates vector V[0..N]
         DimMatrix(A, N, M);      creates matrix A[0..N, 0..M]
                                  where N, M are two integer variables

         If the allocation succeeds, all array elements are initialized
         to zero (for numeric arrays), False (for boolean arrays), or
         the null string (for string arrays). Otherwise, the array is
         initialized to NIL.

         The Dim... procedure may be used again to redimension the array.

  (3) Use arrays as in standard Pascal, noting that:

      (a) You cannot use the assignment operator (:=) to copy the
          contents of an array into another array. Writing B := A
          simply makes B point to the same memory block than A. You
          must use one of the provided Copy... procedures (see their
          documentation in the interface part of the unit).

      (b) All arrays begin at index 0, so that the 0-indexed element
          is always present, even if you don't use it.

      (c) A matrix is declared as an array of vectors, so that A[I]
          denotes the I-th vector of matrix A and may be used as any
          vector.

  (4) To deallocate an array, assign the value NIL

         V := nil;

  For more information, read the comments of each routine in the
  interface part of the unit, and check the demo programs.
  **********************************************************************
  References :
  1) 'Mathematiques et Statistiques' by H. Haut (PSI ed.) : GaussJordan
     and related functions
  2) EISPACK (http://www.netlib.org/eispack) : SV_Decomp
  3) 'Matrix Computations' by Golub & Van Loan : QR_Decomp & QR_Solve
     (Pascal implementation contributed by Mark Vaughan)
  ********************************************************************** }

unit matrices;

interface

uses
  fmath;

{ **********************************************************************
  This section defines some error codes.
  ********************************************************************** }

const
  MAT_OK       =   0;  { No error }
  MAT_SINGUL   = - 1;  { Singular matrix }
  MAT_NON_CONV = - 2;  { Non convergence of iterative procedure }
  MAT_NOT_PD   = - 3;  { Matrix not positive definite }

{ **********************************************************************
  This section defines the vector and matrix types.
  ********************************************************************** }

type
  TVector     = array of Float;
  TIntVector  = array of Integer;
  TCompVector = array of Complex;
  TBoolVector = array of Boolean;
  TStrVector  = array of String;

  TMatrix     = array of TVector;
  TIntMatrix  = array of TIntVector;
  TCompMatrix = array of TCompVector;
  TBoolMatrix = array of TBoolVector;
  TStrMatrix  = array of TStrVector;

{ **********************************************************************
  This section defines a function of several variables
  ********************************************************************** }

type
  TFuncNVar = function(X : TVector) : Float;

{ **********************************************************************
  Memory allocation routines
  ********************************************************************** }

procedure DimVector(var V : TVector; Ubound : Integer); overload;
{ ----------------------------------------------------------------------
  Creates floating point vector V[0..Ubound]
  ---------------------------------------------------------------------- }

procedure DimVector(var V : TIntVector; Ubound : Integer); overload;
{ ----------------------------------------------------------------------
  Creates integer vector V[0..Ubound]
  ---------------------------------------------------------------------- }

procedure DimVector(var V : TCompVector; Ubound : Integer); overload;
{ ----------------------------------------------------------------------
  Creates complex vector V[0..Ubound]
  ---------------------------------------------------------------------- }

procedure DimVector(var V : TBoolVector; Ubound : Integer); overload;
{ ----------------------------------------------------------------------
  Creates boolean vector V[0..Ubound]
  ---------------------------------------------------------------------- }

procedure DimVector(var V : TStrVector; Ubound : Integer); overload;
{ ----------------------------------------------------------------------
  Creates string vector V[0..Ubound]
  ---------------------------------------------------------------------- }

procedure DimMatrix(var A : TMatrix; Ubound1, Ubound2 : Integer); overload;
{ ----------------------------------------------------------------------
  Creates floating point matrix A[0..Ubound1, 0..Ubound2]
  ---------------------------------------------------------------------- }

procedure DimMatrix(var A : TIntMatrix; Ubound1, Ubound2 : Integer); overload;
{ ----------------------------------------------------------------------
  Creates integer matrix A[0..Ubound1, 0..Ubound2]
  ---------------------------------------------------------------------- }

procedure DimMatrix(var A : TCompMatrix; Ubound1, Ubound2 : Integer); overload;
{ ----------------------------------------------------------------------
  Creates complex matrix A[0..Ubound1, 0..Ubound2]
  ---------------------------------------------------------------------- }

procedure DimMatrix(var A : TBoolMatrix; Ubound1, Ubound2 : Integer); overload;
{ ----------------------------------------------------------------------
  Creates boolean matrix A[0..Ubound1, 0..Ubound2]
  ---------------------------------------------------------------------- }

procedure DimMatrix(var A : TStrMatrix; Ubound1, Ubound2 : Integer); overload;
{ ----------------------------------------------------------------------
  Creates string matrix A[0..Ubound1, 0..Ubound2]
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Routines for copying vectors and matrices
  ----------------------------------------------------------------------
  Lbound, Ubound   : indices of first and last vector elements
  Lbound1, Lbound2 : indices of first matrix element in each dimension
  Ubound1, Ubound2 : indices of last matrix element in each dimension
  ********************************************************************** }

procedure SwapRows(I, K : Integer; A : TMatrix; Lbound, Ubound : Integer);
{ ----------------------------------------------------------------------
  Exchanges rows I and K of matrix A
  ---------------------------------------------------------------------- }

procedure SwapCols(J, K : Integer; A : TMatrix; Lbound, Ubound : Integer);
{ ----------------------------------------------------------------------
  Exchanges columns J and K of matrix A
  ---------------------------------------------------------------------- }

procedure CopyVector(Dest, Source : TVector; Lbound, Ubound : Integer);
{ ----------------------------------------------------------------------
  Copies vector Source into vector Dest
  ---------------------------------------------------------------------- }

procedure CopyMatrix(Dest, Source : TMatrix;
                     Lbound1, Lbound2, Ubound1, Ubound2 : Integer);
{ ----------------------------------------------------------------------
  Copies matrix Source into matrix Dest
  ---------------------------------------------------------------------- }

procedure CopyRowFromVector(Dest : TMatrix; Source : TVector;
                            Lbound, Ubound, Row : Integer);
{ ----------------------------------------------------------------------
  Copies vector Source into line Row of matrix Dest
  ---------------------------------------------------------------------- }

procedure CopyColFromVector(Dest : TMatrix; Source : TVector;
                            Lbound, Ubound, Col : Integer);
{ ----------------------------------------------------------------------
  Copies vector Source into column Col of matrix Dest
  ---------------------------------------------------------------------- }

procedure CopyVectorFromRow(Dest : TVector; Source : TMatrix;
                            Lbound, Ubound, Row : Integer);
{ ----------------------------------------------------------------------
  Copies line Row of matrix Source into vector Dest
  ---------------------------------------------------------------------- }

procedure CopyVectorFromCol(Dest : TVector; Source : TMatrix;
                            Lbound, Ubound, Col : Integer);
{ ----------------------------------------------------------------------
  Copies column Col of matrix Source into vector Dest
  ---------------------------------------------------------------------- }

{ **********************************************************************
  Vector and matrix functions
  ********************************************************************** }

function Min(X : TVector; Lbound, Ubound : Integer) : Float; overload;
{ ----------------------------------------------------------------------
  Returns the lowest value of vector X
  ---------------------------------------------------------------------- }

function Min(X : TIntVector; Lbound, Ubound : Integer) : Integer; overload;
{ ----------------------------------------------------------------------
  Returns the lowest value of integer vector X
  ---------------------------------------------------------------------- }

function Max(X : TVector; Lbound, Ubound : Integer) : Float; overload;
{ ----------------------------------------------------------------------
  Returns the highest value of vector X
  ---------------------------------------------------------------------- }

function Max(X : TIntVector; Lbound, Ubound : Integer) : Integer; overload;
{ ----------------------------------------------------------------------
  Returns the highest value of integer vector X
  ---------------------------------------------------------------------- }

procedure Transpose(A : TMatrix; Lbound1, Lbound2, Ubound1, Ubound2 : Integer;
                    A_t : TMatrix); overload;
{ ----------------------------------------------------------------------
  Transposes a real matrix
  ----------------------------------------------------------------------
  Input parameters : A       = original matrix
                     Lbound1,
                     Lbound2 = indices of 1st matrix elem. in each dim.
                     Ubound1,
                     Ubound2 = indices of last matrix elem. in each dim.
  ----------------------------------------------------------------------
  Output parameter : A_t     = transposed matrix
  ---------------------------------------------------------------------- }

procedure Transpose(A : TIntMatrix;
                    Lbound1, Lbound2, Ubound1, Ubound2 : Integer;
                    A_t : TIntMatrix); overload;
{ ----------------------------------------------------------------------
  Transposes an integer matrix
  ---------------------------------------------------------------------- }

function GaussJordan(A              : TMatrix;
                     B              : TVector;
                     Lbound, Ubound : Integer;
                     A_inv          : TMatrix;
                     X              : TVector;
                     var Det        : Float) : Integer; overload;
{ ----------------------------------------------------------------------
  Solves the system A * X = B (where B and X are vectors)
  by the Gauss-Jordan method
  ----------------------------------------------------------------------
  Input parameters  : A      = system matrix
                      B      = constant vector
                      Lbound = index of first matrix element
                      Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameters : A_inv  = inverse matrix
                      X      = solution vector
                      Det    = determinant of A
  ----------------------------------------------------------------------
  Possible results  : MAT_OK
                      MAT_SINGUL
  ---------------------------------------------------------------------- }

function GaussJordan(A, B                     : TMatrix;
                     Lbound, Ubound1, Ubound2 : Integer;
                     A_inv, X                 : TMatrix;
                     var Det                  : Float) : Integer; overload;
{ ----------------------------------------------------------------------
  Solves the system A * X = B (where B and X are matrices)
  by the Gauss-Jordan method
  ----------------------------------------------------------------------
  Input parameters  : A       = system matrix
                      B       = constant matrix
                      Lbound  = index of first matrix element in A and B
                      Ubound1 = index of last matrix element in A
                      Ubound2 = index of last matrix element in B
  ----------------------------------------------------------------------
  Output parameters : A_inv  = inverse matrix
                      X      = solution matrix
                      Det    = determinant of A
  ----------------------------------------------------------------------
  Possible results  : MAT_OK
                      MAT_SINGUL
  ---------------------------------------------------------------------- }

function InvMat(A              : TMatrix;
                Lbound, Ubound : Integer;
                A_inv          : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Computes the inverse of a square matrix by the Gauss-Jordan method
  ----------------------------------------------------------------------
  Parameters : as in Gauss-Jordan
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_SINGUL
  ---------------------------------------------------------------------- }

function InvDet(A              : TMatrix;
                Lbound, Ubound : Integer;
                A_inv          : TMatrix;
                var Det        : Float) : Integer;
{ ----------------------------------------------------------------------
  Computes the inverse and the determinant of a square matrix
  by the Gauss-Jordan method
  ----------------------------------------------------------------------
  Parameters : as in Gauss-Jordan
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_SINGUL
  ---------------------------------------------------------------------- }

function Det(A              : TMatrix;
             Lbound, Ubound : Integer) : Float;
{ ----------------------------------------------------------------------
  Computes the determinant of a square matrix by the Gauss-Jordan method
  ----------------------------------------------------------------------
  Parameters : as in Gauss-Jordan
  ----------------------------------------------------------------------
  Notes : (1) This procedure destroys the original matrix A
          (2) If the matrix is quasi-singular, the value 0 is returned
  ---------------------------------------------------------------------- }

function Cholesky(A : TMatrix; Lbound, Ubound : Integer;
                  L : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Cholesky decomposition. Factors the symmetric positive definite matrix
  A as a product L * L', where L is a lower triangular matrix. This
  procedure may be used as a test of positive definiteness.
  ----------------------------------------------------------------------
  Input parameters : A      = matrix
                     Lbound = index of first matrix element
                     Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameter : L      = Cholesky factor of matrix A
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_NOT_PD
  ---------------------------------------------------------------------- }

 function LU_Decomp(A : TMatrix; Lbound, Ubound : Integer) : Integer; overload;
{ ----------------------------------------------------------------------
  LU decomposition. Factors the square matrix A as a product L * U,
  where L is a lower triangular matrix (with unit diagonal terms) and U
  is an upper triangular matrix. This routine is used in conjunction
  with LU_Solve to solve a system of equations.
  ----------------------------------------------------------------------
  Input parameters : A      = matrix
                     Lbound = index of first matrix element
                     Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameter : A      = contains the elements of L and U
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_SINGUL
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

procedure LU_Solve(A : TMatrix; B : TVector; Lbound, Ubound : Integer;
                   X : TVector); overload;
{ ----------------------------------------------------------------------
  Solves a system of equations whose matrix has been transformed by
  LU_Decomp
  ----------------------------------------------------------------------
  Input parameters : A      = result from LU_Decomp
                     B      = constant vector
                     Lbound,
                     Ubound = as in LU_Decomp
  ----------------------------------------------------------------------
  Output parameter : X      = solution vector
  ---------------------------------------------------------------------- }

function LU_Decomp(A : TCompMatrix; Lbound, Ubound : Integer) : Integer; overload;
{ ----------------------------------------------------------------------
  LU decomposition for complex matrices
  ---------------------------------------------------------------------- }

procedure LU_Solve(A : TCompMatrix; B : TCompVector;