flinalgebkandnvs.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Library for Turbo Pascal: Test LinAlg)
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(09.09.2002)
@lastmod(10.10.2002)
㎞ikolai V. Shokhirev, 2002
}
unit fLinAlgEBKandNVS;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  MathTypes, SpecFunc, StdCtrls, ShellAPI, ComCtrls, LinAlg, ExtCtrls,
  MatrixProc;
type
  TForm1 = class(TForm)
    PageControl1: TPageControl;
    TabSheet3: TTabSheet;
    Button1: TButton;
    Panel1: TPanel;
    Label3: TLabel;
    Button2: TButton;
    TabSheet1: TTabSheet;
    Button3: TButton;
    Button4: TButton;
    TabSheet2: TTabSheet;
    Button6: TButton;
    TabSheet4: TTabSheet;
    Button5: TButton;
    Button7: TButton;
    procedure Label3Click(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure Button3Click(Sender: TObject);
    procedure Button4Click(Sender: TObject);
    procedure Button5Click(Sender: TObject);
    procedure Button6Click(Sender: TObject);
    procedure Button7Click(Sender: TObject);
  private
    { Private declarations }
    signal: IntType;
    m4x3, m4x4, m3x3, C4x4, C3x3,
    m3x4, Old3x3, New3x3, Old4x4, New4x4, EDE4: MtPtr;
    E4, E3, A4, A3: VcPtr;
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

uses
  ReadWrite;

{$R *.DFM}

procedure TForm1.FormCreate(Sender: TObject);
begin
  GetMatrixMemory(m4x3, 4,3);
  GetMatrixMemory(m4x4, 4,4);
  GetMatrixMemory(m3x3, 3,3);
  GetMatrixMemory(m3x4, 3,4);
  GetMatrixMemory(C4x4, 4,4);
  GetMatrixMemory(C3x3, 3,3);
  GetMatrixMemory(Old3x3, 3,3);
  GetMatrixMemory(New3x3, 3,3);
  GetMatrixMemory(Old4x4, 4,4);
  GetMatrixMemory(New4x4, 4,4);
  GetMatrixMemory(EDE4, 4,4);
  GetVectorMemory(E4, 4);
  GetVectorMemory(E3, 3);
  GetVectorMemory(A4, 4);
  GetVectorMemory(A3, 3);
end;

procedure TForm1.FormDestroy(Sender: TObject);
begin
  FreeMatrixMemory(m4x3, 4,3);
  FreeMatrixMemory(m4x4, 4,4);
  FreeMatrixMemory(m3x3, 3,3);
  FreeMatrixMemory(m3x4, 3,4);
  FreeMatrixMemory(C4x4, 4,4);
  FreeMatrixMemory(C3x3, 3,3);
  FreeMatrixMemory(Old3x3, 3,3);
  FreeMatrixMemory(New3x3, 3,3);
  FreeMatrixMemory(Old4x4, 4,4);
  FreeMatrixMemory(New4x4, 4,4);
  FreeMatrixMemory(EDE4, 4,4);
  FreeVectorMemory(E4, 4);
  FreeVectorMemory(E3, 3);
  FreeVectorMemory(A4, 4);
  FreeVectorMemory(A3, 3);
end;

{ Jacobi Diagonalization
       M*C = C*E  or  CT*M*C = E
  test
       M = C*E*CT  or M - C*E*CT = 0 }
procedure TForm1.Button1Click(Sender: TObject);
begin
  ReadMatrix(m4x3^, 4, 3, 'm4x3.txt');

//  m3x3 = Transpose(m4x3)x m4x3
  TransposeMN1N2(4, 3, m3x4^, m4x3^);
  WriteMatrix(m3x4^, 3,4, 'm3x4.txt');
// m3x3 = m3x4 * m4x3
  MN1N2xMN2N3(3, 4, 3,m3x3^, m3x4^, m4x3^);

  WriteMatrix(m3x3^, 3,3, 'm3x3.txt');
  CopyM(3, Old3x3^, m3x3^);

  Jacobi(3, m3x3^, C3x3^, E3^, A3^, Signal);
  if signal > 0 then
    ShowMessage('signal'+IntToStr(signal));
  WriteMatrix(C3x3^, 3,3, 'C3x3.txt');
  WriteVector(E3^, 3, 'E3.txt');

  // test
  MxDxMT(3, New3x3^, C3x3^, C3x3^, E3^);

  CombineM(3, m3x3^, New3x3^, Old3x3^, 1.0, -1.0);
  // Dif3x3 must be 0
  WriteMatrix(m3x3^, 3,3, 'Dif3x3.txt');
end;

{ tred2 - tqli Diagonalization
       M*C = C*E  or  CT*M*C = E
  test
       M = C*E*CT  or M - C*E*CT = 0 }
procedure TForm1.Button2Click(Sender: TObject);
var
  k: IntType;
begin
  ReadMatrix(m4x3^, 4, 3, 'm4x3.txt');

//  m4x4 = m4x3 x Transpose(m4x3)
  TransposeMN1N2(4, 3, m3x4^, m4x3^);
  WriteMatrix(m3x4^, 3,4, 'm3x4.txt');
  MN1N2xMN2N3(4, 3, 4,m4x4^, m4x3^, m3x4^);

  WriteMatrix(m4x4^, 4,4, 'm4x4.txt');
  CopyM(4, Old4x4^, m4x4^);

  tred2(m4x4^, 4, true, E4^, A4^);
  WriteMatrix(m4x4^, 4, 4, 'z4x4.txt');
  WriteVector(E4^, 4, 'diag4.txt');
  WriteVector(A4^, 4, 'offd4.txt');
  // test
  ScalarM(4, EDE4^, 0.0);  // EDE4 = 0
  EDE4^[1]^[1] := E4^[1];
  for k := 2 to 4 do
  begin
    EDE4^[k]^[k] := E4^[k];
    EDE4^[k]^[k-1] := A4^[k];
    EDE4^[k-1]^[k] := A4^[k];
  end;
  // intermediate tred2 test
  MxMxMT(4, New4x4^, EDE4^, m4x4^);
  // New4x4 = New4x4 - Old4x4
  CombineM(4, New4x4^, New4x4^, Old4x4^, 1.0, -1.0);
  // DifTred4x4 must be 0
  WriteMatrix(New4x4^, 4,4, 'DifTred4x4.txt');

  tqli(E4^, A4^, 4, true, m4x4^, signal);

  CopyM(4, C4x4^, m4x4^);
  if signal > 0 then
    ShowMessage('signal'+IntToStr(signal));
  WriteMatrix(C4x4^, 4, 4, 'C4x4.txt');
  WriteVector(E4^, 4, 'E4.txt');

  // total test
  MxDxMT(4, New4x4^, C4x4^, C4x4^, E4^);
  CombineM(4, New4x4^, New4x4^, Old4x4^, 1.0, -1.0);
  // Dif4x4 must be 0
  WriteMatrix(New4x4^, 4,4, 'Dif4x4.txt');
end;

{ SVD Test 4 x 3 test }
procedure TForm1.Button3Click(Sender: TObject);
var
  U4x3, V4x3: MtPtr;
begin
  try
    GetMatrixMemory(U4x3, 4,3);
    GetMatrixMemory(V4x3, 4,3);
    ReadMatrix(m4x3^, 4, 3, 'm4x3.txt');
    SVD(4, 4, 3, m4x3^, U4x3^, V4x3^, E4^, A4^, true, true, signal);
    WriteMatrix(U4x3^, 4,3, 'U4x3.txt');
    WriteMatrix(V4x3^, 4,3, 'V4x3.txt');
    WriteVector(E4^,4, 'W4x3.txt');
  finally
    FreeMatrixMemory(U4x3, 4,3);
    FreeMatrixMemory(V4x3, 4,3);
  end;
end;

{ SVD 3 x 4 test }
procedure TForm1.Button4Click(Sender: TObject);
var
  U4x3, V4x3: MtPtr;
begin
  try
    GetMatrixMemory(U4x3, 4,3);
    GetMatrixMemory(V4x3, 4,3);
    ReadMatrix(m4x3^, 4, 3, 'm4x3.txt');
    TransposeMN1N2(4, 3, m3x4^, m4x3^);
    SVD(3, 4, 3, m3x4^, U4x3^, V4x3^, E4^, A4^, true, true, signal);
    WriteMatrix(U4x3^, 4,3, 'U3x4.txt');
    WriteMatrix(V4x3^, 4,3, 'V3x4.txt');
    WriteVector(E4^,4, 'W3x4.txt');
  finally
    FreeMatrixMemory(U4x3, 4,3);
    FreeMatrixMemory(V4x3, 4,3);
  end;
end;

{ Gram Schmidt Orthogonalization test }
procedure TForm1.Button6Click(Sender: TObject);
var
  Nvec: IntType;
  i, k: IntType;
begin
  Nvec := 4;
  ReadMatrix(m4x3^, 4, 3, 'm4x3.txt');
  Gram_Schmidt(Nvec, 3, m4x3^);
  // Nvec = 3
  ShowMessage('Nvec = '+IntToStr(Nvec));
  WriteMatrix(m4x3^, 4,3, 'Ort_m4x3.txt');
  for i := 1 to 4 do
    for k := 1 to 4 do
      m4x4^[i]^[k] := VTxV(3, m4x3[i]^, m4x3[k]^);
  // test: must be Unit Matrix 3 x 3
  WriteMatrix(m4x4^, 4,4, 'delta4x4.txt');
end;


procedure TForm1.Button7Click(Sender: TObject);
var
  J0: IVcPtr;
  Det: RealType;
begin
  GetIVectorMemory(J0,3);
  ReadMatrix(m3x3^, 3, 3, 'm3x3.txt');
  FastInverse(3, m3x3^, J0^, Det, Signal);
  ReadMatrix(Old3x3^, 3, 3, 'm3x3.txt');
  MxM(3,C3x3^,Old3x3^,m3x3^);
  // must be 1
  WriteMatrix(C3x3^, 3,3, 'U3x3.txt');
  FreeIVectorMemory(J0,3);
end;

procedure TForm1.Label3Click(Sender: TObject);
begin
  ShellExecute(0,'open','http://www.shokhirev.com/nikolai.html','','',SW_SHOW);
end;

procedure TForm1.Button5Click(Sender: TObject);
var
  SL: TStringList;
begin
  try
    SL := TStringList.Create;
    ReadMatrix(  m4x3^, 4, 3,   'm4x3.txt');
    MatrixToHTML(m4x3^, 4, 3,SL,'m4x3.txt');
    ReadMatrix(  m4x3^, 4, 3,   'V4x3.txt');
    MatrixToHTML(m4x3^, 4, 3,SL,'V4x3.txt');
    ReadMatrix(  m4x3^, 4, 3,   'U4x3.txt');
    MatrixToHTML(m4x3^, 4, 3,SL,'U4x3.txt');
    ReadMatrix(  m3x4^, 3, 4,   'U3x4.txt');
    MatrixToHTML(m3x4^, 3, 4,SL,'U3x4.txt');
    ReadMatrix(  m3x4^, 3, 4,   'V3x4.txt');
    MatrixToHTML(m3x4^, 3, 4,SL,'V3x4.txt');
    ReadMatrix(  m3x4^, 3, 4,   'm3x4.txt');
    MatrixToHTML(m3x4^, 3, 4,SL,'m3x4.txt');
    ReadMatrix(  m3x3^, 3, 3,   'C3x3.txt');
    MatrixToHTML(m3x3^, 3, 3,SL,'C3x3.txt');
    ReadMatrix(  m3x3^, 3, 3,   'Dif3x3.txt');
    MatrixToHTML(m3x3^, 3, 3,SL,'Dif3x3.txt');
    ReadMatrix(  m3x3^, 3, 3,   'm3x3.txt');
    MatrixToHTML(m3x3^, 3, 3,SL,'m3x3.txt');
    ReadMatrix(  m4x4^, 4, 3,   'm4x4.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'m4x4.txt');
    ReadMatrix(  m4x4^, 4, 4,   'C4x4.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'C4x4.txt');
    ReadMatrix(  m4x4^, 4, 4,   'C4x4J.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'C4x4J.txt');
    ReadMatrix(  m4x4^, 4, 4,   'Dif4x4.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'Dif4x4.txt');
    ReadMatrix(  m4x4^, 4, 4,   'Dif4x4J.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'Dif4x4J.txt');
    ReadMatrix(  m4x4^, 4, 4,   'DifTred4x4.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'DifTred4x4.txt');
    ReadMatrix(  m4x4^, 4, 4,   'z4x4.txt');
    MatrixToHTML(m4x4^, 4, 4,SL,'z4x4.txt');
    ReadVector(  E3^, 3,   'E3.txt');
    VectorToHTML(E3^, 3,SL,'E3.txt');
    TVectorToHTML(E3^,3,SL,'TE3.txt');
    ReadVector(  E4^, 4,   'E4.txt');
    VectorToHTML(E4^, 4,SL,'E4.txt');
    TVectorToHTML(E4^,4,SL,'TE4.txt');
    ReadVector(  E4^, 4,   'E4J.txt');
    VectorToHTML(E4^, 4,SL,'E4J.txt');
    TVectorToHTML(E4^,4,SL,'TE4J.txt');
    ReadVector(  E4^, 4,   'offd4.txt');
    VectorToHTML(E4^, 4,SL,'offd4.txt');
    ReadVector(  E4^, 4,   'diag4.txt');
    VectorToHTML(E4^, 4,SL,'diag4.txt');
    ReadVector(  E4^, 4,   'W4x3.txt');
    VectorToHTML(E4^, 4,SL,'W4x3.txt');
    TVectorToHTML(E4^,4,SL,'TW4x3.txt');
    ReadVector(  E4^, 4,   'W3x4.txt');
    VectorToHTML(E4^, 4,SL,'W3x4.txt');
    TVectorToHTML(E4^,4,SL,'TW3x4.txt');
    SL.SaveToFile('AllArrays.txt');
  finally
    SL.Free;
  end;
end;

end.