matrixproc.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Pascal-Delphi Math Library: Basic Vector and Matrix Calculations)
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(02.02.1991)
@lastmod(04.04.2002)
This is a temporary publication (reduced variant), will be updated later
check: http://www.shokhirev.com/nikolai/programs/samplecode.html
㎞ikolai V. Shokhirev, 2002
}
unit MatrixProc;

interface

uses
   MathTypes;
{ A = B*R }
procedure VxR(N: IntType; var A: Vector; const B: Vector; R: RealType);
procedure MxR(N: IntType; const A, B: Matrix; R: RealType);
{ A = B }
procedure CopyV(N: IntType; var A: Vector;const B: Vector);
procedure CopyM(N: IntType; const A, B: Matrix);
{ C = rA*A + rB*B}
procedure CombineV(N: IntType; var C: Vector; const A, B: Vector; rA, rB: RealType);
procedure CombineM(N: IntType; const C, A, B: Matrix; rA, rB: RealType);

{ Dot product: VTxV = (A*B) = <A|B> }
function VTxV(N: IntType; const A, B: Vector): RealType;

{ NormalizeV = |B| if |B| > 0 then A := B/|B| }
function NormalizeV(N: IntType; var A: Vector; const B: Vector): RealType;
{ NormV = |B| = sqrt(<B|B>) }
function NormV(N: IntType; const B: Vector): RealType;
{ Norm2V = <B|B> }
function Norm2V(N: IntType; const B: Vector): RealType;

{ A := Transpose(B) = BT }
procedure TransposeM(N: IntType; const A, B: Matrix);
{ D := Diagonal(V) D[i,i] = V[i] }
procedure DiagonalM(N: IntType; const D: Matrix; const V: Vector);
{ D := R*UnitaryMatrix }
procedure ScalarM(N: IntType; const D: Matrix; R: RealType);

procedure MxV(N: IntType; var A: Vector; const B: Matrix; const V: Vector);
procedure VxM(N: IntType; var A: Vector; const V: Vector; const B: Matrix);

procedure MxM(N: IntType; const C, A, B: Matrix);
procedure MTxM(N: IntType; const C, A, B: Matrix);
procedure MxMT(N: IntType; const C, A, B: Matrix);

procedure MTxDxM(N: IntType; const A, B, C: Matrix; const D: Vector);
procedure MxDxMT(N: IntType; const A, B, C: Matrix; const D: Vector);
procedure MTxMxM(N: IntType; const A, M, C: Matrix);
procedure MxMxMT(N: IntType; const A, M, C: Matrix);

{ Rectangualar matrices }

procedure MN1N2xMN2N3(N1, N2, N3: IntType; const A, B, C: Matrix);
procedure TransposeMN1N2(N1, N2: IntType; const A, B: Matrix);

{ Complex calculation emulation }

procedure CMxCMr(N: IntType; s: RealType; const U, V, X, Y, A, B: Matrix);
procedure CMTxCMr(N: IntType; s: RealType; const U, V, X, Y, A, B: Matrix);
procedure CMxDxCMHr(N: IntType; const AR, AI, CR, CI: Matrix; const D: Vector);
procedure CMHxDxCMr(N: IntType; const AR, AI, CR, CI: Matrix; const D: Vector);
procedure CMHxCMxCMr(N: IntType; const AR, AI, MR, MI, CR, CI: Matrix );

{ Other utilities }

{ simple buble sort of N1 to N2 row vectors by E }
procedure SortTVec(N1, N2: IntType; var VCT: Matrix; E: Vector; ascending: boolean);
{ simple buble sort of N1 to N2 row vectors by norm }
procedure SortNornTVec(N1, N2, N: IntType; var VCT: Matrix; ascending: boolean);
{ simple buble sort of N1 to N2 column vectors by E}
procedure SortVec(N1, N2, N: IntType; var VCT: Matrix; E: Vector; ascending: boolean);

implementation

procedure VxR(N: IntType; var A: Vector;const B: Vector; R: RealType);
{ allows A = B }
var
  k: IntType;
begin
  for k:= 1 to N do
    A[k] := R*B[k];
end;{ of VxR }

procedure MxR(N: IntType; const A, B: Matrix; R: RealType);
{ allows A = B }
var
  j, i: IntType;
begin
  for j:= 1 to N do
    for i:= 1 to N do
      A[i]^[j] := R*B[i]^[j];
end;{ of MxR }

procedure CopyV(N: IntType; var A: Vector;const B: Vector);
var
  k: IntType;
begin
  for k:= 1 to N do
    A[k] := B[k];
end;

procedure CopyM(N: IntType; const A, B: Matrix);
var
  j, i: IntType;
begin
  for j:= 1 to N do
    for i:= 1 to N do
      A[i]^[j] := B[i]^[j];
