{Object algebra routines
author: Nikolai Shokhirev )
created: 2002.06.06
last modified: 2004.11.11
ŠNikolai V. Shokhirev, 2002-2004 }
unit uDynObjAlg;

interface

uses
  uMatTypes, uDynArrays, uDynLinAlg;

{ Matrix transposition }
function Transpose(const M: IFArr2D): IFArr2D; overload;
{ Matrix transposition }
function Transpose(const M: ICArr2D): ICArr2D; overload;

{ SumVt = c1*V1 + c2*V2 + . . . }
function SumVt(cc: array of TFloat; const VV: array of IFArr1D): IFArr1D;
{ SumMt = c1*M1 + c2*M2 + . . . }
function SumMt(cc: array of TFloat; const MM: array of IFArr2D): IFArr2D;

{ M := M + d }
function MtShiftDiag(const M: IFArr2D; const d: TFloat): IFArr2D;
{ M := M + D }
function MtAddDiag(const M: IFArr2D; const D: IFArr1D): IFArr2D;

{ Dot (scalar) product of two vectors VtxVt = VL*VR }
function VtxVt(const VL, VR: IFArr1D): TFloat;
{ vector-matrix product: VtxMt = V*M = MT*V }
function VtxMt(const VL: IFArr1D; const MR: IFArr2D): IFArr1D;
{ matrix-vector product: MtxVt = M*V = V*MT }
function MtxVt(const ML: IFArr2D; const VR: IFArr1D): IFArr1D;

{ matrix-matrix product: MtxMt = ML*MR }
function MtxMt(const ML, MR: IFArr2D): IFArr2D;
{ matrix-matrix product: MtTxMt = Transpose(ML)*MR  }
function MtTxMt(const ML, MR: IFArr2D): IFArr2D;
{ matrix-matrix product: MtxMtT = ML*Transpose(MR) }
function MtxMtT(const ML, MR: IFArr2D): IFArr2D;
{ matrix-matrix product: MtTxMtT = Transpose(ML)*Transpose(MR) }
function MtTxMtT(const ML, MR: IFArr2D): IFArr2D;

{ matrix-matrix product: DxMt = DiagonalMatrix(DL)*MR }
function DxMt(const DL: IFArr1D; const MR: IFArr2D): IFArr2D;
{ matrix-matrix product: MtxD = ML*DiagonalMatrix(DR) }
function MtxD(const ML: IFArr2D; const DR: IFArr1D): IFArr2D;

{ matrix-matrix product: MtxDxMtT = M*DiagonalMatrix(D)*Transpose(M) }
function MtxDxMtT(const M: IFArr2D; const D: IFArr1D): IFArr2D;
{ matrix-matrix product: MtxDxMtT = Transpose(M)*DiagonalMatrix(D)*M }
function MtTxDxMt(const M: IFArr2D; const D: IFArr1D): IFArr2D;
{ matrix-matrix product: Mt1xDxMt2T = M*DiagonalMatrix(D)*Transpose(M) }
function Mt1xDxMt2T(const M1: IFArr2D; const D: IFArr1D; const M2: IFArr2D): IFArr2D;
{ matrix-matrix product: Mt1TxDxMt2 = Transpose(M)*DiagonalMatrix(D)*M }
function Mt1TxDxMt2(const M1: IFArr2D; const D: IFArr1D; const M2: IFArr2D): IFArr2D;

{ Complex Algebra routines:  }

{ complex Conjugate }
function Conjugate(const A: ICArr1D): ICArr1D; overload;
{ complex Conjugate }
function Conjugate(const A: ICArr2D): ICArr2D; overload;
{ complex Hermitian Conjugate = Transpose (Conjugate(A)) }
function HConjugate(const A: ICArr2D): ICArr2D;

{ Dot (scalar) product of two vectors CVtxCVt=VL*VR }
function CVtxCVt(const VL, VR: ICArr1D): Complex;
{ vector-Vector product: CVtHxVRt = HConjugate(VL)*VR }
function CVtHxCVt(const VL, VR: ICArr1D): Complex;

