unit xMatrixRN;
{ Feel free to use this unit as you want.
  First developped by Bernard Gisin

  This unit helps manipulating real matrix, of any dimensions,
  just limited by the computer's memory and computation time.
  ----------------------------------------------------------------}

INTERFACE
//"""""""

uses xVectRN, wGlobals2;

//const gglCTEMatrixRNMIndMaxFamilly = gglCTEVecRNIndMaxFamilly;

type
  PRInfoMatrixRNMFloat = ^RInfoMatrixRNMFloat;
  RInfoMatrixRNMFloat  = record
  //===========================
  nSizeRec : int32;
  // Size of memory allocated for this record.
  // By default, the size is  gglnMatRNMRecInfoSize = SizeOf(RInfoMatrixRNMFloat)
  // But it can be changed by program using this unit, and adding
  // information at the end of this information record.
  // Defining a descendent of the matrix with more methodes, will
  // involve the problem that temporary variables will not take these
  // supplementary information into accound, and during assignment will
  // not copy them.
  // Thus it is easier to add information this way than defining a descendent of
  // the matrix

  // In the  LR  decomposition, permuttations.
  panLinePermut       : PTypeArrayInt32;
  nMaxIndexLinePermut : int32; // size of allocation
  nIndexLinePermut    : int32; // max. index defined.


  nMatRank : int32;
  // Rank of the matrix. Determined in the LR decomposition
  end; // RInfoMatrixRNMFloat  = record

type
  PMatrixRNMFloat = ^TMatrixRNMFloat;
  TMatrixRNMFloat = object
  //======================
  bMatFamilly : byte;
  // to which groupe of matrixes this one belongs
  // Memory allocation has been done for index from 1 to nIndMaxLineAbs
     
  fIsTemp : boolean;
  // Normally is FALSE.
  // Is set to TRUE for temporary matrices.

  nIndLine : int32;
  // The matrix has nIndLine lines,

//  nIndCol : int32;
  // Each line vector has nIndCol elements.
  // We have : pavecMat^[nLine]^.nIndMax = nIndCol.

  pavecMat : PArrayRVectRNFloat;
  // The matrix is an array of pointer on vector.
  // Each of these vectors should be considered as a line vector.

  precInfo : PRInfoMatrixRNMFloat;
  // Pointer on a record, storing information on the matrix.

// Private
  // Only used by matrices used for temporary storage.
  // Normally, for internal use. define temporary matrices,
  // for intermediate results.

  fFreeToUse : boolean;
  // defines if this temporary matrix can be used, or is holding an intermediate
  // result.

  pmatNext : PMatrixRNMFloat;
  // Pointer on the next temporary matrix of the same familly

  pOrg : PMatrixRNMFloat;
  // Pointer on the original matrix gglapmatRNMTemp  used by this temporary matrix

  constructor Init(bFamillyP:byte);
  constructor InitWithoutAlloc(bFamillyP:byte);
  // This second constructor is for temporary matrices, without memory alocation

  destructor Done; virtual;

  procedure SetColumns(nNbColP:int32);
  // Define the number of columns of the matrix, i.e.
  // the number of elements of each line vector of the matrix

  procedure SetToIdentity(nLineP,nColP:int32); virtual;
  // set to the identity matrix, having  "1"  on the diagonal,  "0"  else where
  // Has nLineP lines and nColP columns

  procedure CopyFrom(matSrc:TMatrixRNMFloat); virtual;
  // Copy the matrix matSrc into  Self

  procedure Transpose; virtual;
  // Transpose the matrix

  // To define one day ... maybe.
//  function TransposeTo(bNewFamillyP:byte):TMatrixRNMFloat; virtual;
  // Transpose the matrix and give a new matrix as a result.
  // bNewFamillyP  is the familly of the resulting matrix.
  // This familly must have been initialised.
  // It can be the same as  bMatFamilly

  procedure LRDecompose(wFlagsP:word32); virtual;
  // wFlagsP = 0 normally
  // LR decomposition of the matrix.
  // It is the result of a Gaussian elimination.
  // After calling this procedure, we have
  // P*A = L*R  where
  // A is the matrix before the call
  // L is a matrix having one on the diagonal, and zeros above
  // It is stored in pavecMat^[ii].pavV^[jj], for ii > jj
  // R is a matrix having zeroes below the diagonal
  // It is stored in pavecMat^[ii].pavV^[jj], for ii <= jj
  // P is a permuttations, stored in
  // precInfo^.panLinePermut^[1..precInfp^.nIndexLinePermut]
  // The permuttations are :
  // permute line i with line P(i) for i:=1 to precInfp^.nIndexLinePermut
  // If the matrix in not of maximal rank, the decomposition is still valide,
  // but the  R  matrix can have lines with only zeros.
  // Stores the rank of the matrix in  precInfo^.nMatRank
  // If (wFlagsP >= 1) then column exchange are made if necessary
  //    Exchange will be made also on the lines from
  //    nIndLine+1 do nIndLine + wFlagsP-1   i.e.
  //    on  wFlagsP-1  more lines than the original matrix number of lines

  procedure SolvLR(vecSrc:RVectRNFloat; var vecRes:RVectRNFloat; var nErr:int32); virtual;
  // Suppose that an  LR  decomposition has been done before,
  // calling  LRDecompose.
  // If obviously not, because  precInfo^.panLinePermut = NIL, or
  // nIndexLinePermut <>  min (line and column dimension) then
  // call  LRDecompose before.
  // Solves the equation  A*x = b,  where
  // A = the original matrix,
  // b = vecSrc  given as parameter
  // x = vecRes  the result
  // nErr = 0  if there was one solution,
  // nErr = 1 in the rank is less than the min(number of lines, number of cols)
  // nErr = 2 if conditions on  "b"  has to be satisfied because they are more
  //          equations than unkown,
  //          and these conditions don't seem to be satistfied
  //          But a solution is returned. It should be verified by computing
  //          A*x  and comparing with  b.
  // If the numbers of lines is different from the number of columns,
  // there is usually either no solutions, or infinitly many solutions.
  // If the number of columns = Rank then the solution is unique

  end; // TMatrixRNMFloat = object

