{ Testing  xMatrixRN.pp
  See comments in this unit.
  --------------------------------------------}
program tEx061;

// To tell FreePascal that the application is a Graphical User Interface
// i.e. is running in windows, not in text mode.
{$APPTYPE GUI}

Uses

  wApplication0
 ,wWindow0
 ,xMatrixRN
 ,xVectRN
 ,wGlobals2
 ,wDebug0
 ,Windows;

type
  PMyWindow = ^TMyWindow;
  TMyWindow = object(TWindow)
  //"""""""""""""""""""""""""
  matM1, matM2, matM3, matM4, matM5 : TMatrixRNMFloat;
  // Matrix to use for the tests.

  vecV1,vecV2,vecV3 : RVectRNFloat;
  // Vectors to use for the test.

  constructor Init(strWinTitlePrm:ShortString);
  destructor Done; virtual;

  procedure Test1;
  procedure Test2;
  procedure Test3;
  procedure Test4;
  procedure Test5;

  function WMCommand(wPar:WPARAM; lPar:LPARAM):LRESULT; virtual;
  // Handles the menu messages

  function WindowProc(hWindow:HWND; wMsg:UINT; wPar:WPARAM; lPar:LPARAM):LRESULT; virtual;
  end; // TMyWindow = object(TWindow)

var
  oApp        : TApplication;
  poMyMainWin : PMyWindow;
  glf         : boolean; // necessary to use a special assignment operator

constructor TMyWindow.Init(strWinTitlePrm:ShortString);
//=====================================================
begin
inherited Init(strWinTitlePrm);

VNCreateNewVecRNFamilly(9,   75); // 9 = Familly, 75=maximum number of columns
//         ---
VNCreateNewMatRNMFamilly(9,  70); // 9 = Familly, 70=maximum number of lines

// Must initialise the matrices.  9 = familly
matM1.Init(9);  matM2.Init(9);  matM3.Init(9);  matM4.Init(9);  matM5.Init(9);

// Number of lines and columns used are still not defined.

VNInit(vecV1, 9);  VNInit(vecV2, 9);  VNInit(vecV3, 9);

// This line is all what is needed to have a status line.
new(poStatusLine, Init(@Self));
end;  // TMyWindow.Init

destructor TMyWindow.Done;
//=========================
begin
VNDone(vecV1);  VNDone(vecV2);  VNDone(vecV3);

// Free memory. First for the matrices, next for the familly.
matM5.Done;  matM4.Done;  matM3.Done;  matM2.Done;  matM1.Done;
VNDisposeMatRNMFamilly(9);
VNDisposeVecRNFamilly(9); // 9 = Familly

inherited Done;
end;  // TMyWindow.Done

procedure TMyWindow.Test1;
//========================
var
  matML1, matML2, matML3 : TMatrixRNMFloat;
  ii : int32; // to loop on lines and columns
begin
VNCreateNewVecRNFamilly(1,   10); // 1 = Familly, 10=maximum number of columns
VNCreateNewMatRNMFamilly(1,    8); // 1 = Familly, 8=maximum number of lines

// Must initialise the matrices
matML1.Init(1);  matML2.Init(1);  matML3.Init(1);

// Defines the number of lines of the matrix  matML1  and  matML2
matML1.nIndLine:=3;   matML1.SetColumns(3); // Number of lines and columns
matML2.nIndLine:=3;   matML2.SetColumns(3); // Number of lines and columns

// Defining matrix matML1
VNAssign(matML1.pavecMat^[1], 1,3, 1,2,3,0,0); // first line
VNAssign(matML1.pavecMat^[2], 1,3, 2,3,4,0,0); // second line
VNAssign(matML1.pavecMat^[3], 1,3, 3,4,5,0,0); // third line

// Defining matrix matML2
for ii:=1 to 3 do begin // loop on the lines
   VNAssign(matML2.pavecMat^[ii], 1,3, 9-ii, 8-ii, 7-ii, 0,0);
   end;

// Adding both matrices. This defines the number of lines and columns of matML3
glf:=matML3<=matML1 + matML2;

