linalg.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Library for Turbo/Object Pascal: Linear Algebra )
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(09.09.1991)
@lastmod(10.10.2002)
This is a temporary publication (reduced variant), will be updated later
㎞ikolai V. Shokhirev, 2002
}
unit  LinAlg;

interface

uses
  MathTypes;

{ Householder (tridiagonal) reduction of a real, symmetric, n by n matrix a.
@created(09.09.1991)
@lastmod(10.10.2002)
  Input: a, n
  Output: a is replased by the orthogonal transformation matrix Q.
         d returns the diaginal elements of the tridiagonal matrix,
         e returns the off-giagonal elements with e[1] = 0. }
procedure tred2(const a: Matrix; n: integer; vectors: boolean;
                         var d, e: Vector);

{@created(09.09.1991)
 @lastmod(10.10.2002)
 QL diagonalization algorithm with implicit shifts for a real tridiaginal
   symmetric, n by n matrix .
 Input: d contains diagonal elements, e - subdiagonal onses
        e[1] is arbitrary.
        Independent use: a is the identity matrix
        Use in combination with tred2: a is the output from tred2
 Output: d contains eigenvalues,
         a is the eigenvector matrix so that the k-th column (a[i,k])
         is the normalized eigenvector corresponding to d[k] }
procedure tqli(var d, e: Vector; n: integer; vectors: boolean;
               const z: Matrix; var signal: integer);

{@created(09.09.1991)
@lastmod(10.10.2002)
 Sorting of the eigenvalues E and eigenvectors C (if vectors = true)
  in ascending (if Ascending = true) or descending (if Ascending = false) order}
procedure Order(var E: Vector; n: integer; vectors, ascending: boolean;
                const C: Matrix );

{@created(09.09.1991)
@lastmod(10.10.2002)
 combines the following calls
  tred2(a, n, vectors,  d, e);
  tqli( d, e, n, vectors, a, signal);
  Order(e, n, vectors, ascending, a);}
procedure Diag(const a: Matrix; n: integer; vectors, ascending: boolean;
               var d, e: Vector; signal: integer);

{ Claccical Jacobi diagonalization of real symmetric square matrices
@lastmod(09.09.1990)
@lastmod(10.10.2002)
  Input:
    N - dimension of the matrices
    A - the matrix to be diagonalized (the upper triangle is used)
  Output:
    Eigen  - the vector of eigen values in increasing order
    T - the matrix of eigenvectors (by columns) corresponding to the Eigenvalues
    Signal - the error key :
         = 0  <=> OK
         = -1 <=> N < 1
         > 0  <=> convergion has not been reached after Signal iterations
  Key parameters:
    ItMax - the maximum available number of iterations
    Eps1 is used at  SNA  and  CSA  calculations
    Eps2 is the level of the elimination of the non-diagonal matrix elements
    Eps3 - the criterion of the iterations to be stopped.
         The iterations stop if (1-Sigma1/Sigma2)<=Eps3 ,
         where  Sigma1 and  Sigma2 are the L2 norms of diagonal
         (or off-diagonal) elements in two consecutive iteration }
procedure Jacobi(N: integer; const A, T: Matrix; var Eigen, Aik: Vector;
                    var Signal: integer);

{ Jacobi diagonalization of Complex Hermitian square matrices
@lastmod(02.02.1994)
@lastmod(10.10.2002)
  Input:
    N - dimension of the matrices
    A + i*B - the matrix to be diagonalized
  Output:
    E  - the vector of eigen values in increasing order
    U + i*V - the matrix of eigenvectors (by columns) corresponding to E
    Signal - the error key :
        = 0  <=> OK
        = -1 <=> N < 1
        > 0  <=> convergion has not been reached after Signal iterations
  Key parameters:
    ItMax - the maximum available number of iterations
    Eps1 is used at  SNA  and  CSA  calculations
    Eps2 is the level of the elimination of the non-diagonal matrix elements
    Eps3 - the criterion of the iterations to be stopped: OffDsq < Eps3*Norm }
