Up \ObjAlg
{@abstract(Linear algebra routines <br>
Based on EBK&NVS Library for Turbo/Object Pascal <br>
by Nikolai V. Shokhirev and Eugene B. Krissinel )
author: Nikolai Shokhirev)
author: Eugene B. Krissinel http://www.fh.huji.ac.il/~krissinel/)
created: February 02, 2002)
last modified: June 06, 2003)
©Nikolai V. Shokhirev, 2002-2003
}
unit uLinAlg;
interface
uses
uMatTypes, uTensors, uArrays;
{<pre> Pascal version by Nikolai V. Shokhirev and Eugene B. Krissinel <br>
Householder (tridiagonal) reduction of a real, symmetric matrix a.<br>
if (vectors = true) then<br>
a is replaced by the orthogonal transformation matrix Q<br>
d returns the diaginal elements of the tridiagonal matrix<br>
e returns the off-giagonal elements with e[1] = 0.<br>
Epsilon is elimination threshold </pre> }
procedure tred2(const a: ITensor2D; const d, e: ITensor1D;
vectors: boolean = true; Epsilon: TFloat = SqrtMinFloat);
{<pre> Pascal version by Nikolai V. Shokhirev and Eugene B. Krissinel <br>
QL algorithm with implicit shifts<br>
for a real tridiaginal symmetric matrix.<br>
On input, d is the vector of diagonal elements,<br>
e - subdiagonal elements, e[1] is arbitrary.<br>
On input z is identity matrix for initially tridiagonal input matrix<br>
or z is output from tred2 for general real symmetric matrix <br>
On output, d contains eigenvalues <br>
The k-th column z[i,k] returns the normalized eigenvectors,<br>
corresponding to d[k]<br>
IterMax is the maximal number of transformations<br>
if (Iter <= Itermax) then err := 0 else err := Itermax</pre> }
procedure tqli(const d, e: ITensor1D; const z: ITensor2D; var err: TInt;
vectors: boolean = true; IterMax: TInt = 999);
{<pre> Author: Nikolai V. Shokhirev, 2002 <br>
Simple sort algorithm of the vector d and corresponding columns of matrix v.<br>
if (ascending = true) then d[1] <= d[2] <=, . . . <= d[n]<br>
the columns are sorted only if (vectors = true)<br>
if vectors = false then v can be set to nil </pre> }
procedure SortVal(const d: ITensor1D; const v: ITensor2D;
vectors: boolean = true; ascending: boolean = true);
{<pre> Author: Nikolai V. Shokhirev, 2002 <br>
Diagonalization of the real symmetric n by n matrix a.<br>
d returns eigenvalues<br>
if (vectors = true) then a returns eigenvectors by columns<br>
the k-th eigenvector a[i,k] corresponds to the eigenvalue d[k]<br>
if (sort = true) then values (and vectors) are sorted<br>
if (ascending = true) then d[1] <= d[2] <=, . . . <= d[n]<br>
Epsilon is tridiagonalization elimination threshold<br>
IterMax is the maximal number of diagonalization transformations<br>
err = 0 if diagonalization is successful,<br>
err = IterMax if the number of iterations reached IterMax<br>
err < 0 for incorrect limits (see implementation)<br>
Algorithm:<br>
tred2(a, d, e, vectors, Epsilon);<br>
tqli(d, e, a, err, vectors, IterMax);<br>
Order(d, a, vectors, ascending); </pre> }
procedure diag(const a: ITensor2D; const d: ITensor1D; var err: TInt;
sort: boolean = true; vectors: boolean = true; ascending: boolean = true;
Epsilon: TFloat = SqrtMinFloat; IterMax: TInt = 999);
{<pre> Author: Nikolai V. Shokhirev, 1994 <br>
Diagonalization of real symmetric matrix by JACOBI method<br>
Input:<br>
A - the matrix to be diagonalized; A[i,j] for i < j is destroyed<br>
Output:<br>
T - the matrix of eigenvectors by columns: Tn[i] = T[i,n]<br>
Eigen - the vector of eigenvalues: Eigen[1] < Eigen[2] < ... < Eigen[N]<br>
err = 0 O.K.<br>
= IterMax if convergence has not beeb reached after IterMax iterations</pre> }
procedure Jacobi(const A, T: ITensor2D; const Eigen: ITensor1D;
var err: integer; first: boolean = true);
{ <pre>author: Nikolai V. Shokhirev 1994 <br>
Jacoby diagonalization of Hermitian matrix<br>
Input:<br>
A + i*B - the matrix to be diagonalized<br>
Output:<br>
U + i*V - the matrix of eigenvectors by columns<br>
E - the vector of eigenvalues: E[1] < E[2] < ... < E[N]<br>
err = 0 O.K.<br>
= IterMax if convergence has not beeb reached after IterMax iterations<br>
< 0 for incorrect limits </pre> }
procedure CJacobi (const A, B, U, V: ITensor2D; const E: ITensor1D;
var err: TInt; first: boolean = true; order: boolean = true);
{ <pre> author: Nikolai V. Shokhirev 1994<br>
1-Based variant of Jacobi diagonalization of Hermitian matrix<br>
Input:<br>
N - dimension of the matrices<br>
A + i*V - the matrix to be diagonalized<br>
Output:<br>
U + i*V - the matrix of eigenvectors by columns<br>
E - the vector of eigenvalues: E[1] < E[2] < ... < E[N]<br>
err = 0 O.K.<br>
= IterMax if convergence has not beeb reached after IterMax iterations<br>
< 0 incorrect dimensions </pre> }
procedure JacobiC (const A, B, U, V: ITensor2D; const E: ITensor1D;
var err: TInt; first: boolean = true; order: boolean = true);
{ <pre>pascal version by Nikolai V. Shokhirev 2002<br>
This subroutine reduces a complex hermitian matrix to a real symmetric <br>
tridiagonal matrix using unitary similarity transformations.<br>
On input: <br>
ar and ai contain the real and imaginary parts, respectively, of the complex<br>
hermitian input matrix.<br>
Only the lower triangle of the matrix need be supplied.<br>
On output: <br>
ar and ai contain information about the unitary transformations used in the <br>
reduction in their full lower triangles.<br>
Their strict upper triangles and the diagonal of ar are unaltered.<br>
d contains the diagonal elements of the the tridiagonal matrix.<br>
e contains the subdiagonal elements of the tridiagonal matrix in its last n-1 <br>
positions. e(1) is set to zero.<br>
e2 contains the squares of the corresponding elements of e. e2 may coincide<br>
with e if the squares are not needed. <br>
tau contains further information about the transformations.<br>
<br>
calls pythag for dsqrt(a*a + b*b) .<br>
(see the original EISPACK comment for htridi in implementation) </pre> }
procedure htridi(const ar, ai, tau: ITensor2D; const d, e, e2: ITensor1D);
{ <pre>pascal variant by Nikolai V. Shokhirev 2002 <br>
This subroutine forms the eigenvectors of a complex Hermitian matrix by back <br>
transforming those of the corresponding real symmetric tridiagonal matrix<br>
determined by htridi. <br>
On input: <br>
ar and ai contain information about the unitary transformations used in the<br>
reduction by htridi in their full lower triangles except for the<br>
diagonal of ar.<br>
tau contains further information about the transformations.<br>
m is the number of eigenvectors to be back transformed.<br>
zr contains the eigenvectors to be back transformed in its first m columns <br>
(from tqli or tql2).<br>
On output: <br>
zr and zi contain the real and imaginary parts, respectively, of the<br>
transformed eigenvectors in their first m columns.<br>
<br>
Note that the last component of each returned vector is real <br>
and that vector euclidean norms are preserved.<br>
( see the original EISPACK comment for htribk in implementation ). </pre> }
procedure htribk(const ar,ai, tau: ITensor2D; m: TInt; const zr,zi: ITensor2D);
{<pre> Author: Nikolai V. Shokhirev, 2002 <br>
Simple sort algorithm of the vector d <br>
and corresponding columns of matrices vr and vi.<br>
if (ascending = true) then d[1] <= d[2] <=, . . . <= d[n]<br>
the columns are sorted only if (vectors = true)<br>
if vectors = false then vr and vi can be set to nil </pre> }
procedure SortValC(const d: ITensor1D; const vr, vi: ITensor2D; m: TInt;
vectors: boolean = true; ascending: boolean = true);
{<pre> author: Nikolai V. Shokhirev 2002 <br>
Diagonalization of Hermitian matrix a = ar + i*ai.<br>
d returns eigenvalues<br>
if (m > 0) then z = zr + i*zi returns m eigenvectors by columns<br>
if (sort = true) then values (and vectors) are sorted<br>
if (ascending = true) then d[1] <= d[2] <=, . . . <= d[n]<br>
err = 0 if diagonalization is successful,<br>
err = IterMax if the number of iterations reached IterMax<br>
err < 0 for incorrect limits (see implementation) <br>
Algorithm:<br>
htridi(ar, ai, tau, d, e, e2);<br>
tqli(d, e, zr, err); <br>
htribk(ar, ai, tau, m, zr,zi);<br>
SortValC(d, a, m, true, ascending); </pre> }
procedure diagc(const ar, ai, zr, zi: ITensor2D; const d: ITensor1D; m: TInt;
var err: TInt; sort: boolean = true; ascending: boolean = true);
{<pre> Fast Inversion of the matrix A by the Gauss-Jordan elimination algorithm <br>
Input:<br>
A - the matrix to be inversed.<br>
Output:<br>
A - the inversed matrix<br>
Det - the determinant of A<br>
err - the error key :<br>
= 0 - OK <br>
= -1 - matrix is degenerate </pre> }
procedure FastInverse(const A: IMatrix; var Det: TFloat; var err, idx: integer);
{<pre> pascal version by Eugene B. Krissinel, 1991<br>
Singular Value Decomposition - translation from FORTRAN <br>
G.E.Forsythe, M.A.Malcolm, C.B.Moler.<br>
Computer methods for mathematical computations.<br>
Prentice-Hall, 1977.<br>
This subroutine determines the singular value decomposition<br>
A = U * W * VT of a real M by N rectangular matrix.<br>
Householder bidiagonalization and a variant of the QR algorithm are used.<br>
SVD of the matrix A: A = U * W * VT<br>
N1xN2 N1xN2 N2xN2 N2xN2<br>
Input:<br>
Case N1 > N2 <br>
NA = M = N1, N = N2 <br>
Dimensions: A[M,N],W[N],U[M,N],V[M,N],RV1[N] <br>
Case N1 < N2 <br>
NA = N = N1, M = N2<br>
Dimensions: A[N,M],W[N],U[M,N],V[M,N],RV1[N]<br>
*** Always M >= N <br>
A contains the rectangular input matrix to be decomposed.<br>
MatU should be set to true if the U matrix in the <br>
decomposition is desired, and to false otherwise.<br>
MatV should be set to true if the V matrix in the <br>
decomposition is desired, and to false otherwise.<br>
Output: <br>
A is unaltered (unless overwritten by U or V).<br>
W contains the N (non-negative) singular values of A (the diagonal <br>
elements of s). They are unordered. If an error exit is made,<br>
the singular values should be correct for indices ierr+1,ierr+2,...,N.<br>
U contains N orthogonal column U-vectors of the decomposition, <br>
if MatU has been set to true. Otherwise U is used as a temporary array.<br>
U may coincide with A.<br>
If an error exit is made, the columns of U corresponding <br>
to indices of correct singular values should be correct.<br>
*** In the case N1 < N2 the last M-N rows contain zero<br>
V contains N orthogonal column V-vectors of the decomposition,<br>
if MatV has been set to true. Otherwise V is not referenced.<br>
V may also coincide with A if U is not needed.<br>
If an error exit is made, the columns of V corresponding <br>
to indices of correct singular values should be correct. <br>
*** In the case N1 > N2 the last M-N rows are zero <br>
RetCode is set to <br>
0 for normal exit,<br>
k if the k-th singular value has not been<br>
determined after 30 iterations.<br>
RV1[N] is a temporary storage array.<br>
This is a modified version of a routine from the eispack collection by the <br>
NATS project modified to eliminate machep </pre> }
procedure FMMSVD(NA,M,N: integer; const A,U,V: IMatrix; const W,RV1: IVector;
MatU,MatV: boolean; var RetCode: integer);
{<pre> SVD sorting by Eugene B. Krissinel, 1991 </pre> }
procedure OrderSVD(M, N: integer; const U, V: IMatrix; const W: IVector;
MatU,MatV: boolean);
{<pre> Author: Nikolai V. Shokhirev, 2002 <br>
Descending singular value sorting with sorting of U and V<br>
Based on: Sebastian Boßung, Sorting Algorithms in Pascal <br>
http://privat.schlund.de/b/bossung/prog/delphi/sorting.html </pre> }
procedure SortSVD(const W: ITensor1D; const U, V: ITensor2D; MatU, MatV: boolean);
{<pre> used in SVD, based om FMMSVD <pre>}
procedure SVD_(const A, U, V: IMatrix; const W: IVector; var RetCode: integer;
MatU, MatV: boolean; FullMat: boolean = false; ItnLimit: integer = 30);
{<pre> Author: Nikolai V. Shokhirev, 2002. <br>
A is a rectangular D1xD2 matrix: A[L1..H1,L2..H2] <br>
if FullMat = true then <br>
U is a square D1xD1 matrix: U[L1..H1,1..D1], <br>
V is a square D2xD2 matrix: V[L2..H2,1..D2] and <br>
UT*A*V = W <br>
where W is a rectangular D1xD2 matrix: W[1..D1,1..D2]<br>
and only W[i,i] <> 0 i = 1.. N, N = min(D1,D2) <br>
<br>
if FullMat = false then<br>
U is a D1xN matrix: U[L1..H1,1..N], N = min(D1,D2)<br>
V is a D2xN matrix: V[L2..H2,1..N] and<br>
UT*A*V = W<br>
where W is a square NxN diagonal matrix: W[1..N,1..N]<br>
the initial matrix A can be expressed as the following<br>
product: A = U*W*VT (Singular Value Decomposition) </pre> }
procedure SVD(const A, U, V: IMatrix; const W: IVector; var RetCode: integer;
MatU: boolean = true; MatV: boolean = true;
FullMat: boolean = false; ItnLimit: integer = 30);
{<pre>
|b1 c1 0 0 ... | | x1 | | d1 | <br>
|a2 b2 c2 0 ... | | x2 | | d2 | <br>
| 0 a3 b3 c3 ... | | x2 | | d2 | <br>
| ... |*| .. |=| .. | <br>
| ... aN-1 bN-1 cN-1| |xN-1| |dN-1| <br>
| ... 0 aN bN | | xN | | dN | <br> </pre>}
procedure TriDiagSLE(a, b, c, d, x: IVector; var err: TInt);
{<pre> Perturbated Cholesky Decomposition </pre> }
procedure CholDecomp(N: integer; const HDiag: IVector; MaxOff, MachEps: TFloat;
const L: IMatrix; var MaxAdd: TFloat);
{<pre> Cholesky Decomposition A = L*LT </pre> }
procedure Cholesky( const L: IMatrix; var err: TInt);
{<pre> Cholesky L - Solution of L*Y = B (for given B); L[i,j], i >= j </pre> }
procedure LSolve(const L: IMatrix; const B,Y: IVector; var err: TInt);
{<pre> Cholesky LT - Solution of LT*X = Y (for given Y); L[j,i], j >= i </pre> }
procedure LTSolve(const L: IMatrix; const Y,X: IVector; var err: TInt);
{<pre> Solution of the equation L*LT*S = G by the Cholesky method </pre> }
procedure ChSolve(const L: IMatrix; const G,S: IVector; var err: TInt);
{<pre> Cholesky decomposition of the band matrix H[1..p+1,1..N].<br>
The decomposition will be written into L[1..p+1,1..N]. </pre> }
procedure CholBandDec (p: integer; const H,L: IMatrix; var err: TInt);
{<pre> Cholesky L - Solution of L*Y = B (for given B) <br>
for the band matrix L[1..p+1,1..N]. </pre> }
procedure LBandSolve(p: integer; const L: IMatrix; const B,Y: IVector; var err: TInt);
{<pre> Cholesky LT - Solution of LT*X = Y (for given Y) <br>
for the band matrix L[1..p+1,1..N]. </pre> }
procedure LTBandSolve(p: integer; const L: IMatrix; const Y,X: IVector; var err: TInt);
{<pre> Solution of the equation L*LT*S = G by the Cholesky method <br>
for the band matrix L[1..p+1,1..N]. </pre> }
procedure ChBandSolve(p: integer; const L: IMatrix; const G,S: IVector; var err: TInt);
{<pre> QR Decomposition of M </pre> }
procedure QRDecomp(const M: IMatrix; const M1,M2: IVector; var err: TInt);
{<pre> This routine solves R*X = B for X , <br>
where R is placed in the M by the QR-Decomp routine ; <br>
X is returned in B .</pre> }
procedure RSolve(const M: IMatrix; const M2,B: IVector; var err: TInt);
{<pre> This routine solves the equation (QR)*X = B,<br>
where Q and R are placed in the M by the QR-Decomp routine; <br>
the computed X is returned in B. </pre> }
procedure QRSolve(const M: IMatrix; const M1,M2,B: IVector; var err: TInt);
{<pre> solution of a set of linear equations by the Gauss Jordan elimination <br>
with partial pivoting<br>
a = |m|b| - n x M matrix<br>
m - N x N coefficient matrix, b - N x K r.h.s. matrix<br>
m is destroyed, b returns K solutions<br>
if the r.h.s = Identity matrix then b returns the inverse matrix of M<br>
DoBalance = true - balancing so that MaxAbs(row) = 1 and check for degeneracy<br>
- usually not necessary to use this option<br>
err = -1 if a.Dim1 <= a.Dim2 // nothing to solve<br>
err = 0 - OK<br>
err > 0 - err = rank+1 </pre> }
procedure GaussJordan(const a: IMatrix; DoBalance: boolean; var err: TInt);
{<pre> Direct use of GaussJordan for solving m*x =b </pre> }
procedure GJSolve(const m: IMatrix; DoBalance: boolean; const x, b: IVector; var err: TInt);
{<pre> solution of a set of linear equations by the Gauss Jordan elimination <br>
with partial pivoting<br>
a is a M x N matrix of the system of M linear equations of N variables<br>
if M > N then the "best" N equations are uaed<br>
b is a M-vector of the r.h.s., x is a N-vector of the solution<br>
b and x can be the same; a and b are destroyed<br>
err = -2 - a.Lo1 <> b.Lo<br>
err = -1 - a.Dim1 < a.Dim2 (too fiew equations)<br>
err = 0 - OK<br>
err > 0 - err = rank + 1 </pre> }
procedure FastSolve(const a: IMatrix; const b, x: IVector; var err: TInt);
{<pre> solution of a set of linear equations by the Gauss Jordan elimination<br>
with partial pivoting<br>
a is a M x N matrix of the system of M linear equations of N variables<br>
Use is a boolean N-vector: if Use[i] = false then the i-th variable is<br>
excluded from a solution and added to the r.h.s.<br>
if M > Nvar then the "best" Nvar equations are uaed;<br>
Here Nvar is the number of selected variables (Use[i] = true)<br>
b is a M-vector of the r.h.s., x is a N-vector of the solution<br>
a is NOT destroyed, but b IS destroyed<br>
err = -3 - (m.Lo2 <> Use.Lo) or (m.Hi2 > Use.Hi)<br>
err = -2 - a.Lo1 <> b.Lo<br>
err = -1 - a.Dim1 < a.Dim2 (too fiew equations)<br>
err = 0 - OK<br>
err > 0 - err = rank + 1 </pre> }
procedure MaskSolve(const m: IMatrix; const x, b: IVector; const Use: IBVector;
var err: TInt);
{<pre> Gram-Schmidt Orthonormalization of vectors <br>
input: <br>
U - matrix of column-vectors: <br>
U[Lo1..Hi1,Lo2..Hi2] = |Col[Lo2],Col[Lo2+1],...,Col[Hi2]| <br>
m - number of required vectors <br>
output: <br>
U - the first m columns contain orthogonal and normalized vectors <br>
m - actual number of vectors <pre>}
procedure GramSchmidt(const U: ITensor2D; var m: TInt);
implementation
end.
Up \ObjAlg