DEFINITION MODULE SVDLib3;
(*------------------------------------------------------------------------*)
(* Module provides routines for Tageti factorisation and a complex value *)
(* SVD routine . *)
(* *)
(* zSVDc is a Modula-2 translation of the LinPack F77 subroutine ZSVDC. *)
(* LINPACK is a collection of Fortran subroutines that analyze and solve *)
(* linear equations and linear least-squares probles. The package solves *)
(* linear systems whose matrices are general, banded, symmetric *)
(* indefinite, symmetric positive definite, triangular, and tridiagonal *)
(* square. In addition, the package computes the QR and singular value *)
(* decompositions of rect- angular matrices and applies them to *)
(* least-squares problems LINPACK uses column-oriented algorithms to *)
(* increase efficiency by preserving locality of reference. It was *)
(* developed by Jack Dongarra, Jim Bunch, Cleve Moler and Pete Stewartin *)
(* 1984. *)
(* *)
(* Takagi is a M2 translation of Diag library subroutines Takagi *)
(* The Diag libray is developed and maintained by Thomas Hahn, *)
(* Max Planck Institut fuer Physik http:/www.feyarts.de/diag *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SVDLib3.def,v 1.2 2018/07/28 07:06:10 mriedl Exp mriedl $ *)
PROCEDURE zSVDc(VAR X : ARRAY OF ARRAY OF LONGCOMPLEX;
N,P : INTEGER;
VAR S : ARRAY OF LONGCOMPLEX;
VAR E : ARRAY OF LONGCOMPLEX;
VAR U : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR V : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR Work : ARRAY OF LONGCOMPLEX;
Job : INTEGER;
VAR Info : INTEGER);
(*----------------------------------------------------------------*)
(* Translation of LinPack subroutine ZSVDC *)
(* *)
(* X[p,n] = U[p,p] x S[p,n] x V[n,n]^* (Fortran) *)
(* *)
(* The dimension of X is changed from the Fortran original so *)
(* that X[0..P-1,0..N-1] *)
(* *)
(* X[p,n] = V[p,p] x S[p,n] x U[n,n]^* (Modula-2) *)
(* *)
(* HINWEIS: Pruefen ob das alles so in Ordnung ist ... *)
(* ======== *)
(* *)
(* zSVDc is a procedure to reduce a complex N x P matrix X by *)
(* unitary transformations U and V to diagonal form. The *)
(* diagonal elements s(i) are the singular values of X. The *)
(* columns of U are the corresponding left singular vectors, *)
(* and the columns of V the right singular vectors. *)
(* *)
(* On entry *)
(* *)
(* X : X(ldx,p), where ldx >= N *)
(* X contains the matrix whose singular value *)
(* decomposition is to be computed. X is destroyed *)
(* by zSVDc *)
(* ldx : ldx is the leading dimension of the array X *)
(* N : the number of rows of the matrix X *)
(* P : P is the number of columns of the matrix X *)
(* Work : Work(N) is a scratch array *)
(* Job : job controls the computation of the singular *)
(* vectors. It has the decimal expansion ab *)
(* with the following meaning *)
(* a = 0 : do not compute the left singular *)
(* vectors. *)
(* a = 1 : return the n left singular vectors *)
(* in U. *)
(* a > 1 : returns the first min(N,P) left *)
(* singular vectors in U. *)
(* b = 0 : do not compute the right singular *)
(* vectors. *)
(* b = 1 : return the right singular vectors *)
(* in V. *)
(* On return *)
(* *)
(* S : S(mm), where mm=min(N+1,P). *)
(* the first min(n,p) entries of s contain the *)
(* singular values of X arranged in descending *)
(* order of magnitude. *)
(* E : E(P) ordinarily contains zeros. However see the *)
(* discussion of Info for exceptions. *)
(* U : U(ldu,k), where ldu >= N. If joba = 1 then *)
(* k = N, if joba > 1 then k = min(N,P). *)
(* U contains the matrix of left singular vectors. *)
(* U is not referenced if joba = 0. If N <= P or if *)
(* joba > 2, then U may be identified with X in the *)
(* procedure call. *)
(* V : V(ldv,P), where ldv >= P . V contains the matrix *)
(* of right singular vectors. *)
(* V is not referenced if jobb = 0. If P <= N, *)
(* then V may be identified whth X in the *)
(* procedure call. *)
(* Info : The singular values (and their corresponding *)
(* singular vectors) S(info+1),S(info+2),...,S(M) *)
(* are correct (here m=min(N,P)). Thus if info = 0, *)
(* all the singular values and their vectors are *)
(* correct. In any event, the matrix *)
(* B = ctrans(U)*X*V is the bidiagonal matrix with *)
(* the elements of S on its diagonal and the *)
(* elements of E on its super-diagonal (ctrans(U) *)
(* is the conjugate-transpose of U). Thus the *)
(* singular values of X and B are the same. *)
(* *)
(* linpack. This version dated 03/19/79 . *)
(* Correction to shift calculation made 2/85. *)
(* G.W. Stewart, University of Maryland, Argonne National Lab. *)
(* Adopted to Modula-2 by Michael Riedl *)
(*----------------------------------------------------------------*)
PROCEDURE Takagi( N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR D : ARRAY OF LONGREAL;
VAR U : ARRAY OF ARRAY OF LONGCOMPLEX;
sort : INTEGER);
(*----------------------------------------------------------------*)
(* Computes the Takagi factorization of a complex symmetric *)
(* matrix. The code is adapted from the "Handbook" routines *)
(* Wilkinson, Reinsch: Handbook for Automatic Computation, *)
(* pp 202 (Jacobi routine) *)
(* *)
(* The original Fortran code is part of the Diag library by *)
(* T. Hahn (2006) *)
(* *)
(* Adopted to Modula-2 by Michael Riedl *)
(* *)
(* Reference: *)
(* Takagi, T.; Japanese J. Math. 1, pp 83 (1927) *)
(*----------------------------------------------------------------*)
END SVDLib3.