Linear algebra routines
Based on EBK&NVS Library for Turbo/Object Pascal
by Nikolai V. Shokhirev and Eugene B. Krissinel
author: Nikolai Shokhirev, http://www.shokhirev.com/nikolai.html
http://www.chem.arizona.edu/~shokhirn/nikolai.html)
author: Eugene B. Krissinel http://www.fh.huji.ac.il/~krissinel/
created: 2002.02.02
last modified: 2004.12.22
©Nikolai V. Shokhirev, 2002-2004 }
unit uDynLinAlg;
interface
uses
uMatTypes, uDynArrays;
{author Nikolai V. Shokhirev, 2004
| d1 e2 0 0 ... |
| 0 d2 e3 0 ... |
| ... |
| ... 0 dN-1 eN |
| ... 0 0 dN |
M is square dim x dim e.Dim1 = d.Dim1 = dim +}
procedure BiDiagMt(const d, e: IFArr1D; const M: IFArr2D);
{author Nikolai V. Shokhirev, 2004
| b1 c1 0 0 ... |
| a2 b2 c2 0 ... |
| ... |
| ... aN-1 bN-1 cN-1|
| ... 0 aN bN |
M is square dim x dim
a.Dim1 = b.Dim1 = c.Dim1 = dim +}
procedure TriDiagMt(const a, b, c: IFArr1D; const M: IFArr2D);
{ pascal version by Nikolai V. Shokhirev, 2004
Solution of linear system, a*x = b - translation from FORTRAN
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
*** Do not use if decomp has detected singularity. ***
input:
a[1..n,1..n] - triangularized matrix obtained from decomp.
b[1..n] - right hand side vector.
ipvt[1..n] - pivot vector obtained from decomp .
output:
b[1..n] - solution vector, x . +}
procedure FMMSolve(const a: IFArr2D; const b: IFArr1D; const ipvt: IIArr1D);
{ pascal version by Nikolai V. Shokhirev, 2004
Decomposes a matrix by Gaussian elimination and estimates the condition
of the matrix. - translation from FORTRAN
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
*** Use solve to compute solutions to linear systems. ***
input:
a[1..n,1..n] - matrix to be triangularized.
output:
a - contains an upper triangular matrix u and a permuted version of a lower
triangular matrix i-l so that (permutation matrix)*a = l*u .
cond - an estimate of the condition of a .
For the linear system a*x = b, changes in a and b may cause changes cond
times as large in x. If cond+1.0 = cond , a is singular to working
precision. cond is set to 1.0e+32 if exact singularity is detected.
ipvt - the pivot vector.
ipvt[k] = the index of the k-th pivot row
ipvt[n] = (-1)**(number of interchanges)
work space..
the vector work must be declared and included in the call. Its input
contents are ignored. Its output contents are usually unimportant.
The determinant of a can be obtained on output by
det(a) = ipvt[n] * a[1,1] * a[2,2] * ... * a[n,n]. +}
procedure FMMDecomp(const a: IFArr2D; var cond: TFloat; const ipvt: IIArr1D);
{ Author: Nikolai V. Shokhirev, 2004
Solution of linear system a*x = b
The procedure uses pascal translation of decomp and solve
from G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
input:
a[Lo1..Hi1,Lo2..Hi2] - n x n matrix
b[Lo1..Hi1] - right hand side n-vector
output:
b[Lo2..Hi2] - solution vector, x (n vector, the program changes its limits)
algorithm:
FNNDecomp(a, cond, ipvt);
FMMSolve(a, b, ipvt); +}
procedure DecompSolve(const a: IFArr2D; const b: IFArr1D; var cond: TFloat);
{ Pascal version by N. V. Shokhirev and E. B. Krissinel, 1991
Householder (tridiagonal) reduction of a real, symmetric matrix a.
a[Lo..Hi,Lo..Hi]
if (vectors = true) then
a is replaced by the orthogonal transformation matrix Q
d[Lo..Hi] returns the diaginal elements of the tridiagonal matrix
e[Lo..Hi] returns the off-giagonal elements with e[1] = 0.
Epsilon is elimination threshold }
procedure tred2(const a: IFArr2D; const d, e: IFArr1D;
vectors: boolean = true; Epsilon: TFloat = SqrtMinFloat);
{ Pascal version by Nikolai V. Shokhirev and Eugene B. Krissinel
QL algorithm with implicit shifts for a real tridiaginal symmetric matrix.
On input, d[Lo..Hi] is the vector of diagonal elements,
e[Lo..Hi] - subdiagonal elements, e[Lo] is arbitrary.
On input z is identity matrix for initially tridiagonal input matrix
or z is output from tred2 for general real symmetric matrix
On output, d contains eigenvalues
The k-th column z[i,k] returns the normalized eigenvectors,
corresponding to d[k]
IterMax is the maximal number of transformations
if (Iter <= Itermax) then err := 0 else err := Itermax +}
procedure tqli(const d, e: IFArr1D; const z: IFArr2D; var err: TInt;
vectors: boolean = true; IterMax: TInt = 999);
{ Author: Nikolai V. Shokhirev, 1994
Simple sort algorithm of the vector d[Lo..Hi] and corresponding columns
of matrix v[Lo1..Hi1,Lo..Hi]
if (ascending = true) then d[Lo] <= d[Lo+1] <=, ... <= d[Hi]
the columns are sorted only if (vectors = true) +}
procedure SortVal0(const d: IFArr1D; const v: IFArr2D;
vectors: boolean = true; ascending: boolean = true);
{ Author: Nikolai V. Shokhirev, 2002
Quick sort, Functionally the same as SortVal0 +}
procedure SortVal(const d: IFArr1D; const v: IFArr2D;
vectors: boolean = true; ascending: boolean = true);
{ Author: Nikolai V. Shokhirev, 2002
Diagonalization of the real symmetric matrix a[Lo..Hi,Lo..Hi] .
d[Lo..Hi] returns eigenvalues
if (vectors = true) then a returns eigenvectors by columns
the k-th eigenvector a[i,k] corresponds to the eigenvalue d[k]
if (sort = true) then values (and vectors) are sorted
if (ascending = true) then d[Lo] <= d[Lo+1] <=, . . . <= d[Hi]
Epsilon is tridiagonalization elimination threshold
IterMax is the maximal number of diagonalization transformations
err = 0 if diagonalization is successful,
err = IterMax if the number of iterations reached IterMax
Algorithm:
tred2(a, d, e, vectors, Epsilon);
tqli(d, e, a, err, vectors, IterMax);
SortVal(d, a, vectors, ascending); +}
procedure diag(const a: IFArr2D; const d: IFArr1D; sort: boolean = true;
vectors: boolean = true; ascending: boolean = true;
Epsilon: TFloat = SqrtMinFloat; IterMax: TInt = 999);
{ Author: Nikolai V. Shokhirev, 1994
Diagonalization of real symmetric matrix by JACOBI method
Input:
A[Lo..Hi,Lo..Hi] - the matrix to be diagonalized;
A[i,j] for i < j is destroyed
Output:
T[Lo..Hi,Lo..Hi] - the matrix of eigenvectors by columns: Tn[i] = T[i,n]
Eigen[Lo..Hi] - the vector of eigenvalues:
Eigen[Lo] <= Eigen[Lo+1] <= ... <= Eigen[Hi]
err = 0 O.K.
= IterMax if convergence has not beeb reached after IterMax iterations+}
procedure Jacobi(const A, T: IFArr2D; const Eigen: IFArr1D;
var err: integer; first: boolean = true);
{ author: Nikolai V. Shokhirev, 1994
Jacoby diagonalization of Hermitian matrix
Input:
A + i*B - the matrix to be diagonalized ([Lo..Hi,Lo..Hi])
Output:
U + i*V - the matrix of eigenvectors by columns ([Lo..Hi,Lo..Hi])
E[Lo..Hi] - the vector of eigenvalues: E[Lo] <= E[Lo+1] <= ... <= E[Hi]
err = 0 O.K.
= IterMax if convergence has not beeb reached after IterMax iterations
< 0 for incorrect limits
convergence criteria:
Norm2(OffD(A+iB))/2 < sqr(RelEps)*N*Norm2(A+iB)/2
OR NoRotation OR Iter >= IterMax +}
procedure CJacobi (const A, B, U, V: IFArr2D; const E: IFArr1D;
var err: TInt; first: boolean = true; order: boolean = true);
{ author: Nikolai V. Shokhirev, 1994
1-Based variant of Jacobi diagonalization of Hermitian matrix
Input:
N - dimension of the matrices
A + i*V - the matrix to be diagonalized
Output:
U + i*V - the matrix of eigenvectors by columns
E - the vector of eigenvalues: E[1] <= E[2] <= ... <= E[N]
err = 0 O.K.
= IterMax if convergence has not beeb reached after IterMax iterations
< 0 incorrect dimensions +}
procedure JacobiC (const A, B, U, V: IFArr2D; const E: IFArr1D;
var err: TInt; first: boolean = true; order: boolean = true);
{ pascal version by Nikolai V. Shokhirev 2002
This subroutine reduces a complex hermitian matrix to a real symmetric
tridiagonal matrix using unitary similarity transformations.
On input:
ar[1..n,1..n] and ai[1..n,1..n] contain the real and imaginary parts,
respectively, of the complex hermitian input matrix.
Only the lower triangle of the matrix need be supplied (a[i,j], i>=j).
On output:
ar[1..n,1..n] and ai[1..n,1..n] contain information about the unitary
transformations used in the reduction in their full lower triangles.
Their strict upper triangles and the diagonal of ar are unaltered.
d[1..n] contains the diagonal elements of the the tridiagonal matrix.
e[1..n] contains the subdiagonal elements of the tridiagonal matrix in
last n-1 its positions. e[1] is set to zero.
e2[1..n] contains the squares of the corresponding elements of e.
e2 may coincide with e if the squares are not needed.
tau[1..2,1..n] contains further information about the transformations.
calls pythag for dsqrt(a*a + b*b) .
(see the original EISPACK comment for htridi in implementation) +}
procedure htridi(const ar, ai, tau: IFArr2D; const d, e, e2: IFArr1D);
{ pascal variant by Nikolai V. Shokhirev 2002
This subroutine forms the eigenvectors of a complex Hermitian matrix
by back transforming those of the corresponding real symmetric tridiagonal
matrixdetermined by htridi.
On input:
ar[1..n,1..n] and ai[1..n,1..n] contain information about the unitary
transformations used in the reduction by htridi in their full lower
triangles except for the diagonal of ar.
tau[1..2,1..n] contains further information about the transformations.
m is the number of eigenvectors to be back transformed.
zr contains the eigenvectors to be back transformed in its first
m columns (from tqli or tql2).
On output:
zr[1..n,1..n] and zi[1..n,1..n] contain the real and imaginary parts,
respectively, of the transformed eigenvectors in their first m columns.
Note that the last component of each returned vector is real
and that vector euclidean norms are preserved.
( see the original EISPACK comment for htribk in implementation ). +}
procedure htribk(const ar,ai, tau: IFArr2D; m: TInt; const zr,zi: IFArr2D);
{ Author: Nikolai V. Shokhirev, 2002
Simple sort algorithm of the first m components of vector d[Lo..Hi]
and corresponding columns of matrices vr and vi ([Lo1..Hi1,Lo..Hi]).
if (ascending = true) then d[Lo] <= d[Lo+1] <=, . . . <= d[Hi]
the columns are sorted only if (vectors = true) }
procedure SortValC(const d: IFArr1D; const vr, vi: IFArr2D; m: TInt;
vectors: boolean = true; ascending: boolean = true);
{ author: Nikolai V. Shokhirev 2002
Diagonalization of Hermitian matrix a = ar + i*ai.
d returns eigenvalues
if (m > 0) then z = zr + i*zi returns m eigenvectors by columns
if (sort = true) then values (and vectors) are sorted
if (ascending = true) then d[1] <= d[2] <=, . . . <= d[n]
err = 0 if diagonalization is successful,
err = IterMax if the number of iterations reached IterMax
err < 0 for incorrect limits (see implementation)
Algorithm:
htridi(ar, ai, tau, d, e, e2);
tqli(d, e, zr, err);
htribk(ar, ai, tau, m, zr,zi);
SortValC(d, a, m, true, ascending);
Lim checks:
if (ar.Dim1<>ar.Dim2) then exit;
n := ar.Dim1;
if zr.Dim1 < n then exit;
if not SameLim(ar,ai) or not SameLim(zr,zi) then exit;
if (d.Dim1 < n) then exit;
mm := min(ar.Dim2,zr.Dim2);
if m > 0 then
if m > mm then m := mm; +}
procedure diagc(const ar, ai, zr, zi: IFArr2D; const d: IFArr1D; m: TInt;
var err: TInt; sort: boolean = true; ascending: boolean = true);
{authors: Eugene B. Krissinel and Nikolai V. Shokhirev, 1991
Fast Inversion of the matrix A by the Gauss-Jordan elimination algorithm
Input:
A[Lo..Hi,Lo..Hi] - the matrix to be inversed.
Output:
A - the inversed matrix
Det - the determinant of A +}
procedure FastInverse(const A: IFArr2D; var Det: TFloat);
{ pascal version by Eugene B. Krissinel, 1991
Singular Value Decomposition - translation from FORTRAN
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
This subroutine determines the singular value decomposition
A = U * W * VT of a real M by N rectangular matrix.
Householder bidiagonalization and a variant of the QR algorithm are used.
SVD of the matrix A: A = U * W * VT
N1xN2 N1xN2 N2xN2 N2xN2
Input:
Case N1 > N2
NA = M = N1, N = N2
Dimensions: A[M,N],W[N],U[M,N],V[M,N],RV1[N]
Case N1 < N2
NA = N = N1, M = N2
Dimensions: A[N,M],W[N],U[M,N],V[M,N],RV1[N]
*** Always M >= N
A contains the rectangular input matrix to be decomposed.
MatU should be set to true if the U matrix in the
decomposition is desired, and to false otherwise.
MatV should be set to true if the V matrix in the
decomposition is desired, and to false otherwise.
Output:
A is unaltered (unless overwritten by U or V).
W contains the N (non-negative) singular values of A (the diagonal
elements of s). They are unordered. If an error exit is made,
the singular values should be correct for indices ierr+1,ierr+2,...,N.
U contains N orthogonal column U-vectors of the decomposition,
if MatU has been set to true. Otherwise U is used as a temporary array.
U may coincide with A.
If an error exit is made, the columns of U corresponding
to indices of correct singular values should be correct.
*** In the case N1 < N2 the last M-N rows contain zero
V contains N orthogonal column V-vectors of the decomposition,
if MatV has been set to true. Otherwise V is not referenced.
V may also coincide with A if U is not needed.
If an error exit is made, the columns of V corresponding
to indices of correct singular values should be correct.
*** In the case N1 > N2 the last M-N rows are zero
RetCode
= 0 for normal exit,
= k if the k-th singular value has not been
determined after 30 iterations.
RV1[N] is a temporary storage array.
This is a modified version of a routine from the eispack collection
by the NATS project modified to eliminate machep +}
procedure FMMSVD(NA,M,N: integer; const A,U,V: IFArr2D; const W,RV1: IFArr1D;
MatU,MatV: boolean; var RetCode: integer);
{ escending SVD sorting by Eugene B. Krissinel, 1991
W[1..N], U[1..M,1,..N], V[1..M,1,..N]
}
procedure OrderSVD(M, N: integer; const U, V: IFArr2D; const W: IFArr1D;
MatU, MatV: boolean);
{ Author: Nikolai V. Shokhirev, 2002
Descending singular value sorting with sorting of U and V
Based on: Sebastian Boßung, Sorting Algorithms in Pascal
http://privat.schlund.de/b/bossung/prog/delphi/sorting.html }
procedure SortSVD(const W: IFArr1D; const U, V: IFArr2D; MatU, MatV: boolean);
{Author: Nikolai V. Shokhirev, 2003
used in SVD, based om FMMSVD
*** All arrays are 1-based ***
M := max(A.Dim1,A.Dim2);
N := min(A.Dim1,A.Dim2);
if FullMat then mn := M else mn := N;
if (U.Dim1 < M) or (U.Dim2 < mn) or ( V.Dim1 < N) or
(V.Dim2 < N) or (W.Dim1 < mn) then RetCode := -1; +}
procedure SVD_(const A, U, V: IFArr2D; const W: IFArr1D; var RetCode: integer;
MatU, MatV: boolean; FullMat: boolean = false; ItnLimit: integer = 30);
{ Author: Nikolai V. Shokhirev, 2002.
A is a rectangular D1xD2 matrix: A[L1..H1,L2..H2]
if FullMat = true then
U is a square D1xD1 matrix: U[L1..H1,1..D1],
V is a square D2xD2 matrix: V[L2..H2,1..D2] and
UT*A*V = W
where W is a rectangular D1xD2 matrix: W[1..D1,1..D2]
and only W[i,i] <> 0 i = 1.. N, N = min(D1,D2)
if FullMat = false then
U is a D1xN matrix: U[L1..H1,1..N], N = min(D1,D2)
V is a D2xN matrix: V[L2..H2,1..N] and
UT*A*V = W
where W is a square NxN diagonal matrix: W[1..N,1..N]
the initial matrix A can be expressed as the following
product: A = U*W*VT (Singular Value Decomposition)
N1 > N2
Full: A[N1,N2] = U[N1,N1] * W[N1,N2] * VT[N2,N2] V[N2,N2] mn = N1
Small: A[N1,N2] = U[N1,N2] * W[N2,N2] * VT[N2,N2] V[N2,N2] mn = N2
N1 < N2
Full: A[N1,N2] = U[N1,N1] * W[N1,N2] * VT[N2,N2] V[N2,N2] mn = N2
Small: A[N1,N2] = U[N1,N1] * W[N1,N1] * VT[N1,N2] V[N2,N1] mn = N1 +}
procedure SVD(const A, U, V: IFArr2D; const W: IFArr1D; var RetCode: integer;
MatU: boolean = true; MatV: boolean = true;
FullMat: boolean = false; ItnLimit: integer = 30);
{ Author: Nikolai V. Shokhirev, 2004.
|d1 e2 0 ... | | x1 | | b1 |
|0 d2 e3 ... | | x2 | | b2 |
| ... |*| .. |=| .. |
| ... dN-1 eN | |xN-1| |bN-1|
| ... 0 dN | | xN | | bN |
a[1..n], b[1..n], c[1..n], d[1..n], x[1..n]
*** can be used in combination with householder *** +}
procedure BiDiagSLE(d, e, b, x: IFArr1D);
{ Author: Nikolai V. Shokhirev, 2002.
|b1 c1 0 0 ... | | x1 | | d1 |
|a2 b2 c2 0 ... | | x2 | | d2 |
| 0 a3 b3 c3 ... | | x3 | | d3 |
| ... |*| .. |=| .. |
| ... aN-1 bN-1 cN-1| |xN-1| |dN-1|
| ... 0 aN bN | | xN | | dN |
a[1..n], b[1..n], c[1..n], d[1..n], x[1..n] +}
procedure TriDiagSLE(const a, b, c, d, x: IFArr1D);
{ author: Eugene B. Krissinel, 1991
Perturbated Cholesky Decomposition -}
procedure CholDecomp(N: integer; const HDiag: IFArr1D; MaxOff: TFloat;
const L: IFArr2D; var MaxAdd: TFloat);
{ author: Eugene B. Krissinel, 1991
Cholesky Decomposition A = L*LT
Input
L = A - positive-definite symmetric matrix
the upper triangle of L[i,j], j > i is not used
Output
lower triangle L[i,j], i <= j contains L
upper triangle L[i,j], j > i is set to 0 +}
procedure Cholesky(const L: IFArr2D);
{ author: Eugene B. Krissinel, 1991
Solution of L*Y = B;
L is a lower-triangle matrix: L[i,j], i >= j }
procedure LSolve(const L: IFArr2D; const B,Y: IFArr1D);
{ author: Eugene B. Krissinel, 1991
Solution of LT*X = Y L[j,i], j >= i
LT = Transpose(L) is a upper-triangle matrix
L is a lower-triangle matrix: L[i,j], i >= j }
procedure LTSolve(const L: IFArr2D; const Y,X: IFArr1D);
{ author: Eugene B. Krissinel, 1991
Solution of L*LT*X = B by the Cholesky method
L*Y = B // LSolve
LT*X = Y // LTSolve +}
procedure ChSolve (const L: IFArr2D; const B, X: IFArr1D);
{author: Eugene B. Krissinel, 1991
Cholesky decomposition of the band matrix H[1..p+1,1..N].
The decomposition will be written into L[1..p+1,1..N]. -}
procedure CholBandDec (p: integer; const H,L: IFArr2D; var err: TInt);
{ author: Eugene B. Krissinel, 1991
Cholesky L - Solution of L*Y = B (for given B)
for the band matrix L[1..p+1,1..N]. -}
procedure LBandSolve(p: integer; const L: IFArr2D; const B,Y: IFArr1D;
var err: TInt);
{ author: Eugene B. Krissinel, 1991
Cholesky LT - Solution of LT*X = Y (for given Y)
for the band matrix L[1..p+1,1..N]. -}
procedure LTBandSolve(p: integer; const L: IFArr2D; const Y,X: IFArr1D;
var err: TInt);
{ author: Eugene B. Krissinel, 1991
Solution of the equation L*LT*S = G by the Cholesky method
for the band matrix L[1..p+1,1..N]. -}
procedure ChBandSolve(p: integer; const L: IFArr2D; const G,S: IFArr1D;
var err: TInt);
{ authors: Eugene B. Krissinel and Nikolai V. Shokhirev, 1991
QR-decomposition of the n x n matrix M: M = Q*R
Q - orthogonal m x n matrix so that QT*Q = 1 (n x n)
R - upper triangle matrix: R[i,j] = 0 , i > j
Input
M[Lo..Hi,Lo..Hi] - matrix to be decomposed
N[Lo..Hi], D[Lo..Hi]
Oitput
D[k] - the first non-zero component of the k-th transformation vector uk
M[j,k], j =k=1,..,Hi - non-zero components of the k-th transform vector
M[j,k], j <= k - non-zero part of the upper-right R matrix
N[k] = <uk|uk>/2
Algorithm
The program performs the Householder's transformations:
M(k+1) = Q(k)*M(k), M(1) = M
Q =Q(n-1)* .. *Q(2)*Q(1)
Q(k) = 1 - 2*|uk><uk|/<uk|uk> = 1 - |uk><uk|/(<uk|uk>/2)
define v as
v[j] = 0, j = 1,..,k-1
v[j] = M(k)[j,k], j = k,..,n
uk[j] = v[j], j <> k
uk[k] = v[k]+sign(v[k])*|v|, |v| = sqrt(<v|v>)
Q(k) does not change the columns 1,..,k-1;
makes M[j,k] = 0, j = k+1,..,m;
changes the columns k+1,..,n +}
procedure QRDecomp(const M: IFArr2D; const N, D: IFArr1D; tol: TFloat = MinFloat);
{authors: Eugene B. Krissinel and Nikolai V. Shokhirev, 1991
solution of the equation Q*R*X = B
Input
M, D and N are the output arrays from QRDecomp
B - r.h.s. vector
Otput
the computed X is returned in B
M, D and N are not changed +}
procedure QRSolve(const M: IFArr2D; const N,D, B: IFArr1D);
{ author: Nikolai V. Shokhirev, 2002
solution of a set of linear equations by the Gauss Jordan elimination
with partial pivoting
a[Lo1..Hi1,Lo2..Hi2] = |m|b| - n x M matrix
m - N x N coefficient matrix, b - N x K r.h.s. matrix
m is destroyed, b returns K solutions
if the r.h.s = Identity matrix then b returns the inverse matrix of M
DoBalance = true - balancing so that MaxAbs(row) = 1
and check for degeneracy
- usually not necessary to use this option +}
procedure GaussJordan(const a: IFArr2D; DoBalance: boolean; var err: TInt);
{ author: Nikolai V. Shokhirev, 2002
Direct use of GaussJordan for solving m*x =b
m[Lo1..Hi1,Lo2..Hi2], x[Lo2..Hi2], b[Lo1..Hi1], x[Lo2..Hi2],
x = b is OK, in this case the limits of x will be reset to [Lo2..Hi2] +}
procedure GJSolve(const m: IFArr2D; DoBalance: boolean; const x, b: IFArr1D; var err: TInt);
{ authors Eugene B. Krissinel and Nikolai V. Shokhirev, 1991
solution of a set of linear equations a*x = b
by the Gauss Jordan elimination with partial pivoting
a[Lo1..Hi1,Lo2..Hi2] is a M x N matrix of the
system of M linear equations of N variables
if M > N then the "best" N equations are uaed
b[Lo1..Hi1] is a M-vector of the r.h.s.
x is a N-vector of the solution, Lims are set to [Lo2..Hi2]
b and x can be the same
a and b are destroyed +}
procedure FastSolve(const a: IFArr2D; const b, x: IFArr1D);
{ Nikolai V. Shokhirev, 1994
solution of a set of linear equations by the Gauss Jordan elimination
with partial pivoting
a is a M x N matrix of the system of M linear equations of N variables
Use is a boolean N-vector: if Use[i] = false then the i-th variable is
excluded from a solution and added to the r.h.s.
if M > Nvar then the "best" Nvar equations are uaed;
Here Nvar is the number of selected variables (Use[i] = true)
b is a M-vector of the r.h.s., x is a N-vector of the solution
a is NOT destroyed, but b IS destroyed +}
procedure MaskSolve(const m: IFArr2D; const x, b: IFArr1D; const Use: IBArr1D);
{ Nikolai V. Shokhirev, 2001
Gram-Schmidt Orthonormalization of vectors
input:
U - matrix of column-vectors:
U[Lo1..Hi1,Lo2..Hi2] = |Col[Lo2],Col[Lo2+1],...,Col[Hi2]|
m - number of required vectors
output:
U - the first m columns contain orthogonal and normalized vectors
m - actual number of independent vectors +}
procedure GramSchmidt(const U: IFArr2D; var m: TInt);
{ Nikolai V. Shokhirev, 2003
Householder reduction (m >= n) extracted from xsvd
A = U * J * VT UT * A * V = J
m x n m x n n x n n x n n x m m x n n x n n x n
Full (square) U:
A = U * J * VT UT * A * V = J
m x n m x m m x n n x n m x m m x n n x n m x n
|q1 e2 |
| q2 e2 |
| q2 e2 . |
J = | . . en|
| . qn|
|-----------------|
| 0 | +}
procedure householder(const a: IFArr2D; const q, e: IFArr1D; const u, v: IFArr2D;
eps, tol: TFloat; fullmatrix: boolean = false);
{Nikolai V. Shokhirev, 2003
Singular Value Decomposition (m >= n)
J.H.Wilkinson, C.Reinsch, Handbook for Automatic Computation
A = U * Q * VT UT * A * V = Q
m x n m x n n x n n x n n x m m x n n x n n x n
Full (square) U:
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 +}
procedure xsvd(eps, tol: TFloat; const a: IFArr2D; const q: IFArr1D;
const u, v: IFArr2D; withu: boolean = true;
withv: boolean = true; fullmat: boolean = false);
implementation
uses
Math, SysUtils, uDynArrUtils;
end.
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)