procedure JacobiC(N: integer; const A, B, U, V: Matrix; var E: Vector;
                  var Signal: integer);

{ Fast Inversion of the matrix A by the Gauss-Jordan elimination algorithm
@lastmod(02.02.1991)
@lastmod(10.10.2002)
  Input:
    A - the matrix to be inversed.
    N - dimension of the matrix
  Output:
    A - the inversed matrix
    Signal - the error key: = 0 <=> OK
                           <> 0 <=> Signal-1 is the rank of degenerate matix}
procedure FastInverse(N: integer; var A: Matrix; var J0: IVector;
                         var Det: RealType; var Signal: integer);

{ Singular Value Decomposition - translation from FORTRAN
  G.E.Forsythe, M.A.Malcolm, C.B.Moler.
  Computer methods for mathematical computations.
  Prentice-Hall, 1977.
@created(09.09.1990)
@lastmod(10.10.2002)
  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 is set to
                     0  for normal exit,
                     k  if the k-th singular value has not been
                        determined after 30 iterations.
    RV1 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 SVD(NA, M, N: integer; const A, U, V:  Matrix; var W, RV1: Vector;
               MatU, MatV:  boolean; RetCode: integer);

{@created(09.09.1991)
@lastmod(10.10.2002)
 sorting }
procedure OrderSVD(M,N: integer; const U,V: Matrix; var W: Vector;
                   MatU,MatV: boolean );

{@created(09.09.1991)
@lastmod(10.10.2002)
 Perturbated Cholesky Decomposition }
procedure CholDecomp(N: integer; var HDiag: Vector; MaxOff,MachEps: RealType;
                        const L: Matrix; var MaxAdd: RealType);

{@created(09.09.1991)
@lastmod(10.10.2002)
 Cholesky L - Solution  of  L*Y  =  B (for given B) }
procedure LSolve(N: integer; const L: Matrix; var B,Y: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 Cholesky LT - Solution  of LT*X  =  Y (for given  Y) }
procedure LTSolve(N: integer; const L: Matrix; var Y,X: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 Solution of the equation L*LT*S = G by the  Cholesky  method }
procedure ChSolve(N: integer; const L: Matrix; var G,S: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 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(N,p: integer; const H,L: Matrix);

{@created(09.09.1991)
@lastmod(10.10.2002)
 Cholesky L - Solution  of L*Y  =  B (for given  B)
  for the band matrix L[1..p+1,1..N].  }
procedure LBandSolve(N,p: integer; const L: Matrix; var B,Y: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 Cholesky LT - Solution  of LT*X  =  Y  (for given  Y)
  for the  band matrix L[1..p+1,1..N].  }
procedure LTBandSolve(N,p: integer; const L: Matrix; var Y,X: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 Solution of the equation L*LT*S = G by the Cholesky method
  for the band matrix L[1..p+1,1..N].  }
procedure ChBandSolve(N,p: integer; const L: Matrix; var G,S: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 QR Decomposition of  M }
procedure QRDecomp(N: integer; const M: Matrix; var M1,M2: Vector;
                      var Sing: boolean);

{@created(09.09.1991)
@lastmod(10.10.2002)
 This routine solves  R*X = B  for  X ,
   where  R  is placed in the  M  by the QR-Decomp routine ;
   X  is returned in B .}
procedure RSolve(N: integer; const M: Matrix; var M2,B: Vector );

{@created(09.09.1991)
@lastmod(10.10.2002)
 This routine solves the equation  (QR)*X = B,
  where Q and R are placed in the  M  by the QR-Decomp routine;
  the computed X is returned  in B.  }
procedure QRSolve(N: integer; const M: Matrix; var M1,M2,B: Vector );

{@created(10.10.2002)
@lastmod(10.10.2002)
{ Gram-Schmidt Orthogonalization                      }
{ On input:                                           }
{ VT -  row vectors;  Vi[k] = VT[i,k]                 }
{ Nvec - number of vectors to be normalized           }
{ On output:                                          }
{ VT: <Vi|Vi> = 1, <Vi|Vj> = 0, i <> j                }
{ Nvec is the actual number of orthogonalized vectors }
procedure Gram_Schmidt(var Nvec: IntType; N: IntType; var VT: Matrix);
{  ------------------------------------------------------------  }
implementation

uses
  MatrixProc, MatLib_;

{ Householder (tridiagonal) reduction of a real, symmetric, n by n matrix a }
procedure tred2(const a: Matrix; n: integer; vectors: boolean; var d, e: Vector);
var
  i, j, k, l: integer;
  scale, hh, h, g, f: RealType;
begin { of tred2 }
  if n > 1 then
  begin

    for i := n downto 2 do
    begin
      l     := i - 1;
      h     := 0.0;
      scale := 0.0;
      if l > 1 then
      begin
        for k := 1 to l do
        begin
          scale := scale + abs(a[i]^[k]);
        end;{ k }
        
        if scale < MinReal then{ Skip transformation }
        begin
          e[i] := a[i]^[l];
        end  { scale < MinReal }
        else{ transformation }
        begin{ scale > MinReal }

          for k := 1 to l do
          begin
            a[i]^[k] := a[i]^[k]/scale;
            h        := h + sqr(a[i]^[k])
          end;{ k }

          f        :=  a[i]^[l];
          g        := -sign2(sqrt(h), f);
          e[i]     :=  scale*g;
          h        :=  h - f*g;
          a[i]^[l] :=  f - g;
          f        :=  0.0;
          for j := 1 to l do
          begin

            if vectors then
            begin
              a[j]^[i] := a[i]^[j]/h;
            end;
            g := 0.0;
            for k := 1 to j do
            begin
              g := g + a[j]^[k]*a[i]^[k];
            end;{ k }
            if l > j then
            begin
              for k := j + 1 to l do
              begin
                g := g + a[k]^[j]*a[i]^[k];
              end; { k }
            end; { l > j }
            e[j] := g/h;
            f    := f + e[j]*a[i]^[j];

          end;{ j }

          hh := f/(h + h);
          for j := 1 to l do
          begin

            f    := a[i]^[j];
            g    := e[j] -hh*f;
            e[j] := g;
            for k := 1 to j do
            begin
              a[j]^[k] := a[j]^[k] - f*e[k] -g*a[i]^[k];
            end;{ k }

          end;{ j }

        end; { scale > MinReal, transformation }
      end   { l > 1 }
      else
      begin{ l = 1 }
        e[i] := a[i]^[l];
      end; { l = 1 }
      d[i] := h;
    end;{ i }

  end;{ n > 1 }

  d[1] := 0.0;
  e[1] := 0.0;
  for i := 1 to n do
  begin
    if vectors then
    begin
      l := i - 1;
      { if d[i] <> 0.0 }
      if abs(d[i]) > MinReal then
      if i > 0 then
      begin{ This block skipped when i = 1 }

        for j := 1 to l do
        begin

          g := 0.0;
          for k := 1 to l do
          begin
            g := g + a[i]^[k]*a[k]^[j];
          end;{ k }

          for k := 1 to l do
          begin
            a[k]^[j] := a[k]^[j] -g*a[k]^[i];
          end;{ k }

        end;{ j }

      end;{ if d[i] <> 0.0 }

    end; { if vectors }

    d[i] := a[i]^[i];

    if vectors then
    begin
      a[i]^[i] := 1.0;
      if l > 0 then
      begin
        for j := 1 to l do
        begin
          a[i]^[j] := 0.0;
          a[j]^[i] := 0.0;
        end;{ j }

      end;{ l > 0 }

    end;{ if vectors }

  end;{ i }

end;{ of tred2 }

{ QL diagonalization algorithm with implicit shifts for a real tridiaginal }
procedure tqli(var d, e: Vector; n: integer;  vectors: boolean;
               const z: Matrix; var signal: integer);
label
  1, 2, 3;

const
  itermax = 999;

var
  i,  k, l, m, iter: integer;
  s, r, p, g, f, dd, c, b: RealType;
begin { of tqli }
  signal := 0;
  if n > 1 then
  begin
    for i := 2 to n do
    begin
      e[i-1] := e[i];
    end;{ i }
    e[n] := 0.0;

    for l := 1 to n do
    begin
      iter := 0;
1:    for m := l to n - 1 do
      begin
        dd := abs(d[m]) + abs(d[m+1]);
        if(abs(e[m]) + dd = dd) then goto 2;
      end; { m }
      m := n;
2:    if m <> l then
      begin
        if iter = itermax then
        begin
          signal := iter;
          goto 3;
        end;
        iter := iter + 1;
        g := (d[l+1] - d[l])/(2.0*e[l]); { Form shift }
        r := sqrt(sqr(g) + 1.0);
        g := d[m] - d[l] + e[l]/(g + sign2(r,g));
        s := 1.0;
        c := 1.0;
        p := 0.0;
        for i := m - 1 downto l do
        begin
          f := s*e[i];
          b := c*e[i];
          if abs(f) >= abs(g) then
          begin
            c := g/f;
            r := sqrt(sqr(c) + 1.0);
            e[i+1] := f*r;
            s := 1.0/r;
            c := c*s;
          end  { abs(f) >= abs(g) }
          else
          begin{ abs(f) < abs(g) }
            s := f/g;
            r := sqrt(sqr(s) + 1.0);
            e[i+1] := g*r;
            c := 1.0/r;
            s := c*s;
          end; { abs(f) < abs(g) }
          g      := d[i+1] - p;
          r      := (d[i] - g)*s + 2.0*c*b;
          p      := s*r;
          d[i+1] := g + p;
          g      := c*r - b;
          if vectors then
          begin
            for k := 1 to n do
            begin
              f := z[k]^[i+1];
              z[k]^[i+1] := s*z[k]^[i] + c*f;
              z[k]^[i]   := c*z[k]^[i] - s*f;
            end;{ k }
          end; { if vectors }
        end;{ i }
        d[l] := d[l] - p;
        e[l] := g;
        e[m] := 0.0;
        goto 1;
       end; {  m <> l }
    end;{ l }
  end; { n > 1 }
3:;
end;{ of tqli }

{ sorting }
procedure Order(var E: Vector; n: integer; vectors, ascending: boolean;
                const C: Matrix);
var
  i, j, k: IntType;
  t: RealType;
begin
  for i := 1 to (n-1) do
    for j := i+1 to n do
      if    (Ascending and (E[j] < E[i])) or
        (not Ascending and (E[j] > E[i])) then
      begin
        t := E[i];
        E[i] := E[j];
        E[j] := t;
        if vectors then
          for k := 1 to n do
          begin
            t := C[k]^[i];
            C[k]^[i] := C[k]^[j];
            C[k]^[j] := t;
          end;
      end;
end;{ of Order }

procedure Diag(const a: Matrix; n: integer; vectors, ascending: boolean;
               var d, e: Vector; signal: integer);
begin{ of Diag }
  tred2(a, n, vectors,  d, e);
  tqli( d, e, n, vectors, a, signal);
  Order(e, n, vectors, ascending, a);
end;{  of Diag }

procedure Jacobi(N: integer; const A, T: Matrix; var Eigen, Aik: Vector;
                    var Signal: integer);

const
  Eps1 = 6.0e-12;    Eps2  = 1.0e-12;
  Eps3 = 1.0e-12;    ItMax = 9999;

var
  i,   j,   k,   Iter   : integer;
  Sigma1, Sigma2, OffDsq,
  SPQ,    CSA,    SNA,    S,   Q,   P,
  HoldIK, HoldKI        : RealType;

label
  1, 2, 3;

begin