cmatrixproc.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Pascal-Delphi Math Library: Basic Vector and Matrix Complex Calculations)
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(04.04.1994)
@lastmod(04.04.2002)
This is a temporary publication (reduced variant), will be updated later
}
unit cMatrixProc;

interface

uses
   MathTypes, ComplexType;

function CVHxCV(N: IntType; const A, B: CVector): complex;
function NormalizeCV(N: IntType; var A: CVector; const B: CVector): RealType;

procedure CHConjug(N: IntType; const A, C: CMatrix);
procedure DiagonalCM(N: IntType; const D: CMatrix; const V: Vector);
procedure DiagonalCMC(N: IntType; const D: CMatrix; const CV: CVector);
procedure ScalarCM(N: IntType; const D: CMatrix; c: complex);

procedure CMxCM(N: integer; const C, A, B: CMatrix);
procedure CMHxCM(N: integer; const C, A, B: CMatrix);
procedure CMxDxCMH(N: IntType; const A, C: CMatrix; const D: Vector);
procedure CMHxDxCM(N: IntType; const A, C: CMatrix; const D: Vector);
procedure CMHxCMxCM(N: IntType; const A, M, C: CMatrix);

{-----------------------------------------------------------}

implementation

function CVHxCV(N: IntType; const A, B: CVector): complex;
{ result := HTranspose(A)*B }
var
  s: complex;
  i: IntType;
begin
  s := cmplx0;
  for i:= 1 to N do
  begin
    Csum(s,Cmul(A[i],B[i]));
  end;
  result := s;
end;{ of CVHxCV }

function NormalizeCV(N: IntType; var A: CVector; const B: CVector): RealType;
{ result := Norm(B)   }
{ A := B/Norm(B)      }
{ A can be equal to B }
var
  s: RealType;
  i: IntType;
begin
  s := 0.0;
  for i:= 1 to N do
    s := s + Cabs2(B[i]);

  s := 1.0/sqrt(s);

  for i:= 1 to N do
  begin
    A[i] := Cmul(B[i],s);
  end;
  result := sqrt(s);
end;{ of NormalizeCV }


procedure CHConjug(N: IntType; const A, C: CMatrix);
{ A = Conjug(Transpose(C)) }
{ A can be equal to C      }
var
  s: complex;
  i,  j: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := C[j]^[i];
      A[j]^[i] := conjug(C[i]^[j]);
      A[j]^[i] := conjug(s);
    end;
end;{ of CHConjug }

procedure ScalarCM(N: IntType; const D: CMatrix; c: complex);
var
  k, j: IntType;
begin
  for k:= 1 to N do
  begin
    for j:= 1 to N do D[k]^[j] := cmplx0;
    D[k]^[k] := c;
  end;
end;

procedure DiagonalCM(N: IntType; const D: CMatrix; const V: Vector);
var
  k, j: IntType;
begin
  for k:= 1 to N do
  begin
    for j:= 1 to N do D[k]^[j] := cmplx0;
    D[k]^[k] := cmplx(V[k]);
  end;
end;

procedure DiagonalCMC(N: IntType; const D: CMatrix; const CV: CVector);
var
  k, j: IntType;
begin
  for k:= 1 to N do
  begin
    for j:= 1 to N do D[k]^[j] := cmplx0;
    D[k]^[k] := CV[k];
  end;
end;

procedure CMxCM(N: integer; const C, A, B: CMatrix);
{ Product of matrices A and B: C = A*B }
var
  s                    : complex;
  i,  j,   k           : IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := cmplx0;
      for k:= 1 to N do
        Csum(s, Cmul(A[i]^[k], B[k]^[j]));
//        s := s + A[i]^[k]*B[k]^[j];
      C[i]^[j] := s;
    end;
end;{ of CMxCM }

procedure CMHxCM(N: integer; const C, A, B: CMatrix);
{ C = HConjug(A)*B }
var
  s: complex;
  i,  j,   k: IntType;
begin{ of MTxM }
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := cmplx0;
      for k:= 1 to N do
        Csum(s,Cmul(conjug(A[k]^[i]),B[k]^[j]));
//        s := s + A[k]^[i]*B[k]^[j];
      C[i]^[j] := s;
    end;
end;{ of CMHxCM }

procedure CMxDxCMH(N: IntType; const A, C: CMatrix; const D: Vector);
{  A = C*D*HConjug(C) }
{  D - real diagonal  }
var
  s: complex;
  i, k, l: IntType;
begin
  for k:= 1 to N do
    for l:= 1 to N do
    begin
      s := cmplx0;
      for i:= 1 to N do
      begin
        Csum(s,Cmul(D[i],Cmul(C[k]^[i],conjug(C[l]^[i]))));
      end;
      A[k]^[l] := s;
    end;
end;{ of CMxDxCH }

procedure CMHxDxCM(N: IntType; const A, C: CMatrix; const D: Vector);
{ A = HConjug(C)*D*C       }
{ D - real diagonal matrix }
var
  s: complex;
  i, k, j: IntType;
begin
  for i:= 1 to N do
    for j:= 1 to N do
    begin
      s := cmplx0;
      for k:= 1 to N do
      begin
        Csum(s,Cmul(D[k],Cmul(conjug(C[k]^[i]),C[k]^[j])));
      end;
      A[i]^[j] := s;
    end;
end;{ of CMHxDxCM }

procedure CMHxCMxCM(N: IntType; const A, M, C: CMatrix);
{ A = HConjug(C)*M*C   }
var
  Bp: CMtPtr;
begin
  GetCMatrixMemory(Bp, N, N);
  CMHxCM(N, Bp^, C, M);
  CMxCM(N, A, Bp^, C);
  FreeCMatrixMemory(Bp, N, N);
end;{ of CMHxCMxCM }

end.