// displaying the resulting matrix
for ii:=1 to 3 do begin // loop on the lines
   DBW(1, 'Mat3 = ' + VNVectToStr(matML3.pavecMat^[ii], 1,3));
   end;

// Just to see if this makes problems. -> No
glf:=matML3<=matML1 + matML2;
glf:=matML3<=matML1 + matML2;

// Adding a matrix with itself works too.
glf:=matML3<=matML3 + matML2;
for ii:=1 to 3 do begin // loop on the lines
   DBW(1, 'Mat3 = ' + VNVectToStr(matML3.pavecMat^[ii], 1,3));
   end;

// Free memory. First for the matrices, next for the familly.
matML3.Done;  matML2.Done;  matML1.Done;
VNDisposeMatRNMFamilly(1);
VNDisposeVecRNFamilly(1); // 1 = Familly
end;  // TMyWindow.Test1

procedure TMyWindow.Test2;
//========================
// Test the transposition of a matrix.
// Then the multiplication
var ii : integer;
begin
// Defines the number of lines of the matrix  matM1
matM1.nIndLine:=2;   matM1.SetColumns(3); // Number of lines and columns

// Defining matrix matM1
VNAssign(matM1.pavecMat^[1], 1,3, 1,2,3,0,0); // first line
VNAssign(matM1.pavecMat^[2], 1,3, 4,5,6,0,0); // second line

matM1.Transpose;

// After transposing, display the matrix.
DBW(1, 'Mat1, line,Col = 3, 2 =' + ITS5(matM1.nIndLine) +
                                   ITS5(matM1.pavecMat^[1].nIndMax));
for ii:=1 to 3 do begin // loop on the lines
   DBW(1, 'Mat1 = ' + VNVectToStr(matM1.pavecMat^[ii], 1,2));
   end;

// The number of lines and colums do no match, but are completed with zeros
// Thus return a 3 by 2 matrix
glf:=matM2<=matM1*matM1;
DBW(1, 'Mat2, line,Col = 3, 2 =' + ITS5(matM1.nIndLine) +
                                   ITS5(matM1.pavecMat^[1].nIndMax));
for ii:=1 to 3 do begin // loop on the lines
   DBW(1, 'Mat2 = ' + VNVectToStr(matM2.pavecMat^[ii], 1,2));
   end;

// Retranspose matM1
matM1.Transpose;

// The number of lines and colums do no match, but are completed with zeros
// Thus return a 2 by 3 matrix
glf:=matM2<=matM1*matM1;
DBW(1, 'Mat2, line,Col = 2, 3 =' + ITS5(matM1.nIndLine) +
                                   ITS5(matM1.pavecMat^[1].nIndMax));
for ii:=1 to 2 do begin // loop on the lines
   DBW(1, 'Mat2 = ' + VNVectToStr(matM2.pavecMat^[ii], 1,3));
   end;
end;  // TMyWindow.Test2

procedure TMyWindow.Test3;
//========================
// Test the  LR  decomposition and solving equations
var
   ii         : integer;
   nErr       : int32;
begin
// Defining matrix matM1. Dimensions are set after
VNAssign(matM1.pavecMat^[1], 1,5, 1000/1.0, 1000/2.0, 1000/3.0, 1000/4.0, 1000/5.0); // first line
VNAssign(matM1.pavecMat^[2], 1,5, 1000/2.0, 1000/3.0, 1000/4.0, 1000/5.0, 1000/6.0);
VNAssign(matM1.pavecMat^[3], 1,5, 1000/3.0, 1000/4.0, 1000/5.0, 1000/6.0, 1000/7.0);
VNAssign(matM1.pavecMat^[4], 1,5, 1000/4.0, 1000/5.0, 1000/6.0, 1000/7.0, 1000/8.0);
VNAssign(matM1.pavecMat^[5], 1,5, 1000/5.0, 1000/6.0, 1000/7.0, 1000/8.0, 1000/9.0);

// This only fixe values, but does not set the dimension
VNAssign(vecV1, 1,5, 3,5,7,9,11);