end;

{ C = rA*A + rB*B}
procedure CombineV(N: IntType; var C: Vector; const A, B: Vector;
                                                   rA, rB: RealType);
var
  i: IntType;
begin
  for i:= 1 to N do
    C[i] := rA*A[i] + rB*B[i];
end;

{ C = rA*A + rB*B}
procedure CombineM(N: IntType; const C, A, B: Matrix; rA, rB: RealType);
var
  j, i: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
      C[i]^[j] := rA*A[i]^[j] + rB*B[i]^[j];
end;

function VTxV(N: IntType; const A, B: Vector): RealType;
var
  s: RealType;
  k: IntType;
begin
  s := 0.0;
  for k:= 1 to N do
    s := s + A[k]*B[k];
  result := s;
end;{ of VTxV }

function NormalizeV(N: IntType; var A: Vector; const B: Vector): RealType;
{ allows A = B }
var
  s: RealType;
  k: IntType;
begin
  s := sqrt(VTxV(N, B, B));
  if s > HalfMinReal then
    for k:= 1 to N do
      A[k] := B[k]/s;
  result := s;
end;{ of NormalizeV }

function NormV(N: IntType; const B: Vector): RealType;
begin
  result := sqrt(VTxV(N, B, B));
end;{ of NormV }

function Norm2V(N: IntType; const B: Vector): RealType;
begin
  result := VTxV(N, B, B);
end;{ of Norm2V }

procedure TransposeM(N: IntType; const A, B: Matrix);
var
  s: RealType;
  i, j: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := B[i]^[j];
      A[i]^[j] := B[j]^[i];
      A[j]^[i] := s;
    end;
end;{ of TransposeM }

procedure DiagonalM(N: IntType; const D: Matrix; const V: Vector);
var
  i, j: IntType;
begin
  for i:= 1 to N do
  begin
    for j:= 1 to N do
      D[i]^[j] := 0.0;
  D[i]^[i] := V[i];
  end;    
end;{ of MxM }

procedure ScalarM(N: IntType; const D: Matrix; R: RealType);
var
  i, j: IntType;
begin
  for i:= 1 to N do
  begin
    for j:= 1 to N do
      D[i]^[j] := 0.0;
  D[i]^[i] := R;
  end;    
end;{ of MxM }

 // allows A = V
procedure MxV(N: IntType; var A: Vector; const B: Matrix; const V: Vector);
var
  s: RealType;
  i: IntType;
  Vc: VcPtr;
begin
  GetVectorMemory(Vc, N);
//  CopyV(N, Vc^, V);
  for i:= 1 to N do
  begin
    s := VTxV(N, B[i]^,Vc^);
//    s := 0.0;
//    for k:= 1 to N do
//      s := s + B[i]^[k]*Vc^[k];
    A[i] := s;
  end;
  FreeVectorMemory(Vc, N);
end;{ of MxV }

procedure VxM(N: IntType; var A: Vector; const V: Vector; const B: Matrix);
var
  s: RealType;
  j, k: IntType;
  Vc: VcPtr;                    // allows A = V
begin
  GetVectorMemory(Vc, N);       // allows A = V 
  CopyV(N, Vc^, V);                // allows A = V
  for j:= 1 to N do
  begin
    s := 0.0;
    for k:= 1 to N do
      s := s + Vc^[k]*B[k]^[j]; // allows A = V
//      s := s + V[k]*B[k]^[j];
    A[j] := s;
  end;
  FreeVectorMemory(Vc, N);     // allows A = V
end;{ of MxM }

{ }
 
procedure MxM(N: IntType; const C, A, B: Matrix);
{ Product of matrices A and B: C = A*B  }
var
  s: RealType;
  i, j, k: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := 0.0;
      for k:= 1 to N do
        s := s + A[i]^[k]*B[k]^[j];
      C[i]^[j] := s;
    end;
end;{ of MxM }

 
procedure MTxM(N: IntType; const C, A, B: Matrix);
{ C = AT*B  }
var
  s                    : RealType;
  i,  j,   k           : IntType;
begin{ of MTxM }
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := 0.0;
      for k:= 1 to N do
        s := s + A[k]^[i]*B[k]^[j];
      C[i]^[j] := s;
    end;
end;{ of MTxM }

procedure MxMT(N: IntType; const C, A, B: Matrix);
{ C = A*BT }
var
  s                    : RealType;
  i,  j,   k           : IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := 0.0;
      for k:= 1 to N do
        s := s + A[i]^[k]*B[j]^[k];
      C[i]^[j] := s;
    end;
end;{ of MxMT }

