DEFINITION MODULE SVDLib1;
(*------------------------------------------------------------------------*)
(* Dieses Modul enthaellt SVD routinen fuer realwertige Matrizen *)
(* *)
(* This module provides SVD routines for real matrices *)
(* *)
(* This Module provides two simple SVD routines based on the power *)
(* (Lanczos) method (PowSVD) or on plane (Jacobi) rotations (JacobiSVD). *)
(* Furthermore an adpoted of the least square solver (CalcLQsol) found in *)
(* ref [2] and a translation of the complete least squares solver GivSVD *)
(* *)
(* GivSVD can be used as an alternative to the Golub/Reinsch routine *)
(* MinFit found in SVDLib1 *)
(*------------------------------------------------------------------------*)
(* Literatur / references: *)
(* *)
(* [1] : Nash, John C.; Shlien, Seymour; "Simple Algorithms for the *)
(* Partial Singular Value Decomposition", Computer J., 30,3 (1987) *)
(* [2] : Nash, John C.; "Compact numerical methods for computers: *)
(* linear algebra and function minimisation", 2. Edition, *)
(* Adam Hilger, Bristol, UK (1990) *)
(* [3] : Golub, Gene H.; Van Loan, Charles F.; *)
(* "Matrix Computations", Third Ed., J. Hopkins Univ. Press, *)
(* Baltimore, US (1996) - Chapter 9 "Lanczos Methods" *)
(* *)
(* *)
(*------------------------------------------------------------------------*)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SVDLib1.def,v 1.4 2018/08/26 14:25:52 mriedl Exp mriedl $ *)
PROCEDURE PowSVD( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR sv : ARRAY OF LONGREAL;
VAR U,V : ARRAY OF ARRAY OF LONGREAL;
klimit : INTEGER;
VAR rank : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Calculates the SVD of matrix A by power iterations *)
(* This routine is intended to be used if not the full rank SVD *)
(* is needed and the singular values are clearly separated *)
(* *)
(* M : Number of equations *)
(* N : Number of unknowns *)
(* A : M x N Matrix for which the SVD shall be computed *)
(* sv : Singular values *)
(* U : M x N Matrix of left hand vectors *)
(* tol : threshold for regaring a singular value as zero *)
(* klimit : on entry the required number of singular values *)
(* rank : on exit the number of calculated singular values *)
(* iFehl : Error indicator *)
(* 0 : No error occured *)
(* 1 : matrix A is rank deficite (see klimit) *)
(* 2 : no convergence whilst calculating one SVD plane *)
(* 3 : Both error 1 & 2 occured *)
(* *)
(* Please note that the order for storing the matrix elements is *)
(* "uncommon". One equation is stored row-wise, not column-wise *)
(* (if we refer to a set of linear equations, par example). *)
(* That way we gain speed in the calculations. *)
(*----------------------------------------------------------------*)
PROCEDURE JacobiSVD( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
job : CARDINAL;
VAR rank : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Calculates the SVD of matrix A by Jacobi iterations *)
(* *)
(* M : number of equations *)
(* N : number of unknowns *)
(* A : on entry the N x M Matrix for which the SVD shall be *)
(* computed *)
(* on output the N x M Matrix of left hand vectors U *)
(* W : Singular values *)
(* V : N x N Matrix of right hand vectors *)
(* job : 0 - neither left hand V nor right hand vectors are *)
(* calculated *)
(* 1 - only left hand vectors V are calculated *)
(* 2 - only right hand vectors U are calculated *)
(* 3 - both right and left hand vectors are calculated *)
(* iFehl : Error indicator *)
(* 0 : No error occured *)
(* 1 : matrix A is rank deficite (see klimit) *)
(* 2 : too many sweeps for calculating one SVD plane *)
(* 3 : Both error 1 & 2 occured *)
(* *)
(* If V is not needed you can give the matrix A as argument to V *)
(* to fill up the argument list *)
(* *)
(* Please note that the order for storing the matrix elements is *)
(* "uncommon". One equation is stored row-wise, not column-wise *)
(* (if we refer to a set of linear equations, par example). *)
(* That way we gain speed in the calculations. *)
(* *)
(* This is essentialy the same routine as described by Nash *)
(* and the Pascal source provided (NashSVD) with another *)
(* organization of A[0..N-1,0..M-1] and all loops organized in a *)
(* cache optimal manner. In addition, U is normalized *)
(*----------------------------------------------------------------*)
PROCEDURE CalcLQsol1( M,N : INTEGER;
rank : INTEGER;
IfT : CHAR;
VAR U : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
tol : LONGREAL);
(*----------------------------------------------------------------*)
(* Least squares solution via singular value decomposition. *)
(* On entry, U and V must have the working matrix resulting from *)
(* the operation of JacobiSVD on a real matrix A with job=3. *)
(* Note that A could be omitted if residuals were not wanted. *)
(* However, the user would then lose the ability to interact with *)
(* the problem by changing the tolerance tol. *)
(* *)
(* M : number of equations *)
(* N : numnber of unknowns *)
(* rank : numnber of unknowns *)
(* IfT : if IfT="{T|t} it is assumed that both U and V are *)
(* to be accessed in transposed order, else classical *)
(* order as used e.g. by SVDLib1.dSVD. *)
(* U : the matrix of right hand SVD vectors calculated *)
(* W : the the singular values *)
(* V : the matrix of left hand SVD vectors calculated *)
(* X : on return the vector of (estimated) parameters *)
(* B : the vector to be approximated *)
(* tol : threshold for singular values. Values smaller than *)
(* tol are considered to be zero *)
(* If tol <=0 on entry the default value 100*MachEps *)
(* will be used *)
(* *)
(* Slightly modified and adopted to Modula-2 by M.Riedl *)
(* Please read the notes for organizing the equations made in the *)
(* heading of JacobiSVD (as with the other SVD routines) *)
(* In effect the SVD can also be provided by other routines than *)
(* JacobiSVD. *)
(* *)
(* Copyright 1988 John C. Nash (Pascal procedure SVDlss) *)
(* Published under the LGPL with written permission of the author *)
(*----------------------------------------------------------------*)
PROCEDURE CalcLQres1( M,N : INTEGER;
IfT : CHAR;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR R : ARRAY OF LONGREAL;
VAR r1 : LONGREAL;
VAR r2 : LONGREAL);
(*----------------------------------------------------------------*)
(* Calculate the residuals of the eqations A x X_i = B_i *)
(* *)
(* M : number of equations *)
(* N : numnber of unknowns *)
(* IfT : if IfT="{T|t} it is assumed that A is to be accessed *)
(* in transposed order or classical order as used by *)
(* e.g. SVDLib1.dSVD. *)
(* A : original matrix of coefficients *)
(* X : the vector of (estimated) parameters *)
(* B : the vector to be approximated *)
(* R : calculated residuals *)
(* r1 : calculated sum of residuals *)
(* r2 : calculated sum of squares of residuals devided by M *)
(* *)
(* Getestet mit PowSVD, JacobiSVD, dSVD & dSVDc *)
(* *)
(*----------------------------------------------------------------*)
PROCEDURE GivSVD( M,N : INTEGER;
nRHS : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF ARRAY OF LONGREAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR rss : ARRAY OF LONGREAL;
VAR svs : ARRAY OF LONGREAL;
VAR W : ARRAY OF ARRAY OF LONGREAL;
tol : LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Givens reduction, svd, and least squares solution *)
(* *)
(* List of arguments (bitte Ueberarbeiten !!!) *)
(* ======================= *)
(* *)
(* N : Order of least squares problem *)
(* M : number of obversions *)
(* nRHS : no. of right hand sides *)
(* A : Matrix of coefficients of the linear system *)
(* X : Solution vector(s) *)
(* B : Right hand side of equations (B[1..nRHS,1..M]) *)
(* rss : Calculate uncorrelated residual(s) *)
(* trss : Calculated sum of squares of residuals *)
(* svs : Calculated singular values *)
(* W : On input the matrix of coefficients, on output the *)
(* computed SVD decomposition (together with Z) *)
(* W must be n+1 by n+nRHS in size *)
(* tol : Threshhold for singular values. Values smaller as tol *)
(* are considered to be zero. *)
(* iFehl : Error indicator *)
(* 0 : No error occured *)
(* 1 : matrix A is rank deficite (see klimit) *)
(* 2 : too many sweeps for calculating one SVD plane *)
(* 3 : Both error 1 & 2 occured *)
(* *)
(* Please note that the memory allocation scheme of this routine *)
(* is not changed to be in line with JacobiSVD, PowSVD and dSVDc *)
(* *)
(* A simple change could be to provide a routine which returns *)
(* on line of the linear equation per call as an argument to the *)
(* routine instead of explicitly handing over A and B *)
(* *)
(* Translated to Modula-2 and slighly modified by M.Riedl (2016) *)
(* *)
(* Copyright 1988 J. C. Nash *)
(* *)
(* Published under the GPL with written authorisation of J.C.Nash *)
(*----------------------------------------------------------------*)
(* Explenation of the algorithem (taken from ref [2]) *)
(* *)
(* In this program, which is designed to use a very small *)
(* working array yet solve least squares problems with large *)
(* numbers of observations, we do not explicitly calculate the *)
(* U matrix of the singular value decomposition. *)
(* One could save the rotations and carefully combine them to *)
(* produce the U matrix. However, this algorithm uses plane *)
(* rotations not only to zero elements in the data in the *)
(* Givens��� reduction and to orthogonalize rows of the work array *)
(* in the svd portion of the code, but also to move the data *)
(* into place from the (n+1)st row of the working array into *)
(* which the data is read. *)
(* These movements i.e. of the observation number nobs, *)
(* would normally move the data to row number nobs of the *)
(* original matrix A to be decomposed. *)
(* However, it is possible, as in the array given by the *)
(* following example with 3 columns in the matrix A and 1 right *)
(* hand side, *)
(* *)
(* Col 1 Col 2 Col 3 RHS *)
(* *)
(* 5.0 0.000001 1.0 1.0 *)
(* 6.0 0.999999 1.0 2.0 *)
(* 7.0 2.000010 1.0 3.0 *)
(* 8.0 2.999900 1.0 4.0 *)
(* *)
(* that this movement does not take place. This is because we *)
(* use a complete cycle of Givens��� rotations using the diagonal *)
(* elements W[i,j], j:=1 to n, of the work array to zero the *)
(* first n elements of row nobs of the (implicit) matrix A. In *)
(* the example, row 1 is rotated from row 4 to row 1 since W is *)
(* originally null. Observation 2 is loaded into row 4 of W, but *)
(* the first Givens��� rotation for this observation will zero the *)
(* first TWO elements because they are the same scalar multiple *)
(* of the corresponding elements of observation 1. Since W[2,2] *)
(* is zero, as is W[4,2], the second Givens��� rotation for *)
(* observation 2 is omitted, whereas we should move the data to *)
(* row 2 of W. Instead, the third and last Givens��� rotation for *)
(* observation 2 zeros element W[4,3] and moves the data to row 3 *)
(* of W. In the least squares problem such permutations are *)
(* irrelevant to the final solution or sum of squared residuals. *)
(* However, we do not want the rotations which are only used to *)
(* move data to be incorporated into U. Unfortunately, as shown *)
(* in the example above, the exact form in which such rotations *)
(* arise is not easy to predict. Therefore, we do not recommend *)
(* that this algorithm be used via the rotations to compute the *)
(* svd unless the process is restructured as in Algorithm 3. *)
(* Note that in any event a large data array is needed. *)
(* The main working matrix W must be n+l by n+nRHS in size. *)
(* *)
(* Please note that the memory organization of Matrix B is *)
(* different to that used in MinFit *)
(*----------------------------------------------------------------*)
PROCEDURE dSVD(VAR A : ARRAY OF ARRAY OF LONGREAL;
m,n : INTEGER;
MatU : CARDINAL;
VAR W : ARRAY OF LONGREAL;
MatV : CARDINAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR IErr : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet die SVD-Zerlegung der Matrix A = U x W x V *)
(* *)
(* A : A[M,N] (eingabeseitig), U[M,N] ausgabeseitig) *)
(* m : Spaltendimension (column) *)
(* n : Zeilendimension (row) *)
(* W : Singulearwerte der SVD-Zerlegung *)
(* MatU : Wenn 1 werden die Linksvektoren der SVD-Zerlegung *)
(* berechnet *)
(* MatV : Wenn 1 werden die Rechtsvektoren der SVD-Zerlegung *)
(* berechnet *)
(* V : V[N,N], Rechtsvektoren der SVD-Zerlegung *)
(* IErr : Fehlercode *)
(* > 1 : Der Singulaerwert mit der Ordnungszahl *)
(* IErr+1 ist nicht konvergiert, Werte und *)
(* Vektoren kleinergleich IErr sollten in *)
(* Ordnung sein *)
(* Errors.Fehlerflag wird entsprechend gesetzt *)
(* = 0 : Alles in Ordnung *)
(* = -1 : m < n, Errors.Fehlerflag entsprechend *)
(* gesetzt *)
(* < -1 : Wie IErr > 1, zusaetzlich ist m < n *)
(* *)
(* Computation of the singular values and complete orthogonal *)
(* decomposition of a real rectangular matrix A, *)
(* *)
(* A = U diag(W) V^T, U^T U = V^T V = I, *)
(* *)
(* where the arrays A[1:m,1:n], U[1:m,1:n], V[1:n,1:n], W[1:n] *)
(* represent A, U, V, W respectively. The actual parameters *)
(* corresponding to A (U) and V may all be identical unless *)
(* MatU = MatV = 1. In this case, the actual parameters *)
(* corresponding to A/U and V must differ. m >= n is assumed *)
(* *)
(* Takes an mxn matrix A and decomposes it into UxWxV, where U *)
(* and V are left and right orthogonal transformation matrices, *)
(* and W is a diagonal matrix of singular values. *)
(* *)
(* Input to dsvd is as follows: *)
(* *)
(* A = m x n matrix to be decomposed. *)
(* Contains the matrix U (orthogonal column vectors) of the *)
(* decomposition if MatU has been set to 1 *)
(* m = row dimension of A *)
(* n = column dimension of A *)
(* W = returns the vector of singular values of A *)
(* V = returns the right orthogonal transformation matrix *)
(* if MatV has been set to 1 *)
(* *)
(* Biggest partes of the code originate form Algol 60 procedure *)
(* svd published in Golub, G.H.; Reinsch, C.; Num. Math. 14, *)
(* 403 (1970) respectivly the Fortran 77 adaption in Eispack *)
(*----------------------------------------------------------------*)
PROCEDURE OrderSVD(VAR U : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
m,n : CARDINAL);
(*----------------------------------------------------------------*)
(* Ordnet die SVD-Zerlegung (A = UWV^+) nach Gr"o3e der *)
(* Singulaerwerte W. Geeignet fuer dSVD und dSVDc (Test OB) *)
(* *)
(* Orders the SVD factorisation accoding to the size of the *)
(* singular values. Can be used for dSVD and dSVDc *)
(*----------------------------------------------------------------*)
PROCEDURE SVDSolv(VAR U : ARRAY OF ARRAY OF LONGREAL; (* m,n *)
Utr : CHAR;
VAR W : ARRAY OF LONGREAL; (* n *)
VAR V : ARRAY OF ARRAY OF LONGREAL; (* n,n *)
VAR X : ARRAY OF LONGREAL; (* n *)
VAR C : ARRAY OF LONGREAL;
M,N : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet das lineare Gleichungssystem A X = C, wobei *)
(* A mit der Routine dSVD nach A = U W V zerlegt wurde. *)
(* Wenn Utr = {t|T} wird abgeommen das U^{tr} uebergeben wurde *)
(* Geeignet fuer PowSVD, JacobiSVD, dSVD und dSVDc (Tests OB) *)
(* *)
(* U : Linksvektoren der SVD-Zerlegung (V[1..M,1..N]) *)
(* Utr : Wenn {"T"|"t"} wird transponiert auf U zugegriffen *)
(* W : Singulaerwerte (W[1..N]) *)
(* V : Rechtsvektoren der SVD-Zerlegung (V[1..N,1..N]) *)
(* X : Gesuchter Loesungsvektor (X[1..N]) *)
(* C : Rechte Seite des Gleichungssystems (C[1..M]) *)
(* M : Anzahl der Gleichungen *)
(* N : Anzahl der Unbekannten *)
(* *)
(* Calculates the solution of the linear system A x X = C where *)
(* A had been SVD decomposited according to A = U x W x V *)
(* If CAP(Utr) = T then U is accessed in transposed order *)
(*----------------------------------------------------------------*)
PROCEDURE MinFit( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
Nb : INTEGER;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* This is a translation of the algol procedure minfit *)
(* num. math. 14, 403-420(1970) by golub and reinsch. *)
(* handbook for auto. comp., vol ii-linear algebra, 134-151(1971) *)
(* *)
(* This subroutine determines, towards the solution of the linear *)
(* t *)
(* system ax=b, the singular value decomposition a=usv of a *)
(* t *)
(* real m by n rectangular matrix, forming u b rather than u. *)
(* householder bidiagonalization and a variant of the QR *)
(* algorithm are used. *)
(* *)
(* on input *)
(* *)
(* nm must be set to the row dimension of two-dimensional *)
(* array parameters as declared in the calling program *)
(* dimension statement. note that nm must be at least *)
(* as large as the maximum of m and n. *)
(* M is the number of rows of a and b. *)
(* N is the number of columns of a and the order of v. *)
(* A contains the rectangular coefficient matrix of the *)
(* system *)
(* Ip is the number of columns of b. ip can be zero. *)
(* B contains the constant column matrix of the system *)
(* if ip is not zero. otherwise b is not referenced. *)
(* *)
(* on output *)
(* *)
(* A has been overwritten by the matrix v (orthogonal) of *)
(* the decomposition in its first n rows and columns. if *)
(* an error exit is made, the columns of v corresponding *)
(* to indices of correct singular values should be correct. *)
(* 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. *)
(* t *)
(* B has been overwritten by u b. if an error exit is made, *)
(* t *)
(* the rows of u b corresponding to indices of correct *)
(* singular values should be correct. *)
(* ierr is set to *)
(* zero for normal return, *)
(* k if the k-th singular value has not been *)
(* determined after 30 iterations. *)
(* rv1 is a temporary storage array. *)
(* *)
(* this version dated august 1983. *)
(* *)
(* Please not that the memory organization for matrix B is *)
(* B[1..M,1..IP] *)
(* *)
(* Adaption nach Modula-2 von Michael Riedl, July 2016 *)
(*----------------------------------------------------------------*)
PROCEDURE BerLsg(VAR M,N : INTEGER;
IP : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL; (* [1..N,1..N] *)
VAR A0 : ARRAY OF ARRAY OF LONGREAL; (* [1..M,1..N] *)
VAR B,B0 : ARRAY OF ARRAY OF LONGREAL; (* [1..M,IP] *)
VAR X : ARRAY OF LONGREAL; (* [1..N] *)
VAR W : ARRAY OF LONGREAL; (* [1..N] *)
VAR R : ARRAY OF LONGREAL; (* [1..M] *)
VAR ifR : BOOLEAN);
(*----------------------------------------------------------------*)
(* Calculate the solution of the least squares problem A x X = B *)
(* after a call to MinFit. *)
(* *)
(* M : Anzahl der Gleichungen *)
(* N : Anzahl der Unbekannten *)
(* IP : Anzahl der rechten Seiten des Gleichungssystems *)
(* A : Linksvektoren der SVD-Zerlegung wie von MinFit zu- *)
(* rueckgegeben A[1..M,1..N]) *)
(* A0 : Orginale Koeffizientenmatrix der Gleichungssystems zur *)
(* Berechnung der Residuen *)
(* B : Rechte Seite des Gleichungssystems (B[1..M,1..IP]) *)
(* wie von MinFit zurueckgegeben *)
(* B0 : Originale rechte Seite des Gleichungssystems *)
(* fuer die berechnung der Residuen (B[1..M,1..IP]) *)
(* X : Gesuchter Loesungsvektor (X[1..N]) *)
(* W : Singulaerwerte (W[1..N]) wie von MinFit berechnet *)
(* R : Residuen der Loesung (R[1..M]) *)
(* ifR : Wenn wahr werden mithilfe von A0 und B0 die Residuen *)
(* der Loesung berechnet *)
(* *)
(* (Tests OB) *)
(*----------------------------------------------------------------*)
PROCEDURE dSVDc(VAR A : ARRAY OF ARRAY OF LONGREAL;
m : INTEGER;
n : INTEGER;
VAR S : ARRAY OF LONGREAL;
VAR E : ARRAY OF LONGREAL;
VAR U,V : ARRAY OF ARRAY OF LONGREAL;
job : INTEGER;
VAR info : INTEGER);
(*----------------------------------------------------------------*)
(* Modula-2 Version der LinPack SUBROUTINE DSVDC *)
(* *)
(* Berechnet die SVD-Zerlegung der Matrix A = U x W x V *)
(* *)
(* A : A[n,m] *)
(* m : Zeilendimension (row) *)
(* n : Spaltendimension (column) *)
(* S : Singulaerwerte der SVD-Zerlegung *)
(* E : Subdiagonalelemente der Bidiagonalgestalt von A *)
(* U : Linksvektoren der SVD-Zerlegung *)
(* U[M,M] wenn job > 1 *)
(* U[min(M,N),M] wenn job > 1 *)
(* V : V[N,N], Rechtsvektoren der SVD-Zerlegung *)
(* job : kontrolliert die Berechnung der Singulaervektoren *)
(* Job ist als dezimale Erweiterung AB mit der folgenden *)
(* Bedeutung implementiert: *)
(* A = 0 : keine Berechnung der Linksvektoren *)
(* A = 1 : berechne m Linksvektoren in U *)
(* A > 1 : berechne die ersten min(M,N) Linksvektoren *)
(* B = 0 : keine Berechnung der Rechtsvektoren *)
(* B = 1 : berechne n Rechtsvektoren in V *)
(* Also fuer job = 21 werden alle Vektoren berechnet *)
(* info : siehe unten *)
(* *)
(* Purpose: *)
(* *)
(* DSVDC computes the singular value decomposition of a real *)
(* rectangular matrix. *)
(* *)
(* Discussion: *)
(* *)
(* This routine reduces an M by N matrix A to diagonal form by *)
(* orthogonal transformations U and V. The diagonal elements *)
(* S[i] are the singular values of A. The columns of U are the *)
(* corresponding left singular vectors, and the columns of V *)
(* the right singular vectors. *)
(* *)
(* The form of the singular value decomposition is then *)
(* *)
(* A(MxN) = U(MxM) * S(MxN) * V(NxN)' *)
(* *)
(* Licensing: *)
(* *)
(* This code is distributed under the GNU LGPL license. *)
(* *)
(* Modified: *)
(* *)
(* 03 May 2007 *)
(* *)
(* Author: *)
(* *)
(* Original FORTRAN77 version by *)
(* Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, *)
(* *)
(* Reference: *)
(* *)
(* Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, *)
(* LINPACK User's Guide, *)
(* SIAM, (Society for Industrial and Applied Mathematics), *)
(* 3600 University City Science Center, *)
(* Philadelphia, PA, 19104-2688. *)
(* ISBN 0-89871-172-X *)
(* *)
(* Parameters: *)
(* *)
(* Input/output, A[LDA,N]. On input, the M by N matrix *)
(* whose singular value decomposition is to be computed. On *)
(* output, the matrix has been destroyed. Depending on the *)
(* user's requests, the matrix may contain other useful *)
(* information. *)
(* *)
(* Input, LDA, the leading dimension of the array A. *)
(* LDA must be at least M. *)
(* *)
(* Input, M, the number of rows of the matrix. *)
(* *)
(* Input, N, the number of columns of the matrix A. *)
(* *)
(* Output, S[MM], where MM = min(M+1,N). The first *)
(* min(M,N) entries of S contain the singular values of A *)
(* arranged in descending order of magnitude. *)
(* *)
(* Output, E[MM], where MM = min(M+1,N), ordinarily *)
(* contains zeros. However see the discussion of INFO for *)
(* exceptions. *)
(* *)
(* Output, U[LDU,K]. If JOBA = 1 then K = M; *)
(* if 2 <= JOBA, then K = min(M,N). U contains the M by M *)
(* matrix of left singular vectors. U is not referenced if *)
(* JOBA = 0. If M <= N or if JOBA = 2, then *)
(* U may be identified with A in the subroutine call. *)
(* *)
(* Input, LDU, the leading dimension of the array U. *)
(* LDU must be at least M. *)
(* *)
(* Output, V[LDV,N], the N by N matrix of right *)
(* singular vectors. V is not referenced if JOB is 0. *)
(* If N <= M, then V may be identified with A in the *)
(* subroutine call. *)
(* *)
(* Input, LDV, the leading dimension of the array V. *)
(* LDV must be at least N. *)
(* *)
(* Input, 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 M left singular vectors in U. *)
(* A >= 2, return the first min(M,N) singular vectors in U. *)
(* B = 0, do not compute the right singular vectors. *)
(* B = 1, return the right singular vectors in V. *)
(* *)
(* Info : On output status indicator *)
(* The singular values (and their corresponding singular *)
(* vectors) S(INFO+1), S(INFO+2),...,S(MN) are correct. *)
(* Here MN = min ( M, N ). *)
(* Thus if INFO is 0, all the singular values and their *)
(* vectors are correct. In any event, the matrix *)
(* B = U' * A * V is the bidiagonal matrix with the elements *)
(* of S on its diagonal and the elements of E on its *)
(* superdiagonal. Thus the singular values of A and B are *)
(* the same. *)
(* *)
(* Note that the singular values are ordered in decending order *)
(* *)
(* Anpassung an Modula-2 : Michael Riedl *)
(* Lizenz : GNU public licence *)
(*----------------------------------------------------------------*)
END SVDLib1.