{Linear algebra utilities and basic routines
author: Nikolai Shokhirev
http://www.shokhirev.com/nikolai.html,
created: (February 02, 2001)
last modified: June 06, 2004)
©Nikolai V. Shokhirev, 2001-2004 }
unit uDynArrUtils;
interface
uses
uMatTypes, uDynArrays;
{ Replaces the last positions of Prefix with Lo1, Lo1+1,.., Hi1 }
procedure SetStringArr(const A: ISArr1D; Prefix: string; Truncate: boolean = false);
{ A := Diag(c) }
procedure AssignScalar(const A: IFArr2D; c: TFloat);
{ A := B
If not MatchLim :
|* * * * * *| |* * * * * *|
|* * * * * *| |x x x x x x| |* *|x x x x|x x|
A = |* * * * * *| B = |x x x x x x| A := B |* *|x x x x|x x|
|* * * * * *| |x x x x x x| |* *|x x x x|x x|
|* * * * * *| |* * * * * *| }
procedure AssignMatrix(const A, B: IFArr2D; MatchLim: boolean);
{ AL.Col[c1] := AR.Col[c2] }
procedure CopyCC(const AL: IFArr2D; c1: TInt; AR: IFArr2D; c2: TInt);
{ AL.Row[r1] := AR.Row[r2] }
procedure CopyRR(const AL: IFArr2D; r1: TInt; AR: IFArr2D; r2: TInt);
{ random boolean vector }
procedure RandArr1D(const V: IBArr1D); overload;
{ fills with r cMin <= r <= cMax }
procedure RandArr1D(const V: IIArr1D; cMin, cMax: TInt); overload;
{ random vector: cMin <= V[i] <= cMax }
procedure RandArr1D(const V: IFArr1D; cMin, cMax: TFloat); overload;
{ random vector: cMin <= V.Re[i] <= cMax, cMin <= V.Im[i] <= cMax }
procedure RandArr1D(const V: ICArr1D; cMin, cMax: TFloat); overload
{<pre>random boolean matrix <br>
if sym = 1 then M[i,j] = M[j,i] <br>
if sym <> 1 - no symmetry </pre>}
procedure RandArr2D(const M: IBArr2D; sym: TInt = 0); overload;
{<pre>fills with r cMin <= r <= cMax <br>
if sym = -1 then M[i,j] = - M[j,i] <br>
if sym = 1 then M[i,j] = M[j,i] <br>
if sym = 0 - no symmetry </pre>}
procedure RandArr2D(const M: IIArr2D; cMin, cMax: TInt; sym: TInt = 0); overload;
{<pre>fills with r cMin <= r < cMax <br>
if sym = -1 then M[i,j] = - M[j,i] <br>
if sym = 1 then M[i,j] = M[j,i] <br>
if sym = 0 - no symmetry </pre>}
procedure RandArr2D(const M: IFArr2D; cMin, cMax: TFloat; sym: TInt = 0); overload;
{ fills with r cMin <= r < cMax }
procedure RandArr3D(const M: IIArr3D; cMin, cMax: TInt); overload;
{ fills with r cMin <= r < cMax }
procedure RandArr3D(const M: IFArr3D; cMin, cMax: TFloat); overload;
{ random boolean vector - Same as RandArr1D, backward compatibility}
procedure RandBVt(const V: IBArr1D);
{ random boolean matrix - Same as RandArr2D, backward compatibility}
procedure RandBMt(const M: IBArr2D);
{ fills with r iMin <= r <= iMax - Same as RandArr1D, backward compatibility}
procedure RandIVt(iMin, iMax: TInt; const V: IIArr1D);
{<pre>fills with r iMin <= r <= iMax <br>
if sym = -1 then M[i,j] = - M[j,i] <br>
if sym = 1 then M[i,j] = M[j,i] <br>
if sym = 0 - no symmetry - Same as RandArr2D, backward compatibility</pre>}
procedure RandIMt(iMin, iMax: TInt; const M: IIArr2D; sym: TInt);
{ random vector: cMin <= V[i] <= cMax - Same as RandArr1D, backward compatibility}
procedure RandVt(cMin, cMax: TFloat; const V: IFArr1D);
{<pre>fills with r cMin <= r < cMax <br>
if sym = -1 then M[i,j] = - M[j,i] <br>
if sym = 1 then M[i,j] = M[j,i] <br>
if sym = 0 - no symmetry - Same as RandArr2D, backward compatibility </pre>}
procedure RandMt(cMin, cMax: TFloat; const M: IFArr2D; sym: TInt);
{ fills with r cMin <= r < cMax - Same as RandArr3D, backward compatibility}
procedure Rand3Mt(cMin, cMax: TFloat; const M: IFarr3D);
{ if |B| > SqrtMinReal then A := B/|B| , allows A = B, Norm := |B| }
procedure NormalizeV(const A, B: IFArr1D; var Norm: TFloat);
{ NormV := |B| }
function NormV(const B: IFArr1D; Normalize: boolean = false): TFloat;
{ NormV := |B| }
function NormCV(const B: ICArr1D; Normalize: boolean = false): TFloat;
{ Norm2V := (B|B) }
function Norm2V(const B: IFArr1D): TFloat;
{ Norm1MCol := Max(Norm1(M.Col[j])) }
function Norm1MCol(const M: IFArr2D): TFloat;
{ Norm1MCol := Max(Norm1(M.Row[i])) }
function Norm1MRow(const M: IFArr2D): TFloat;
{ A = c1*A1 + c2*A2 + . . . ; A can coincide with any Vi }
procedure SumArr1D(cc: array of TFloat; const AA: array of IFArr1D; const A: IFArr1D);
{ A = c1*A1 + c2*A2 + . . . ; A can coincide with any Vi }
procedure SumArr2D(cc: array of TFloat; const AA: array of IFArr2D; const A: IFArr2D);
{ C = rA*A + rB*B C = A or C = B is OK}
procedure CombineArr(const A, B: IFArr1D; rA, rB: TFloat; const C: IFArr1D); overload;
{ C = rA*A + rB*B C = A or C = B is OK}
procedure CombineArr(const A, B: IFArr2D; rA, rB: TFloat; const C: IFArr2D); overload;
{ MaxAbsV = max( |V[i]| }
function MaxAbsV(const V: IFArr1D): TFloat;
{ MaxAbsIndexV = IndexOf(MaxAbsV), val = V[MaxAbsIndexV] }
function MaxAbsIndexV(const V: IFArr1D; var val: TFloat):TInt; overload;
{ MaxAbsIndexV = IndexOf(MaxAbsV), val = V[MaxAbsIndexV] }
function MaxAbsIndexV(const V: IFArr1D; var val: TFloat; imin, imax: TInt):TInt; overload;
{ MaxAbsCol = IndexOf(Max(Abs(A.Col[c]))), val = Abs(Max(A.Col[c])) }
function MaxAbsIndexCol(const A: IFArr2D; var val: TFloat; c, imin, imax: TInt):TInt; overload;
{ MaxAbsCol = IndexOf(Max(Abs(A.Col[c]))), val = Abs(Max(A.Col[c])) }
function MaxAbsIndexCol(const A: IFArr2D; var val: TFloat; c: TInt):TInt; overload;
{ MaxAbsCol = IndexOf(Max(Abs(A.Row[r]))), val = Abs(Max(A.Row[r])) }
function MaxAbsIndexRow(const A: IFArr2D; var val: TFloat; r, imin, imax: TInt):TInt; overload;
{ MaxAbsCol = IndexOf(Max(Abs(A.Row[r]))), val = Abs(Max(A.Row[r])) }
function MaxAbsIndexRow(const A: IFArr2D; var val: TFloat; r: TInt):TInt; overload;
{ MaxAbsCol = A.Col[c].MaxAbs }
function MaxAbsCol(const A: IFArr2D; c: TInt):TFloat;
{ MaxAbsCol = A.Row[r].MaxAbs }
function MaxAbsRow(const A: IFArr2D; r: TInt):TFloat;
{ V1 := V1 + V2*h }
procedure VAddV(const V1, V2: IFArr1D; h: TFloat; imin, imax: TInt); overload;
{ V1 := V1 + V2*h }
procedure VAddV(const V1, V2: IFArr1D; h: TFloat); overload;
{ V1[i] := V2[i] + h }
procedure ShiftV(const V1, V2: IFArr1D; var err: TInt; h: TFloat; imin, imax: TInt); overload;
{ V1[i] := V2[i] + h }
procedure ShiftV(const V1, V2: IFArr1D; var err: TInt; h: TFloat); overload;
{ VDotV := (V1|V2) }
function VDotV(const V1, V2: IFArr1D; imin, imax: TInt): TFloat; overload;
{ VDotV := (V1|V2) }
function VDotV(const V1, V2: IFArr1D): TFloat; overload;
{ InnerProd := c + Sum(i=low,high,step; A[i]*B[i]) }
function InnerProd(low, step, high: integer; const A, B: IFArr1D; c: TFloat = 0.0): TFloat;
{ M1 := M1 + M2*h }
procedure MAddM(const M1, M2: IFArr2D; h: TFloat);
{ Mt := M + Diag(d); Mt = M - OK }
procedure MShiftDiag(const M: IFArr2D; const d: TFloat; const Mt: IFArr2D;
RequireSquare: boolean = true);
{ Mt := M + D; Mt = M - OK }
procedure MAddDiag(const M: IFArr2D; const D: IFArr1D; const Mt: IFArr2D);
{ Norm2Col = U.Col[i].GetNorm2 }
function Norm2Col(const U: IFArr2D; i: TInt):TFloat;
{ dot product: ColxCol = (U.Col[i]|U.Col[i]) }
function ColxCol(const U: IFArr2D; i, j: TInt):TFloat;
{ U.Col[i] := U.Col[i] + C*U.Col[j] }
procedure ColAddCol(const U: IFArr2D; i, j: TInt; C: TFloat);
{ U.Col[i] := U.Col[i] + f*U.Col[j] }
procedure CxCol(const U: IFArr2D; i: TInt; C: TFloat);
{ U.Col[i] := U.Col[i]/U.Col[i].norm }
procedure NormalizeCol(const U: IFArr2D; i: TInt; var err: TInt);
{ |U.Col[i]| if Normalize then U.Col[i] := U.Row[i]/U.Row[i].norm }
function NormCol(const U: IFArr2D; i: TInt; Normalize: boolean = false): TFloat;
// M.Col[j].AssignValues(A); (regardless of V.Base)
procedure ValuesToCol(const M: IFArr2D; j: TInt; const V: IFArr1D);
// V.AssignValues(M.Col[j]); (regardless of V.Base)
procedure ColValuesToV(const M: IFArr2D; j: TInt; const V: IFArr1D);
{ dot product: RowxCol = (U.Row[i]|U.Col[j]) }
function RowxCol(const U: IFArr2D; i, j: TInt):TFloat;
{ dot product: RowxCol = (U1.Row[i]|U2.Col[j]) }
function Row1xCol2(const U1, U2: IFArr2D; i, j: TInt):TFloat;
{ dot product: VxCol = (V |U.Col[j]) }
function VxCol(const V: IFArr1D; const U: IFArr2D; j: TInt):TFloat;
{ dot product: VxCol = (V |U.Row[j]) }
function VxRow(const V: IFArr1D; const U: IFArr2D; j: TInt):TFloat;
{ Norm2Row = U.Row[i].GetNorm2 }
function Norm2Row(const U: IFArr2D; i: TInt):TFloat;
{ dot product: RowxRow = (U.Row[i]|U.Row[i]) }
function RowxRow(const U: IFArr2D; i, j: TInt):TFloat;
{ U.Row[i] := U.Row[i] + C*U.Row[j] }
procedure RowAddRow(const U: IFArr2D; i, j: TInt; C: TFloat); overload;
{ U.Row[i] := U.Row[i] + C*U.Row[j] }
procedure RowAddRow(const U: IFArr2D; i, j: TInt; C: TFloat; imin, imax: TInt); overload;
{ U.Row[i] := U.Row[i] + f*U.Row[j] }
procedure CxRow(const U: IFArr2D; i: TInt; C: TFloat); overload;
{ U.Row[i] := U.Row[i] + f*U.Row[j] }
procedure CxRow(const U: IFArr2D; i: TInt; C: TFloat; imin, imax: TInt); overload;
{ U.Row[i] := U.Row[i]/U.Row[i].norm }
procedure NormalizeRow(const U: IFArr2D; i: TInt; var err: TInt);
{ |U.Row[i]| if Normalize then U.Row[i] := U.Row[i]/U.Row[i].norm }
function NormRow(const U: IFArr2D; i: TInt; Normalize: boolean = false): TFloat;
// M.Row[j].AssignValues(A); (regardless of A.Base)
procedure ValuesToRow(const M: IFArr2D; j: TInt; const A: IFArr1D);
// V.AssignValues(M.Row[j]); (regardless of V.Base)
procedure RowValuesToV(const M: IFArr2D; j: TInt; const V: IFArr1D);
{ VL := VR*C}
procedure VxC(const VL, VR: IFArr1D; C: TFloat);
{ Dot (scalar) product of two vectors VtxVt = VL*VR without limit checks}
function VxV(const VL, VR: IFArr1D): TFloat;
{ vector-matrix product: VxM = V*M = MT*V }
procedure VxM(const VL: IFArr1D; const MR: IFArr2D; const Vt: IFArr1D);
{ matrix-vector product: MxV = M*V = V*MT }
procedure MxV(const ML: IFArr2D; const VR: IFArr1D; const Vt: IFArr1D);
{ matrix-matrix product: Mt = ML*MR }
procedure MxM(const ML, MR: IFArr2D; const Mt: IFArr2D);
{ matrix-matrix product: MtTxMt = Transpose(ML)*MR }
procedure MTxM(const ML, MR: IFArr2D; const Mt: IFArr2D);
{ matrix-matrix product: MtxMtT = ML*Transpose(MR) }
procedure MxMT(const ML, MR: IFArr2D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = Transpose(ML)*Transpose(MR) }
procedure MT_MT(const ML, MR: IFArr2D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = DiagonalMatrix(DL)*MR }
procedure DxM(const DL: IFArr1D; const MR: IFArr2D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = ML*DiagonalMatrix(DR) }
procedure MxD(const ML: IFArr2D; const DR: IFArr1D; const Mt: IFArr2D);
{ Mt = Transpose(ML)*DiagonalMatrix(DR) }
procedure MTxD(const ML: IFArr2D; const DR: IFArr1D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = MxDxMT = M*DiagonalMatrix(D)*Transpose(M) }
procedure MxDxMT(const M: IFArr2D; const D: IFArr1D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = MxDxMT = Transpose(M)*DiagonalMatrix(D)*M }
procedure MTxDxM(const M: IFArr2D; const D: IFArr1D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = M1xDxM2T = M*DiagonalMatrix(D)*Transpose(M) }
procedure M1xDxM2T(const M1, M2: IFArr2D; const D: IFArr1D; const Mt: IFArr2D);
{ matrix-matrix product: Mt = M1TxDxM2 = Transpose(M)*DiagonalMatrix(D)*M }
procedure M1TxDxM2(const M1, M2: IFArr2D; const D: IFArr1D; const Mt: IFArr2D);
{ Ascending 1-Darray sorting <br>
Based on: Sebastian Boßung, Sorting Algorithms in Pascal<br>
@link(http://privat.schlund.de/b/bossung/prog/delphi/sorting.html) }
procedure SortFArr1D(const A: IFArr1D; Ascending: boolean);
{ Ascending 1D-array sorting with sorting of the associated 2D-array <br>
by columns (idx = 2) or by rows (idx = 1) <br>
Based on: Sebastian Boßung, Sorting Algorithms in Pascal <br>
http://privat.schlund.de/b/bossung/prog/delphi/sorting.html }
procedure SortFArr2D(const A: IFArr1D; const B: IFArr2D; Ascending: boolean; st: TSliceType);
{ MT = Transpose(M) }
procedure MtT(const M, MT: IBArr2D); overload;
{ MT = Transpose(M) }
procedure MtT(const M, MT: IIArr2D); overload;
{ MT = Transpose(M) }
procedure MtT(const M, MT: IFArr2D); overload;
{ Complex Algebra - backward compatibility support }
{ MT = Transpose(M) }
procedure CMT(const MR, MI, MTR, MTI: IFArr2D);
{ MC = Conjugate(M); M = MC is OK }
procedure CMC(const MR, MI, MCR, MCI: IFArr2D);
{ MH = HTranspose(M) }
procedure CMH(const MR, MI, MHR, MHI: IFArr2D);
{ Dot (scalar) product of two vectors }
procedure CVxCV(const V1R, V1I, V2R, V2I: IFArr1D; var ReResult, ImResult: TFloat);
{ HConjugate(V1) Dot V2 }
procedure CVHxCV(const V1R, V1I, V2R, V2I: IFArr1D; var ReResult, ImResult: TFloat);
{ (V3R + i*V3I) = (M1R + i*M1I)*(V2R + i*V2I)
V3R = M1R*V2R - M1I*V2I; M3I = M1R*V2I + M1I*V2R }
procedure CMxCV(const M1R, M1I: IFArr2D; const V2R, V2I, V3R, V3I: IFArr1D);
{ (V3R + i*V3I) = (V1R + i*V1I)*(M2R + i*M2I)
V3R = V1R*M2R - V1I*M2I; V3I = V1R*M2I + V1I*M2R }
procedure CVxCM(const V1R, V1I: IFArr1D; const M2R, M2I: IFArr2D; const V3R, V3I: IFArr1D);
{ (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: IFArr2D);
{ (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: IFArr2D);
{ (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: IFArr2D);
{ (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: IFArr2D; const D: IFArr1D; const M3R, M3I: IFArr2D);
{ (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: IFArr2D; const D: IFArr1D; const M3R, M3I: IFArr2D);
{ ========================================================================== }
implementation
uses
Math, sysutils, classes, uDisplay;
end.
Back
Generated by Lore's Source to HTML Converter ( http://www.newty.de/lsc/index.html )