procedure MTxDxM(N: IntType; const A, B, C: Matrix; const D: Vector);
{ A = Transpose(B)*D*C     }
{ D - real diagonal  matrix }
var
  Re                  : RealType;
  i,   j,   k         : IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      Re := 0.0;
      for k:= 1 to N do
      begin
        Re := Re + D[k]*B[k]^[i]*C[k]^[j] ;
      end;
      A[i]^[j] := Re;
    end;
end;{ of MTxDxM }

{ A = B*D*Transpose(C)      }
{ D - real diagonal  matrix }
{ B = C is allowed          }
procedure MxDxMT(N: IntType; const A, B, C: Matrix; const D: Vector);
var
  Re                  : RealType;
  i,   j,   k         : IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      Re := 0.0;
      for k:= 1 to N do
      begin
        Re := Re + D[k]*B[i]^[k]*C[j]^[k] ;
      end;
      A[i]^[j] := Re;
    end;
end;{ of MxDxMT }

procedure MTxMxM(N: IntType; const A, M, C: Matrix);
{ A = Transpose(C)*M*C   }
var
  Bp: MtPtr;
begin
  GetMatrixMemory(Bp, N, N);
  MTxM(N, Bp^, C, M);
  MxM(N, A, Bp^, C);
  FreeMatrixMemory(Bp, N, N);
end;{ of MTxMxM }

procedure MxMxMT(N: IntType; const A, M, C: Matrix);
{ A = C*M*Transpose(C)   }
var
  Bp: MtPtr;
begin
  GetMatrixMemory(Bp, N, N);
  MxM(N, Bp^, C, M);
  MxMT(N, A, Bp^, C);
  FreeMatrixMemory(Bp, N, N);
end;{ of MTxMxM }


{ Rectangualar matrices }

{   A  =  B  *  C   }
{ N1xN3 N1xN2 N2xN3 }
procedure MN1N2xMN2N3(N1, N2, N3: IntType; const A, B, C: Matrix);
var
  i1, i2, i3: IntType;
  x: RealType;
begin
  for i1 := 1 to N1 do
  begin
    for i3 := 1 to N3 do
    begin
      x := 0.0;
      for i2 := 1 to N2 do
      begin
        x := x + B[i1]^[i2]*C[i2]^[i3]
      end;
      A[i1]^[i3] := x;
    end;
  end;
end;

{   A  = Transpose( B )  }
{ N2xN1           N1xN2 }
procedure TransposeMN1N2(N1, N2: IntType; const A, B: Matrix);
var
  i1, i2: IntType;
begin
  for i1 := 1 to N1 do
  begin
    for i2 := 1 to N2 do
    begin
      A[i2]^[i1] := B[i1]^[i2];
    end;
  end;
end;


{ Complex calculation emulation }

procedure CMxCMr(N: IntType; s: RealType; const U, V, X, Y, A, B: Matrix);
{  U + i*V = (X + i*Y*s )*(A + i*B )    }
var
  Re, Im: RealType;
  i,   j,   k         : IntType;
begin

  for i:= 1 to N do
    for j:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for k:= 1 to N do
      begin
        Re := Re + X[i]^[k]*A[k]^[j] - s*Y[i]^[k]*B[k]^[j];
        Im := Im + X[i]^[k]*B[k]^[j] + s*Y[i]^[k]*A[k]^[j];
      end;
      U[i]^[j] := Re;
      V[i]^[j] := Im;
    end;
end;{ of CMxCMr }


procedure CMTxCMr(N: IntType; s: RealType; const U, V, X, Y, A, B: Matrix);
{  U + i*V = Transpose(X + i*Y*s ) *(A + i*B )   }
var
  Re, Im: RealType;
  i,   j,   k: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for k:= 1 to N do
      begin
        Re := Re + X[k]^[i]*A[k]^[j] - s*Y[k]^[i]*B[k]^[j];
        Im := Im + X[k]^[i]*B[k]^[j] + s*Y[k]^[i]*A[k]^[j];
      end;
      U[i]^[j] := Re;
      V[i]^[j] := Im;
    end;
end;{ of CMTxCMr }

procedure CMxDxCMHr(N: IntType; const AR, AI, CR, CI: Matrix; const D: Vector);
{ A = C*D*HConjugate(C)                      }
{ A = AR + i*Ai     C = CR + i*CI            }
{ D - real diagonal  matrix                  }
{ AR + i*AI = (CR+i*CI)*D*Transpose(CR-i*CI) }
var
  Re,  Im             : RealType;
  i,   k,   l         : IntType;
begin
  for k:= 1 to N do
    for l:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for i:= 1 to N do
      begin
        Re := Re + D[i]*(CR[k]^[i]*CR[l]^[i] + CI[k]^[i]*CI[l]^[i]);
        Im := Im + D[i]*(CI[k]^[i]*CR[l]^[i] - CR[k]^[i]*CI[l]^[i]);
      end;
      AR[k]^[l] := Re;
      AI[k]^[l] := Im;
    end;