const
   gglnMatRNMRecInfoSize : int32 = SizeOf(RInfoMatrixRNMFloat);

//  gglvVNEpsilon = 1.0e-12;  See XVectRN.pp

var
  gglanMatRNMLineMax : array[1..gglCTEVecRNIndMaxFamilly] of int32;
  // Memorise the maximum Number of line of each matrx
  // Used to allocate memory for the pointers on the lines

  gglapmatRNMTemp : array[1..gglCTEVecRNIndMaxFamilly] of PMatrixRNMFloat;
  // Memory allocation is done automatically, if temporary matrices are needed.
  // For example in adding matrices an multiplying matrices,
  // a temporary matrix is needed.

procedure VNCreateNewMatRNMFamilly(bFamillyP:byte; nIndLineMaxP:int32);
//=====================================================================
// This creates a new familly of matrices, which has all the same
// amount of memory allocated for lines.
// This must be called before using the matrices.
// xVectRn.VNCreateNewVecRNFamilly(bFamillyP, nIndMaxP)  must have been called
// before, with the same number designing the familly.
// nIndMaxP  will define the number of columns of the matrix.
// bFamilly, in [1..10] define which groupe of matrices we want to define.
// bFamilly defines also the size of the columns since it will
// defines the familly of each line vector of the matrix.
// for example,  bFamilly=1  and nIndLineMaxP=10
// will define a familly of matrices having at maximum 10 lines.
// The number of columns depend on the familly of vectors.
// pavecMat^[0] is allocated for a line vector, but is not part of the matrix.
// It will be used for intermediate calculs.

procedure VNDisposeMatRNMFamilly(bFamillyP:byte);
//==============================================
// If  CreateNewMatRNMFamilly(bFamillyP)  has been called then memory might have
// been allocated, and must be disposed by calling this procedure.
//  Memory could have been allocated, to store intermediate matrix result,
// as for example, in addition of matrices and multiplication of matrices.
// Thus call this procedure when this matrix familly is no more needed.

procedure SetTemporaryMat(var matM1,matM2,matM3:TMatrixRNMFloat);
//===============================================================
// matM3 = temporary matrix.
// free the vectors matM1 and matM2 if they are temporary
// Used in matrix additions, substraction and mult..
// Normally not used by an external unit.

//operator := (matIn:TMatrixRNMFloat) matOut:TMatrixRNMFloat;
// The assignment operator doesn't work as I would like.
// I would like the assignment operator := to be defined as any comparaison
// operator, but without returning any results. For example :
// operator := (var rOut:TMatrixRNMFloat; rIn:TMatrixRNMFloat);
// rOut.blabla:=rIn.blabla
// Thus I will use the  <=  operator as an assignment operator,
// ignoring the boolean result.

operator <= (var matRes:TMatrixRNMFloat; matSrc:TMatrixRNMFloat) fRes:boolean;
//============================================================================
// It is used as an assignment operator, matRes:=matSrc,
// ignoring the result.
// The :  operator := (matSrc:TMatrixRNMFloat) matRes:TMatrixRNMFloat;
// cannot be redefined, and should not be used to assign one matrix to another.
// DO NOT USE mat1:=mat2
//---------------*****----
// mat1:=mat2;  would assign the pointer  mat2.pavecMat  to  mat1.pavecMat
// the  mat1.pavecMat  pointer would be lost.
// If matSrc is the result of an operation which has been stored in
// gglapmatRNMTemp,  then this exanges the two pointers
// matSrc.pavecMat and matRes.pavecMat
// Else the content of matSrc is copied in matRes, but not its bFamilly.