// Defines the number of lines of the matrix  matM1
matM1.nIndLine:=5;   matM1.SetColumns(5); // Number of lines and columns
vecV1.nIndMax:=matM1.pavecMat^[1].nIndMax; // =number of columns

// We will use vecV2 to solve the equation matM1 * X = vecV2
glf:=vecV2<=matM1 * vecV1;  // Dim = number of lines

// Less or more unknown than equations
matM1.nIndLine:=4;   matM1.SetColumns(5); // Number of lines and columns
vecV2.nIndMax :=matM1.nIndLine;

// Memorise the matrix
glf:=matM2<=matM1;
matM1.LRDecompose(0);

DBW(1, 'Permuttations :');
for ii:=1 to matM1.precInfo^.nIndexLinePermut do
   DBW(1, ITS5(ii) + ITS5(matM1.precInfo^.panLinePermut^[ii]));

DBW(1, 'After LR decomp : Rank =' + ITS5(matM1.precInfo^.nMatRank));
for ii:=1 to matM1.nIndLine do begin // loop on the lines
   DBW(1, 'Mat1 = ' + VNVectToStr(matM1.pavecMat^[ii], 1, matM1.pavecMat^[1].nIndMax));
   end;

// Solve  MatM1 * X = vecV2,  X=vecV3 as a result
matM1.SolvLR(vecV2, vecV3, nErr);

DBW(1, 'After resolution, nErr=' + ITS5(nErr));
DBW(1, 'sol=' + VNVectToStr(vecV3,1, vecV3.nIndMax));

// Test the solution
glf:=vecV1<=matM2 * vecV3;
DBW(1, 'V2=' + VNVectToStr(vecV2,1, vecV2.nIndMax));
DBW(1, 'V1=' + VNVectToStr(vecV1,1, vecV1.nIndMax));
end;  // TMyWindow.Test3

procedure TMyWindow.Test4;
//========================
// Test the  LR  decomposition and solving equations
var
   ii,jj      : integer;
   nErr       : int32;
begin
// Defines the number of lines of the matrix  matM1
matM1.nIndLine:=16;   matM1.SetColumns(13); // Number of lines and columns
vecV1.nIndMax:=matM1.pavecMat^[1].nIndMax; // =number of columns

// Defining matrix matM1
// Theoretically it is not singular, but numericaly it is.
for ii:=1 to matM1.nIndLine do begin
   for jj:=1 to matM1.pavecMat^[1].nIndMax do begin
      matM1.pavecMat^[ii].pavV^[jj]:= 1.0e3 / (ii+jj-1);
      end;
   end;

// This only fixe values, but does not set the dimension
for jj:=1 to vecV1.nIndMax do vecV1.pavV^[jj]:=(2*jj+1);

// We will use vecV2 to solve the equation matM1 * X = vecV2
glf:=vecV2<=matM1 * vecV1;  // Dim = number of lines

vecV2.pavV^[10]:=vecV2.pavV^[10] + 1e-3;
DBW(1, 'V2=' + VNVectToStr(vecV2,1, vecV2.nIndMax));

// Less or more unknown than equations
//matM1.nIndLine:=matM1.nIndLine-0;
//matM1.SetColumns(vecV1.nIndMax-0);
//vecV2.nIndMax :=matM1.nIndLine;

// If more precision is wanted
//gglvVNEpsilon:=1e-15;  // Or multiply all number by a big factor (1000 for example)

// Memorise the matrix
glf:=matM2<=matM1;
matM1.LRDecompose(0);

DBW(1, 'Permuttations :');
for ii:=1 to matM1.precInfo^.nIndexLinePermut do
   DBW(1, ITS5(ii) + ITS5(matM1.precInfo^.panLinePermut^[ii]));

DBW(1, 'After LR decomp : Rank =' + ITS5(matM1.precInfo^.nMatRank));
for ii:=1 to matM1.nIndLine do begin // loop on the lines
   DBW(1, 'Mat1 = ' + VNVectToStr(matM1.pavecMat^[ii], 1, matM1.pavecMat^[1].nIndMax));
   end;