end;{ of CMxDxCMHr }

procedure CMHxDxCMr(N: IntType; const AR, AI, CR, CI: Matrix; const D: Vector);
{ A = HConjugate(C) *D*C                     }
{ A = AR + i*Ai     C = CR + i*CI            }
{ D - real diagonal  matrix                  }
{ AR + i*AI = Transpose(CR-i*CI)*D*(CR+i*CI) }
var
  Re,  Im: RealType;
  i,   j,   k: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for k:= 1 to N do
      begin
        Re := Re + D[k]*(CR[k]^[i]*CR[k]^[j] + CI[k]^[i]*CI[k]^[j]);
        Im := Im + D[k]*(CR[k]^[i]*CI[k]^[j] - CI[k]^[i]*CR[k]^[j]);
      end;
      AR[i]^[j] := Re;
      AI[i]^[j] := Im;
    end;

end;{ of CMHxDxCMr }

procedure CMHxCMxCMr(N: IntType; const AR, AI, MR, MI, CR, CI: Matrix );
{ A = HConjugate(C)*M*C                      }
{ M = MR + i*Mi         C = CR + i*CI        }
{ A = AR + i*Ai                              }
{ AR + i*AI = Transpose(CR-i*CI)*M*(CR+i*CI) }
var
  Re,      Im                   : RealType;
  i,       j,      k,       l   : IntType;
  BR,      BI                   : MtPtr;
begin
  GetMatrixMemory(BR, N, N);
  GetMatrixMemory(BI, N, N);
  for k:= 1 to N do
    for i:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for j:= 1 to N do
      begin
        Re := Re + MR[i]^[j]*CR[j]^[k] - MI[i]^[j]*CI[j]^[k];
        Im := Im + MR[i]^[j]*CI[j]^[k] + MI[i]^[j]*CR[j]^[k];
      end;
      BR^[i]^[k] := Re;
      BI^[i]^[k] := Im;
    end;
  for k:= 1 to N do
    for l:= 1 to N do
    begin
      Re := 0.0;    Im := 0.0;
      for j:= 1 to N do
      begin
        Re := Re + CR[j]^[l]*BR^[j]^[k] + CI[j]^[l]*BI^[j]^[k];
        Im := Im + CR[j]^[l]*BI^[j]^[k] - CI[j]^[l]*BR^[j]^[k];
      end;
      AR[l]^[k] := Re;
      AI[l]^[k] := Im;
    end;
  FreeMatrixMemory(BI, N, N);
  FreeMatrixMemory(BR, N, N);
  
end;{ of CMHxCMxCMr }

{ Other utilities }

{ simple buble sort of N1 to N2 row vectors by E }
procedure SortTVec(N1, N2: IntType; var VCT: Matrix; E: Vector; ascending: boolean);
var
  i, j: IntType;
  vp: VcPtr;
  t: RealType;
begin
  for i := N1 to (N2-1) do
    for j := i+1 to N2 do
      if    (Ascending and (E[j] < E[i])) or
        (not Ascending and (E[j] > E[i])) then
      begin
        t := E[i];
        E[i] := E[j];
        E[j] := t;
        vp := VCT[i];
        VCT[i] := VCT[j];
        VCT[j] := vp;
      end;
end;

{ simple buble sort of N1 to N2 row vectors by norm }
procedure SortNornTVec(N1, N2, N: IntType; var VCT: Matrix; ascending: boolean);
var
  i, j: IntType;
  vp: VcPtr;
  Ei, Ej: RealType;
begin
  for i := N1 to (N2-1) do
  begin
    Ei := VTxV(N, VCT[i]^, VCT[i]^);
    for j := i+1 to N2 do
    begin
      Ej := VTxV(N, VCT[i]^, VCT[i]^);
      if    (Ascending and (Ej < Ei)) or
        (not Ascending and (Ej > Ei)) then
      begin
        vp := VCT[i];
        VCT[i] := VCT[j];
        VCT[j] := vp;
      end;
    end;
  end;
end;

{ simple buble sort of N1 to N2 column vectors by E}
procedure SortVec(N1, N2, N: IntType; var VCT: Matrix; E: Vector; ascending: boolean);
var
  i, j, k: IntType;
  t: RealType;
begin
  for i := N1 to (N2-1) do
    for j := i+1 to N2 do
      if    (Ascending and (E[j] < E[i])) or
        (not Ascending and (E[j] > E[i])) then
      begin
        t := E[i];
        E[i] := E[j];
        E[j] := t;
        for k := 1 to N do
        begin
          t := VCT[k]^[i];
          VCT[k]^[i] := VCT[k]^[j];
          VCT[k]^[j] := t;
        end;
      end;
end;

end.