operator + (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
operator - (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;

operator * (vScal:Float; matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
// number time a Matrix

operator * (matM1:TMatrixRNMFloat; vecV2:RVectRNFloat) vecV3:RVectRNFloat;
// Matrix time a vector

operator * (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
// Matrix Multiplication

IMPLEMENTATION
//""""""""""""

Uses //wUtiles0,
  wObjectCreatorLOW2
 ,wDebug0;

const fFirstInit : boolean = TRUE; // Should the first initialisation be done ?

var glf : boolean; // To use  <=  as an assignment operator
//var nNotUsed : integer; // only to avoid an error message while compliling in DevPascal 1.9

procedure VNCreateNewMatRNMFamilly(bFamillyP:byte; nIndLineMaxP:int32);
//=====================================================================
var
  ii : int32;
begin
if fFirstInit then begin
   // Indicate that no familly has been defined yet,
   // and that no temporary matrices have been allocated yet.
   for ii:=1 to gglCTEVecRNIndMaxFamilly do begin
      gglanMatRNMLineMax[ii] :=0;
      gglapmatRNMTemp[ii]   :=NIL; // No temporary matrices have been defined yet
      end;
   fFirstInit:=FALSE;
   end;

// One familly of matrices is defines.
// The corresponding familly of vectors must also be defined
gglanMatRNMLineMax[bFamillyP]:=nIndLineMaxP;
end;  // VNCreateNewMatRNMFamilly

procedure VNDisposeMatRNMFamilly(bFamillyP:byte);
//===============================================
var
  pmatM,
  pmatMLast : PMatrixRNMFloat; // used to search a free temp. vector
begin
// Get the first matrix of the familly
pmatM:=gglapmatRNMTemp[bFamillyP];

// Loop on all chained matrices of the same familly
while (pmatM <> NIL) do begin
   pmatMLast :=pmatM;
   pmatM     :=pmatM^.pmatNext;
   Dispose(pmatMLast, Done);
   end;

// Just to indicate that this vector has not been initialised
gglapmatRNMTemp[bFamillyP]    :=NIL;
gglanMatRNMLineMax[bFamillyP]  :=0; // this familly has not been definedend;  // VNCreateNewMatRNFamilly
end;  // VNDisposeMatRNMFamilly

//****************************************************************
// TMatrixRNMFloat = object
//****************************************************************

constructor TMatrixRNMFloat.Init(bFamillyP:byte);
//===============================================
var
  ii : int32; // To loop over the lines of the matrix
begin
bMatFamilly :=bFamillyP;
nIndLine    :=0;
//nIndCol     :=0;

fIsTemp     :=FALSE;
fFreeToUse  :=FALSE;
pmatNext    :=NIL;
pOrg        :=@Self;

// Memory allocation for the array of line vectors of the matrix
GetMem(pavecMat, (gglanMatRNMLineMax[bFamillyP]+1)*SizeOf(RVectRNFloat));

// The line 0 is also initialised, but use for temporary storage
// Init each vector corresponding to one line of the matrix
for ii:=0 to gglanMatRNMLineMax[bFamillyP] do VNInit(pavecMat^[ii], bFamillyP);

if (gglnMatRNMRecInfoSize > 0) then begin
   GetMem(precInfo, gglnMatRNMRecInfoSize);
   precInfo^.nSizeRec:=gglnMatRNMRecInfoSize;
   precInfo^.panLinePermut:=NIL;
   precInfo^.nMaxIndexLinePermut:=0;
   precInfo^.nMatRank:=0;
   end
else precInfo:=NIL;
end;  //  TMatrixRNMFloat.Init

constructor TMatrixRNMFloat.InitWithoutAlloc(bFamillyP:byte);
//===========================================================
// Used for temporary matrices.  Destructor Done will not be called for them.
begin
bMatFamilly :=bFamillyP;
nIndLine    :=0;
//nIndCol     :=0;

fIsTemp     :=TRUE;
fFreeToUse  :=FALSE;
//pmatNext    :=NIL;
//pOrg        :=NIL;
//precInfo    :=NIL;

// No memory allocation.
end;  //  TMatrixRNMFloat.InitWithoutAlloc

destructor TMatrixRNMFloat.Done;
//==============================
var
  ii : int32; // To loop over the lines of the matrix
begin
if fIsTemp then begin
   pOrg^.fFreeToUse:=TRUE;
   exit; // temporary matrices has nothing to free
   end;

// Free memory for the record information
if (precInfo <> NIL) then begin
   if (precInfo^.panLinePermut <> NIL) then
      FreeMem( precInfo^.panLinePermut,
               (precInfo^.nMaxIndexLinePermut+1)*SizeOf(int32));

   FreeMem(precInfo, precInfo^.nSizeRec);
   end;

// Free memory for each vector corresponding to one line of the matrix
for ii:=gglanMatRNMLineMax[bMatFamilly] downto 0 do VNDone(pavecMat^[ii]);

// Free the matrix
FreeMem(pavecMat, (gglanMatRNMLineMax[bMatFamilly]+1)*SizeOf(RVectRNFloat));
pavecMat:=NIL;
end;  //  TMatrixRNMFloat.Done

procedure TMatrixRNMFloat.SetColumns(nNbColP:int32);
//==================================================
// Define the number of columns of the matrix, i.e.
// the number of elements of each line vector of the matrix
var
  ii : int32; // To loop over the lines of the matrix
begin
// Loop on each line and define the number of element they have which are defined
for ii:=0 to gglanMatRNMLineMax[bMatFamilly] do pavecMat^[ii].nIndMax:=nNbColP;
end;  // TMatrixRNMFloat.SetColumns

procedure TMatrixRNMFloat.SetToIdentity(nLineP,nColP:int32);
//==========================================================
// set to the identity matrix, having  "1"  on the diagonal,  "0"  else where
// Has nLineP lines and nColP columns
var
  ii : int32; // To loop over the lines of the matrix
begin
// No range checking...

// Loop on each line and define set their length and value
for ii:=1 to nLineP do VNSetToBasis(pavecMat^[ii], nColP, ii);

nIndLine:=nLineP;
end;  // TMatrixRNMFloat.SetToIdentity

procedure TMatrixRNMFloat.CopyFrom(matSrc:TMatrixRNMFloat);
//=========================================================
// Copy the matrix matSrc into  Self
var ii : int32;
begin
// Copy the lines
for ii:=1 to matSrc.nIndLine do VNCopyFromTo(matSrc.pavecMat^[ii], pavecMat^[ii]);

nIndLine:=matSrc.nIndLine;
precInfo^.nMatRank:=matSrc.precInfo^.nMatRank;

// To copy this, memory allocation should be done.
//  panLinePermut
//  nMaxIndexLinePermut
//  nIndexLinePermut
end;  // TMatrixRNMFloat.CopyFrom

procedure TMatrixRNMFloat.Transpose;
//==================================
// Transpose the matrix
// Not all overflow checking are made.
// The number of lines is the minimum length of the line vectors.
// They should normally be all of the same length
var
  ii,jj,      // To loop over the lines and columns of the matrix
  nColLen,    // Minimum line length, i.e. number of columns of the matrix
  nMinLine    : int32;  // take the minimum length of the line vectors
  vTemp       : Float; // to swap lines and columns
begin
// maximum number of lines the transposed matrix can have
nColLen:=gglanMatRNMLineMax[bMatFamilly];

// First loop to see what is the shortest line length
for ii:=1 to nIndLine do
   if (nColLen > pavecMat^[ii].nIndMax) then nColLen:=pavecMat^[ii].nIndMax;

// The transposed matrix can not have line length longer than
// gglanVectRNIndMax[bMatFamilly], thus cut the lines if necessary
if (nIndLine > gglanVectRNIndMax[bMatFamilly]) then
    nIndLine:= gglanVectRNIndMax[bMatFamilly];

// To loop first on the square part of the matrix
if (nColLen <= nIndLine) then nMinLine:=nColLen
else                          nMinLine:=nIndLine;
for ii:=1 to nMinLine do begin
   for jj:=ii+1 to nMinLine do begin
      // swap Aij with Aji
      vTemp                  :=pavecMat^[jj].pavV^[ii];
      pavecMat^[jj].pavV^[ii]:=pavecMat^[ii].pavV^[jj];
      pavecMat^[ii].pavV^[jj]:=vTemp;
      end;
   end; // for ii:=1 to nMinLine do begin

if (nColLen <= nIndLine) then begin
   // Matrix has more lines than colums   nMinLen = nColLen
   for ii:=nColLen+1 to nIndLine do begin
      for jj:=1 to nColLen do pavecMat^[jj].pavV^[ii]:=pavecMat^[ii].pavV^[jj];
      end;
   end
else begin
   // Matrix has less lines than colums   nMinLen = nIndLine
   for ii:=1 to nIndLine do begin
      for jj:=nIndLine+1 to nColLen do
          pavecMat^[jj].pavV^[ii]:=pavecMat^[ii].pavV^[jj];
      end;
   end;

// length of lines of the transposed matrix
for jj:=1 to nColLen do pavecMat^[jj].nIndMax:=nIndLine;

// Number of lines of the transposed matrix = the minimum length of the
// line vectors.
nIndLine:=nColLen;
end;  // TMatrixRNMFloat.Transpose

//function TMatrixRNMFloat.TransposeTo(bNewFamillyP:byte):TMatrixRNMFloat;
//======================================================================
// Transpose the matrix and give a new matrix as a result.
// bNewFamillyP  is the familly of the resulting matrix.
// This familly must have been initialised.
// It can be the same as  bMatFamilly
//begin
// TO define ...
//end;  // TMatrixRNMFloat.TransposeTo

procedure TMatrixRNMFloat.LRDecompose(wFlagsP:word32);
//====================================================
var
  k1,k2,      // (k1,k2) = diagonal of the matrix, to do the elimination
              // Or, if the rank of the square sub matrix is not maxial
              // then the elimination will continue next to the diagonal
  k3,         // Used if column permuttation is wanted
  ii,jj,      // Line and column index
  nDimColL,   // number of columns
  nLinePermL  // Line permutted with the line k1
              : int32;
  vMaxColL,   // Maximum value of a column, during Gaussian elimination
  vPivotL,
  vMat_ik     // Element of the  "L"  matrix
              : Float;
  pavV1L,     // to permute lines and to do gaussian elimination
  pavV2L      : PTypeArrayFloat;

begin
// Suppose that all line vectors has the same length
nDimColL:=pavecMat^[1].nIndMax;

// Memory allocation to store the permuttations if necessary
with precInfo^ do if (nMaxIndexLinePermut < nIndLine) then begin
   // Need to allocate more memory for the permuttations
   if (nMaxIndexLinePermut > 0) then // first free what was allocated
      FreeMem(panLinePermut, (nMaxIndexLinePermut+1)*SizeOf(int32));

   // Then allocate memory
   nMaxIndexLinePermut:=nIndLine;
   GetMem(panLinePermut, (nMaxIndexLinePermut+1)*SizeOf(int32));
   end; // with precInfo^ do if ....

// Loop over the "diagonal" (k1,k2), to do Gaussian elimination
// (k1, k2) = line, column   (k1 <= k2)
k3:=nDimColL;
k2:=0;
k1:=0;
repeat
   inc(k1);  inc(k2);  // k1=k2 if no rank limitation occured

   // Look for a "pivot" and memorise line permuttation
   vMaxColL:=0.0;
   for ii:=k1 to nIndLine do begin
      if (vMaxColL < abs(pavecMat^[ii].pavV^[k2])) then begin
         nLinePermL:=ii;
         vMaxColL:=abs(pavecMat^[ii].pavV^[k2]);
         end;
      end; // for ii:=k1 to nIndLine do begin

   // Test if the column (k1..nIndLin, k2) = 0
   if (vMaxColL <= gglvVNEpsilon) then begin
      // The rank of the square sub matrix is not maximal
      if (wFlagsP = 0) then begin
         // Indicate that the column is zero
         pavecMat^[k1].pavV^[k2]:=0.0;

         // Not necessary, but to see, set the column to zero
         for ii:=k1+1 to nIndLine do pavecMat^[ii].pavV^[k2]:=0.0;
         end
      else begin
         // Columns permuttation to continue the decomposition, until the
         // maximum rank has been found
         // Exange column  k2  with  column  k3
         // Do the exchange also on  wFlagsP-1  more lines.
         if (k2 < k3) then for ii:=k1 to nIndLine + int32(wFlagsP-1) do begin
             vPivotL                 :=pavecMat^[ii].pavV^[k2];
             pavecMat^[ii].pavV^[k2] :=pavecMat^[ii].pavV^[k3];
             pavecMat^[ii].pavV^[k3] :=vPivotL;
             end;

         // Next time, will exchange with another column
         dec(k3);

         // Try with this column to do gaussian elimination
         dec(k2);
         end; // else of if (wFlagsP = 0) then begin

      // Try the elimination with the next column.
      dec(k1);
      end
   else begin
      // memorise the line permuttation
      precInfo^.panLinePermut^[k1] :=nLinePermL;

      // Do the line permuttation. Only necessary if  nLinePermL <> k1
      if (nLinePermL <> k1) then for jj:=k2 to nDimColL do begin
         vPivotL                         :=pavecMat^[k1].pavV^[jj];
         pavecMat^[k1].pavV^[jj]         :=pavecMat^[nLinePermL].pavV^[jj];
         pavecMat^[nLinePermL].pavV^[jj] :=vPivotL;
         end;

      // Line to multiply and substract to the Line  "ii"
      pavV1L:=pavecMat^[k1].pavV;

      // One step of the LR decomposition
      vPivotL:=pavV1L^[k2];  // = A[k1,k2] = pavecMat^[k1].pavV^[k2]

      // Loop on the lines to combine eliminate the element A[ii,k2], k2 fixed here
      for ii:=k1+1 to nIndLine do begin
         pavV2L :=pavecMat^[ii].pavV;  // line to be modified
         vMat_ik:=pavV2L^[k2] / vPivotL; // = A[ii,k2] / A[k1,k2]

         // Store the "L"  part of the decomposition
         pavV2L^[k2]:=vMat_ik; // in  A[ii,k2] = pavecMat^[ii].pavV^[k2]

         // Loop on the columns
         for jj:=k2+1 to nDimColL do begin
            // A[ii,jj] := A[ii,jj]   -  A[k1,jj]   * A[ii,k2] / A[k1,k2] }
            pavV2L^[jj] :=pavV2L^[jj] - pavV1L^[jj] * vMat_ik;
            end; // for jj:=k2+1 to nDimColL do begin
         end;  // for ii:=k1+1 to nIndLine do begin
      end; // else of if (vMaxColL <= gglvVNEpsilon) then begin
   until (k1 >= nIndLine) or (k2 >= k3);  // k3=nDimColL  if (wFlagsP=0)

// Number of indexes having a meaning in the permuttation array
precInfo^.nIndexLinePermut:=k1;

// Store the Rank of the matrix.
precInfo^.nMatRank:=k1
end;  // TMatrixRNMFloat.LRDecompose

procedure TMatrixRNMFloat.SolvLR(vecSrc:RVectRNFloat; var vecRes:RVectRNFloat; var nErr:int32);
//=============================================================================================
var
  k1,k2,      // (k1,k2) = diagonal of the matrix, to do the elimination
              // Or, if the rank of the square sub matrix is not maxial
              // then the elimination will continue next to the diagonal
  ii,jj,      // Line and column index
  nDimColL,   // number of columns
  nMax_k2,    // X[kk] = 0  for kk > nMax_k2
  nLinePermL  // Line permutted with the line jj
              : int32;
  vStoreL,
  vPivotL,
  vSumL       // For the multipl.  R * x
              : Float;
  pavV1L,     // One line of the  R  matrix
  pavVBL,     // the B vector
  pavVXL      // The X vector  (They can point to the same value)
              : PTypeArrayFloat;
  anColSkippedL : array[1..100] of int32;
              // Must memorise the columns which where skipped.
              // 100 => suppose that nIndLine - Rank <= 100
  nCS         : int32; // =number of columns skipped
  vEpsilonL   : Float; // To test the existance of a solution
begin
// It is supposed that nIndLine = vecSrc.nIndMax

// Test to see if the LRDecompose has been done
if (precInfo^.panLinePermut = NIL) or
   (precInfo^.nIndexLinePermut > nIndLine) then begin
   // The LR decomposition must be done
   LRDecompose(0);
   end;

// Suppose that all line vectors has the same length
nDimColL:=pavecMat^[1].nIndMax;

// Dimension of the result vector
vecRes.nIndMax:=nDimColL;

pavVBL:=vecSrc.pavV;
pavVXL:=vecRes.pavV;

// If the matrix Rank is less than the number of lines and columns,
// No numerical solutions can be found
if (precInfo^.nMatRank < nIndLine) and (precInfo^.nMatRank < nDimColL) then begin
   // No solutions
   for ii:=1 to nDimColL do pavVXL^[ii]:=0;
   nErr:=1;
   exit;
   end;

if (precInfo^.nMatRank < nIndLine) then begin
   // Compute the max abs value of  pavVBL  to test if a sol. existe or not
   vEpsilonL:=0.0;
   for ii:=1 to vecSrc.nIndMax do
      if (vEpsilonL < abs(pavVBL^[ii])) then vEpsilonL := abs(pavVBL^[ii]);

   vEpsilonL:=vEpsilonL * gglvVNEpsilon;
   end
else vEpsilonL:=gglvVNEpsilon; // not used

if (vecSrc.pavV <> vecRes.pavV) then begin
   // Copy the content of the source vector into the result vector,
   // in order to not modify the source vector.
   for ii:=1 to vecSrc.nIndMax do pavVXL^[ii]:=pavVBL^[ii];
   end;
// Else they are the same, and nothing to do.

pavVBL:=pavVXL; // they are maintained for the reading of the code

// Loop over the "diagonal" (k1,k2), to do the equivalent of Gaussian elimination
// (k1, k2) = line, column   (k1 <= k2)
k2:=0;
k1:=0;
nCS:=0;
repeat
   inc(k1);  inc(k2);  // k1=k2 if no rank limitation occured

   if (pavecMat^[k1].pavV^[k2] = 0.0) then begin
      // This column was skipped
      dec(k1);

      // Memorise the skipped column
      if (nCS < 100) then inc(nCS); // else the result will be wrong, but...
      anColSkippedL[nCS]:=k2;
      end
   else begin
      // Do the permuttation equivalent of that in the LRDecompose
      nLinePermL :=precInfo^.panLinePermut^[k1];
      vStoreL             :=pavVBL^[nLinePermL];
      pavVBL^[nLinePermL] :=pavVBL^[k1];
      pavVBL^[k1]         :=vStoreL;

      // Do the equivalent of the LRDecomposition
      vPivotL:=pavVBL^[k1];
      for ii:=k2+1 to nIndLine do begin
         // b[ii]   := b[ii]      - A[ii,k2]                * b[k1]
         pavVBL^[ii]:=pavVBL^[ii] - pavecMat^[ii].pavV^[k2] * vPivotL;
         end;
      end; // else of if ...
   until (k1 >= nIndLine) or (k2 >= nDimColL);

// Test if the system has a solution.
// If Dim(Lines) <= Rank (=k1), the loop is not entered
nErr:=0;
for ii:=k1 + 1 to vecSrc.nIndMax do begin  // vecSrc.nIndMax = nIndLine is supposed
//DBW(1, 'b=' + FTS(Ln(abs(pavVBL^[ii]) + 1e-30) / Ln(10)));
   if (abs(pavVBL^[ii]) > vEpsilonL) then begin
      // Numerically, no solutions exist, but an answer will be given.
      // The test that A*x=b can be made by the user.
      nErr:=2;
      end;
   end;

nMax_k2:=k2;

// Sets to zero elements of X if the matrix has more columns than its Rank
for ii:=nMax_k2+1 to nDimColL do pavVXL^[ii]:=0;

// Solve the system
//""""""""""""""""""
// pavecMat^[k1].pavV^[k2] <> 0  unless k1=0
while (k1 > 0) and (k2 > 0) do begin
   pavV1L:=pavecMat^[k1].pavV;

   if (pavV1L^[k2] = 0.0) then begin
     // This column was skipped
      // pavXL^[k2] can have any value, but set it to zero
      pavVXL^[k2]:=0;
      inc(k1); // go to the previous column
      end
   else begin
      // Compute R * X, for the known X
      vSumL:=0.0;
      for jj:=k2+1 to nMax_k2 do begin
         //sum:=sum   + A[k1,jj]    * X[jj] }
         vSumL:=vSumL + pavV1L^[jj] * pavVXL^[jj];
         end;

      // x[k2]   :=(      b[k1] -   sum) /    A[k1,k2]
      pavVXL^[k2]:=(pavVBL^[k1] - vSumL) / pavV1L^[k2];
      end;

   // Skip the columns which where skipped while doing the elimination
   while (nCS > 0) and (anColSkippedL[nCS] = k2-1) do begin
      dec(k2); dec(nCS);

      // pavXL^[k2] can have any value, but set it to zero
      pavVXL^[k2]:=0;
      end; // while

   dec(k1);  dec(k2);
   end; // while (k1 > 0) and (k2 > 0) do begin

end;  // TMatrixRNMFloat.SolvLR

//****************************************************************
// Operators definitions
//****************************************************************

procedure SetTemporaryMat(var matM1,matM2,matM3:TMatrixRNMFloat);
//===============================================================
// matM3 = temporary matrix.
// free the vectors matM1 and matM2 if they are temporary
var
  pmatMNew,  // If a new temporary matrix should be created
  pmatMLast : PMatrixRNMFloat; // used to search a free temp. matrix
begin
// Initialise the new vector matM3
matM3.InitWithoutAlloc(matM1.bMatFamilly);  // = matM2.bMatFamilly  normally
// matM3.fIsTemp      :=TRUE;

// if matM1 or matM2 are temporary matrices then they are no more used
if (matM1.fIsTemp) then begin
   // Use matM1 for the temporary result
   matM3.pavecMat :=matM1.pavecMat;
   matM3.pOrg     :=matM1.pOrg; // Pointer on  gglapmatRNMTemp[..] or one of the list
   matM3.precInfo :=matM1.precInfo;
   end
else if (matM2.fIsTemp) then begin
   // Use matM2 for the temporary result
   matM3.pavecMat :=matM2.pavecMat;
   matM3.pOrg     :=matM2.pOrg; // Pointer on  gglapmatRNMTemp[..]  or one of the list
   matM3.precInfo :=matM2.precInfo;
   end
else begin
   // Need a temporary matrix to store the result
   // look for the first temporary free matrix of the familly
   pmatMLast :=NIL;
   pmatMNew  :=gglapmatRNMTemp[matM1.bMatFamilly];
   while (pmatMNew <> NIL) and (not pmatMNew^.fFreeToUse) do begin
         pmatMLast :=pmatMNew;
         pmatMNew  :=pmatMNew^.pmatNext;
         end;

   if (pmatMNew = NIL) then begin
      // No free temporary matrix, a new one must be created
      // This is seldom necessary
      pmatMNew:=gglpoObjCreate^.CreateMatrixRNMFloat(matM1.bMatFamilly);
//      new(pmatMNew, Init(matM1.bMatFamilly));
//      pmatMNew^.fIsTemp    :=FALSE;  it is used by temp matrices, but it is not a temp matrix
      pmatMNew^.pmatNext   :=NIL;
      pmatMNew^.pOrg       :=pmatMNew;

      // Add the new matrix to the chained list of temporary matrix of this familly
      if (pmatMLast <> NIL) then pmatMLast^.pmatNext:=pmatMNew
      else                       gglapmatRNMTemp[matM1.bMatFamilly]:=pmatMNew;
      end; // if (pmatMNew = NIL) then begin

   // Use for matM3, the array of the temporary matrix pamatMNew
   matM3.pavecMat :=pmatMNew^.pavecMat;
   matM3.pOrg     :=pmatMNew;
   matM3.precInfo :=pmatMNew^.precInfo;

   // Mark that  pmatMNew  is not free to use
// matM3.pOrg^.fFreeToUse:=FALSE;  must be done in the operator +,-,*
   end; // if .. Need a temporary vector to store the result
end;  // SetTemporaryMat

operator <= (var matRes:TMatrixRNMFloat; matSrc:TMatrixRNMFloat) fRes:boolean;
//============================================================================
// This defines the assignment operator, since the normal assignment
// operator cannot be defined.
begin
// For assignment of the kind :
// glf:=mat1 <= mat2 + mat3, the test will be TRUE.
// glf is a side effect product to ignore.
if (matSrc.fIsTemp) then begin
   // we know that the  matSrc  matrix use the data of
   // matSrc.pOrg^.pavecMat
   // Exchange simply  matSrc.pOrg^.pavecMat  with  matRes.pavecMat
   // Exchange also the familly, in case they are different
   matSrc.pOrg^.pavecMat    :=matRes.pavecMat;
   matRes.pavecMat          :=matSrc.pavecMat;
   matSrc.pOrg^.precInfo    :=matRes.precInfo;
   matRes.precInfo          :=matSrc.precInfo;
   matSrc.pOrg^.bMatFamilly :=matRes.bMatFamilly;
   matRes.bMatFamilly       :=matSrc.bMatFamilly;
   matRes.nIndLine          :=matSrc.nIndLine;
//   matRes.nIndCol           :=matSrc.nIndCol;

   // It is free for another use
   MatSrc.Done;

   // matSrc itself is a temporary variable, which will disapear.
   // It has not allocated memory, thus does not call its destructor Done

   fRes:=FALSE; // not important
   exit;
   end;

// Here we copy simply one matrix into another

// Copy the content of the vecSrc vector in the vecRes vector
matRes.CopyFrom(matSrc);
fRes:=TRUE; // not important
end; // operator <=

operator + (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
//==============================================================
var
   ii,niMin,niMax : int32;
begin
// matM3 = temporary matrix.
// free the matrices matM1 and matM2 if they are temporary
SetTemporaryMat(matM1,matM2,matM3);

// Selecting the shortest number of lines
if (matM1.nIndLine < matM2.nIndLine) then begin
   niMin:=matM1.nIndLine;
   niMax:=matM2.nIndLine;

   // The matrix having the least number of lines
   // is taken as if it was completed with lines of zeros.
   for ii:=niMin+1 to niMax do glf:=matM3.pavecMat^[ii]<=matM2.pavecMat^[ii];
   end
else begin
   niMin:=matM2.nIndLine;
   niMax:=matM1.nIndLine;

   // The shortest number of lines is taken as if it was completed with zeros.
   for ii:=niMin+1 to niMax do glf:=matM3.pavecMat^[ii]<=matM1.pavecMat^[ii];
   end;

// Addition of the matrices, i.e. addition of the line vectors
for ii:=1 to niMin do
   glf:=matM3.pavecMat^[ii]<=matM1.pavecMat^[ii] + matM2.pavecMat^[ii];

matM3.nIndLine  :=niMax; // Length of the resulting matrix

// destructor are called, they do little for temp matrices
if (matM1.fIsTemp) then MatM1.Done;
if (matM2.fIsTemp) then MatM2.Done;

// Indicate that this is not free to use
matM3.pOrg^.fFreeToUse:=FALSE;
end; // operator +

operator - (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
//==============================================================
var
   ii,niMin,niMax : int32;
begin
// matM3 = temporary matrix.
// free the matrices matM1 and matM2 if they are temporary
SetTemporaryMat(matM1,matM2,matM3);

// Selecting the shortest vector length
if (matM1.nIndLine < matM2.nIndLine) then begin
   niMin:=matM1.nIndLine;
   niMax:=matM2.nIndLine;

   // The matrix having the least number of lines
   // is taken as if it was completed with lines of zeros.
   for ii:=niMin+1 to niMax do glf:=matM3.pavecMat^[ii]<=-1.0*matM2.pavecMat^[ii];
   end
else begin
   niMin:=matM2.nIndLine;
   niMax:=matM1.nIndLine;

   // The shortest line number is taken as if it was completed with zeros.
   for ii:=niMin+1 to niMax do glf:=matM3.pavecMat^[ii]<=matM1.pavecMat^[ii];
   end;

// Substraction of the matrices, i.e. substraction of the line vectors
for ii:=1 to niMin do
   glf:=matM3.pavecMat^[ii]<=matM1.pavecMat^[ii] - matM2.pavecMat^[ii];

matM3.nIndLine  :=niMax; // Length of the resulting matrix

// destructor are called, they do little for temp matrices
if (matM1.fIsTemp) then MatM1.Done;
if (matM2.fIsTemp) then MatM2.Done;

// Indicate that this is not free to use
matM3.pOrg^.fFreeToUse:=FALSE;
end; // operator -

operator * (vScal:Float; matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
//====================================================================
// number time a Matrix
var
   ii : int32;
begin
// matM3 = temporary matrix.
// free the matrice matM2 if it is temporary
SetTemporaryMat(matM2,matM2,matM3);

// Multiplication of each line by the scalar
for ii:=1 to matM2.nIndLine do begin  // Loop on lines
   matM3.pavecMat^[ii]:=vScal * matM2.pavecMat^[ii];
   end; // for ii  lines

// Number of lines of the resulting matrix
matM3.nIndLine  :=matM2.nIndLine;

// destructor are called, they do little for temp matrices
if (matM2.fIsTemp) then MatM2.Done;

// Indicate that this is not free to use
matM3.pOrg^.fFreeToUse:=FALSE;
end; // operator *  number time matrix

operator * (matM1:TMatrixRNMFloat; vecV2:RVectRNFloat) vecV3:RVectRNFloat;
//========================================================================
// Matrix time a vector
var ii : int32;     // to loop over lines
begin
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
SetTemporaryVect(matM1.pavecMat^[1], matM1.pavecMat^[1], vecV3);
// In fact, this uses a temporary vector  vecV3.
// Cannot use vecV2, even if it is temporary

// do the multiplication
for ii:=1 to matM1.nIndLine do vecV3.pavV^[ii]:=matM1.pavecMat^[ii] * vecV2;
vecV3.nIndMax:=matM1.nIndLine;

// Free temporary matrix and vector
if (matM1.fIsTemp) then MatM1.Done;
if (vecV2.fIsTemp) then vecV2.pOrg^.fFreeToUse:=TRUE;
end; // operator *  matrix time vector

operator * (matM1, matM2:TMatrixRNMFloat) matM3:TMatrixRNMFloat;
//==============================================================
var
   ii,jj,kk,nkMax : int32;
   vSumL  : Float; // sums needed in matrix multiplication
begin
// matM3 = temporary matrix.
// free the matrices matM1 and matM2 if they are temporary
SetTemporaryMat(matM1,matM2,matM3);

// Multiplication of the matrices
for ii:=1 to matM2.nIndLine do begin  // Loop on lines
   // Select the number of termes in the sum
   if (matM2.nIndLine <= matM1.pavecMat^[ii].nIndMax) then
        nkMax :=matM2.nIndLine
   else nkMax :=matM1.pavecMat^[ii].nIndMax;

   // It is possible that each line of the matrix has a different length
   matM3.pavecMat^[ii].nIndMax:=matM2.pavecMat^[ii].nIndMax;

   for jj:=1 to matM2.pavecMat^[ii].nIndMax do begin  // Loop on columns
      vSumL:=0.0;

      // Line Col multiplication
      for kk:=1 to nkMax do
         vSumL:=vSumL + matM1.pavecMat^[ii].pavV^[kk] *
                        matM2.pavecMat^[kk].pavV^[jj];

      matM3.pavecMat^[ii].pavV^[jj]:=vSumL;
      end; // for jj  columns
   end; // for ii  lines

// Number of lines of the resulting matrix
matM3.nIndLine  :=matM1.nIndLine;

// destructor are called, they do little for temp matrices
if (matM1.fIsTemp) then MatM1.Done;
if (matM2.fIsTemp) then MatM2.Done;

// Indicate that this is not free to use
matM3.pOrg^.fFreeToUse:=FALSE;
end; // operator *  matrix multiplication


end. // unit xMatrixRN