// Solve  MatM1 * X = vecV2,  X=vecV3 as a result
matM1.SolvLR(vecV2, vecV3, nErr);

DBW(1, 'After resolution, nErr=' + ITS5(nErr));
DBW(1, 'sol=' + VNVectToStr(vecV3,1, vecV3.nIndMax));

// Test the solution
glf:=vecV1<=matM2 * vecV3;
DBW(1, 'V2=' + VNVectToStr(vecV2,1, vecV2.nIndMax));
DBW(1, 'V1=' + VNVectToStr(vecV1,1, vecV1.nIndMax));
end;  // TMyWindow.Test4

procedure TMyWindow.Test5;
//========================
// Test the  LR  decomposition in the special case of having
// more columns than lines, but a rank less than the number of lines
var
   ii,jj      : integer;
   nErr       : int32;
begin
// Defines the number of lines of the matrix  matM1
matM1.nIndLine:=6;   matM1.SetColumns(10); // Number of lines and columns
vecV1.nIndMax:=matM1.pavecMat^[1].nIndMax; // =number of columns

// Defining matrix matM1
// Theoretically it is not singular, but numericaly it is.
for ii:=1 to matM1.nIndLine do begin
   for jj:=1 to matM1.pavecMat^[1].nIndMax do begin
      matM1.pavecMat^[ii].pavV^[jj]:= 1.0e3 / (ii+jj-1);
      end;
   end;

// Line 6 = line 4 + line 5  =>  Rank = 5
glf:=matM1.pavecMat^[6]<= matM1.pavecMat^[4] + matM1.pavecMat^[5];

// We will use vecV2 to solve the equation matM1 * X = vecV2
//glf:=vecV2<=matM1 * vecV1;  // Dim = number of lines

// Memorise the matrix
glf:=matM2<=matM1;
matM1.LRDecompose(1);

DBW(1, 'Permuttations :');
for ii:=1 to matM1.precInfo^.nIndexLinePermut do
   DBW(1, ITS5(ii) + ITS5(matM1.precInfo^.panLinePermut^[ii]));

DBW(1, 'After LR decomp : Rank =' + ITS5(matM1.precInfo^.nMatRank));
for ii:=1 to matM1.nIndLine do begin // loop on the lines
   DBW(1, 'Mat1 = ' + VNVectToStr(matM1.pavecMat^[ii], 1, matM1.pavecMat^[1].nIndMax));
   end;

end;  // TMyWindow.Test5

function TMyWindow.WMCommand(wPar:WPARAM; lPar:LPARAM):LRESULT;
//=============================================================
// Handles the menu messages
Var wNrMenu     : word;
begin
WMCommand := 0; // Default return value when processing a message

if (lPar <> 0) then begin
   // call the parent WMCommand method.
   WMCommand:=inherited WMCommand(wPar, lPar);
   exit;
   end;

// for the menus,  lPar = 0
//=========================
wNrMenu   := LoWord(wPar);
Case wNrMenu Of
   101 : Test1;
   102 : Test2;
   103 : Test3;
   104 : Test4;
   105 : Test5;

   else WMCommand:=inherited WMCommand(wPar, lPar);
   end;  // Case NrMenu Of
end;  // TMyWindow.WMCommand

function TMyWindow.WindowProc(hWindow:HWND; wMsg:UINT; wPar:WPARAM; lPar:LPARAM):LRESULT;
//======================================================================================
begin
// Default treatment of a message. Does nearly nothing
WindowProc := inherited WindowProc(hWindow, wMsg, wPar, lPar);
end;  // TMyWindow.WindowProc

Begin // MAIN ******
// Name under which the main window is registered. Spaces are not allowed
oApp.Init('Example061');

// Creation of the main window
new(poMyMainWin, Init('tEx061. Testing xMatrixRN.pp.'));

oApp.Run;
oApp.Done;  // The main window object is disposed here.

glf:=glf; // To avoid a Hint
End.