{ vector-matrix product: CVtxCMt = V*M = MT*V  }
function CVtxCMt(const VL: ICArr1D; const MR: ICArr2D): ICArr1D;
{ matrix-vector product: CMtxCVt = M*V = V*MT }
function CMtxCVt(const ML: ICArr2D; const VR: ICArr1D): ICArr1D;

{ matrix-matrix product: CMtxCMt = ML*MR }
function CMtxCMt(const ML, MR: ICArr2D): ICArr2D;

{ matrix-matrix product: DxCMt = DiagonalMatrix(DL)*MR }
function DxCMt(const DL: IFArr1D; const MR: ICArr2D): ICArr2D;
{ matrix-matrix product: CMtxD = ML*DiagonalMatrix(DR) }
function CMtxD(const ML: ICArr2D; const DR: IFArr1D): ICArr2D;

{ Diagonalization of the real symmetric matrix a. }
function DiagMt(const a: IFArr2D; var code: TInt;
  sort: boolean = true; vectors: boolean = true; Ascending: boolean = true;
  Epsilon : TFloat = SqrtMinFloat; IterMax: TInt = 999): IEigenSys;

{ result := MtxDxMtT(ES.Vectors,ES.Values)  }
function EigenSysToMt(const ES: IEigenSys): IFArr2D;

{ SVD for an arbitrary matrix a:
  A[aLo1..aHi1, aLo2..aHi2],  m = A.Dim1, n = A.Dim2, nm = min(n,m)
  fullmat = false (default)
    A   =  U    *   Q    *  VT        UT  *  A  *   V    =   Q
  m x n  m x nm  nm x nm  nm x n   nm x m  m x n  n x nm  nm x nm
  U[aLo1..aHi1, 1..nm], V[aLo2..aHi2, 1..nm2], Q[1..nm]
  fullmat = true
    A   =  U  *   Q   *  VT      UT  *  A  *   V   =  Q
  m x n  m x m  m x n  n x n   m x m  m x n  n x n  m x n
  U[aLo1..aHi1, 1..m], V[aLo2..aHi2, 1..n]
  Q[1, nm] - first nm values                                       }
function SVDMt(const a: IFArr2D; fullmat: boolean = false;
               eps : TFloat = cMachEps; tol : TFloat = MinFloat): ISVDSys;

{ Inverse to SVDMt: A = U * Q * VT  }
function SVDSysToMt(const SVDS: ISVDSys): IFArr2D;

{ author: Nikolai V. Shokhirev 2002
  Diagonalization of complex Hermitian matrix }
function DiagHMt(const a: ICArr2D; var code: TInt; sort: boolean = true;
                vectors: boolean = true; Ascending: boolean = true): IHEigenSys;

{ result := MtxDxMtH(HES.Vectors,HES.Values)  }
function HEigenSysToMt(const HES: IHEigenSys): ICArr2D;

{ author: Nikolai V. Shokhirev 2004
  Diagonalization of complex Hermitian matrix
 Moore-Penrose Matrix Inverse (Pseudoinverse)
Input:
  a - m x n arbitrary rectangular matrix; mn = min(m,n)
  SVD: a = Sum(i=1,nm)[wi*|ui><vi|]
Pseudoinverse:
  a~ = Sum(j=1,k)[1/(wi+s)*|ui><vi|]  wi > e;  e, s >= 0;  k <= mn

  Moore-Penrose Inverse takes place when  e = s = 0.
  Properties: a*(a~*a)=a  (a*a~)*a = a
  For non-singular matrix a~ = a^(-1)  (or a**(-1) )

  The parameters e, s > 0;  k < mn provide regularization
                                           (stabilization)          }
function PseudoinverseMt(const a: IFArr2D; k: TInt = _; e : TFloat = MinFloat;
                                                    s : TFloat = 0.0): IFArr2D;
implementation

uses
  SysUtils, Math, uDynArrUtils;

end.

Back

Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)