{@abstract(Linear algebra vector and matrix procedures )
author: Nikolai Shokhirev )
created: February 02, 2001)
last modified: June 06, 2003)
ŠNikolai V. Shokhirev, 2001-2003
}
unit  uAlgUtils;

interface

uses
  Math, uMatTypes, uTensors, uArrays, uCompoundArrays;

{ V[i] := h }
procedure FillV(const V: ITensor1D; h: TFloat);

{ V1 := V2*C ; allows A = B,}
procedure VxC(const V1, V2: ITensor1D; c: TFloat);

{ swaps V[i] and V[j] }
procedure SwapV(const V: ITensor1D; i, j: TInt);

{ V1[i] := V2[i] }
procedure VAssignV(const V1, V2: ITensor1D; var err: TInt);

{ MaxAbsV = max( |V[i]| }
function MaxAbsV(const V: ITensor1D): TFloat;

{ Dot (scalar) product: VxV := (A|B) }
function VxV(const A, B: ITensor1D): TFloat;

{ InnerProd := c + Sum(i=low,high,step; A[i]*B[i]) }
function InnerProd(low, step, high: integer; const A, B: ITensor1D;
                                                       c: TFloat = 0.0): TFloat;
{ Scalar ( dot ) product with check of limits }
function VDotV(const V1, V2: IVector; var err: TInt;
                             imin: TInt = LowInt; imax: TInt = HighInt): TFloat;

{ if |B| > SqrtMinReal then A := B/|B| , allows A = B, Norm := |B| }
procedure NormalizeV(const A, B: ITensor1D; var Norm: TFloat);

{ NormV := |B| }
function NormV(N: TInt; const B: ITensor1D): TFloat;

{ Norm2V := (B|B) }
function Norm2V(const B: ITensor1D): TFloat;

{ C = rA*A + rB*B }
procedure CombineV(const A, B: ITensor1D; rA, rB: TFloat; const C: ITensor1D;
                                                                 var err: TInt);

{ MaxAbsIndexV = IndexOf(MaxAbsV), val = V[MaxAbsIndexV] }
function MaxAbsIndexV(const V: ITensor1D; var val: TFloat; imin: TInt = LowInt;
                                                     imax: TInt = HighInt):TInt;

{ V1 := V1 + V2*h }
procedure VAddV(const V1, V2: ITensor1D; h: TFloat; var err: TInt;
                                     imin: TInt = LowInt; imax: TInt = HighInt);

{ M1 := M1 + MV2*h }
procedure MAddM(const M1, M2: ITensor2D; h: TFloat; var err: TInt);

{ Matrix transposition MT := M^T - same as MtT in uTensorUtils }
procedure Transpose(const M, MT: ITensor2D; var err: TInt);

{ Vt = c1*V1 + c2*V2 + . . . ; Vt can coincide with any Vi }
procedure SumV(cc: array of TFloat; const VV: array of ITensor1D;
                                            const Vt: ITensor1D; var err: TInt);

{ Mt = c1*M1 + c2*M2 + . . . ; Mt can coincide with any Mi}
procedure SumM(cc: array of TFloat; const MM: array of ITensor2D;
                                            const Mt: ITensor2D; var err: TInt);

{ Mt := M + Disg(d) Mt = M }
procedure MShiftDiag(const M: ITensor2D; const d: TFloat;
                                            const Mt: ITensor2D; var err: TInt);
{ Mt := M + D; Mt = M - OK }
procedure MAddDiag(const M: ITensor2D; const D: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);

{ ColNorm2 = U.Col[i].GetNorm2 }
function ColNorm2(const U: ITensor2D; i: TInt):TFloat;
{ dot product: ColxCol = (U.Col[i]|U.Col[i]) }
function ColxCol(const U: ITensor2D; i, j: TInt):TFloat;
{ U.Col[i] := U.Col[i] + C*U.Col[j] }
procedure ColAddCol(const U: ITensor2D; i, j: TInt; C: TFloat);
{ U.Col[i] := U.Col[i] + f*U.Col[j] }
procedure CxCol(const U: ITensor2D; i: TInt; C: TFloat);
{ U.Col[i] := U.Col[i]/U.Col[i].norm }
procedure NormalizeCol(const U: ITensor2D; i: TInt; var err: TInt);

{ dot product: RowxCol = (U.Row[i]|U.Col[j]) }
function RowxCol(const U: ITensor2D; i, j: TInt):TFloat;
{ dot product: RowxCol = (U1.Row[i]|U2.Col[j]) }
function Row1xCol2(const U1, U2: ITensor2D; i, j: TInt):TFloat;

{ RowNorm2 = U.Row[i].GetNorm2 }
function RowNorm2(const U: ITensor2D; i: TInt):TFloat;
{ dot product: RowxRow = (U.Row[i]|U.Row[i]) }
function RowxRow(const U: ITensor2D; i, j: TInt):TFloat;
{ U.Row[i] := U.Row[i] + C*U.Row[j] }
procedure RowAddRow(const U: ITensor2D; i, j: TInt; C: TFloat);
{ U.Row[i] := U.Row[i] + f*U.Row[j] }
procedure CxRow(const U: ITensor2D; i: TInt; C: TFloat);
{ U.Row[i] := U.Row[i]/U.Row[i].norm }
procedure NormalizeRow(const U: ITensor2D; i: TInt; var err: TInt);


{ vector-matrix product: VxM = V*M = MT*V }
procedure VxM(const VL: ITensor1D; const MR: ITensor2D;
                                            const Vt: ITensor1D; var err: TInt);
{ matrix-vector product: MxV = M*V = V*MT }
procedure MxV(const ML: ITensor2D; const VR: ITensor1D;
                                            const Vt: ITensor1D; var err: TInt);

{ matrix-matrix product: Mt = ML*MR }
procedure MxM(const ML, MR: ITensor2D; const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: MtTxMt = Transpose(ML)*MR  }
procedure MTxM(const ML, MR: ITensor2D; const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: MtxMtT = ML*Transpose(MR) }
procedure MxMT(const ML, MR: ITensor2D; const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: Mt = Transpose(ML)*Transpose(MR) }
procedure MT_MT(const ML, MR: ITensor2D; const Mt: ITensor2D; var err: TInt);

{ matrix-matrix product: Mt = DiagonalMatrix(DL)*MR }
procedure DxM(const DL: ITensor1D; const MR: ITensor2D;
                                            const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: Mt = ML*DiagonalMatrix(DR) }
procedure MxD(const ML: ITensor2D; const DR: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);

{ matrix-matrix product: Mt = MxDxMT = M*DiagonalMatrix(D)*Transpose(M) }
procedure MxDxMT(const M: ITensor2D; const D: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: Mt = MxDxMT = Transpose(M)*DiagonalMatrix(D)*M }
procedure MTxDxM(const M: ITensor2D; const D: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: Mt = M1xDxM2T = M*DiagonalMatrix(D)*Transpose(M) }
procedure M1xDxM2T(const M1, M2: ITensor2D; const D: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);
{ matrix-matrix product: Mt = M1TxDxM2 = Transpose(M)*DiagonalMatrix(D)*M }
procedure M1TxDxM2(const M1, M2: ITensor2D; const D: ITensor1D;
                                            const Mt: ITensor2D; var err: TInt);


{ Old variant support }

{ (M3R + i*M3I) = (M1R + i*M1I)*(M2R + i*M2I)
   M3R = M1R*M2R - M1I*M2I;  M3I = M1R*M2I + M1I*M2R }
procedure CMxCM(const M1R, M1I, M2R, M2I, M3R, M3I: ITensor2D; var err: TInt);

{ (M3R + i*M3I) = (M1RT + i*M1IT)*(M2R + i*M2I)
   M3R = M1RT*M2R - M1IT*M2I;  M3I = M1RT*M2I + M1IT*M2R }
procedure CMHxCM(const M1R, M1I, M2R, M2I, M3R, M3I: ITensor2D; var err: TInt);

{ (M3R + i*M3I) = (M1R + i*M1I)*(M2RT - i*M2IT)
   M3R = M1R*M2RT + M1I*M2IT;  M3I = -M1R*M2IT + M1I*M2RT }
procedure CMxCMH(const M1R, M1I, M2R, M2I, M3R, M3I: ITensor2D; var err: TInt);

{ (M3R + i*M3I) = (M1RT - i*M1IT)*D*(M2R + i*M2I)
   M3R = M1RT*D*M1R + M1IT*D*M1I;  M3I = M1RT*D*M2I - M1IT*D*M1R }
procedure CMHxDxCM(const M1R, M1I: ITensor2D; const D: ITensor1D;
                                        const M3R, M3I: ITensor2D; var err: TInt);

{ (M3R + i*M3I) = (M1R + i*M1I)*D*(M2RT - i*M2IT)
   M3R = M1R*D*M1RT + M1I*D*M1IT;  M3I = -M1R*D*M1IT + M1I*D*M1RT }
procedure CMxDxCMH(const M1R, M1I: ITensor2D; const D: ITensor1D;
                                        const M3R, M3I: ITensor2D; var err: TInt);

implementation


end.