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.
