DEFINITION  MODULE SVDLib2;

  (*------------------------------------------------------------------------*)
  (* Modul zur Singul"arwertzerlegung von dynamisch erzeugten Matrizen      *)
  (* Siehe hierzu auch die Module DynMat und PMatLib.                       *)
  (*                                                                        *)
  (* Modul providing routines to calculate the singual value decomposition  *)
  (* of dynamically allocated Matrices (see modules DynMat and PMatLib)     *)
  (* The SVD routines are based on the Eispack subroutine SVD and LinPack   *)
  (* subroutine dsvdc with pointer MATRIX definitions.                      *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: SVDLib2.def,v 1.3 2018/08/25 15:04:43 mriedl Exp mriedl $ *)

FROM Deklera IMPORT PMATRIX;

PROCEDURE PdSVD(VAR A   : PMATRIX; (* Dynamische Matrix A[1..m,1..n] *)
                    m,n : INTEGER;
                VAR W   : ARRAY OF LONGREAL; (* n *)
                VAR V   : PMATRIX);          (* n,n *) 

          (*----------------------------------------------------------------*)
          (* Berechnet die Sigul"arwertzerlegung von A = U W V^+            *)
          (* wobei U in A zur"uckgegeben wird.                              *)
          (* Dimensionierungen : A[1..m,1..n],W[1..m],V[1..m,1..m].         *)
          (*                                                                *)
          (* Translation of the Eispack SVD and original Algol60 routines   *)
          (*                                                                *)
          (* Lit.: Golub,G.H.; Reinsch,C.; Numer. Math. 14, 403 (1970)      *)
          (*----------------------------------------------------------------*)

PROCEDURE PSVDSolv(VAR U   : PMATRIX;           (* m,n *)
                       Utr : CHAR;
                   VAR W   : ARRAY OF LONGREAL; (* n   *)
                   VAR V   : PMATRIX;           (* n,n *)
                   VAR X   : ARRAY OF LONGREAL; (* n   *)
                   VAR C   : ARRAY OF LONGREAL; (* n   *)
                       m,n : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet das lineare Gleichungssystem A X = C, wobei          *)
          (* A mit der Routine SVD nach A = U W V zerlegt wurde.            *)
          (* Wenn Utr = {t|T} wird angenommen das U^{tr} uebergeben wurde   *)
          (*----------------------------------------------------------------*)

PROCEDURE POrderSVD(VAR U    : PMATRIX;
                    VAR W    : ARRAY OF LONGREAL;
                    VAR V    : PMATRIX;
                        m,n  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Ordnet die SVD-Zerlegung (A = UWV^+) nach Gr"o3e der           *)
          (* Singulaerwerte W.                                              *)
          (*----------------------------------------------------------------*)

PROCEDURE PdSVDc(VAR A    : PMATRIX;
                     m    : INTEGER;
                     n    : INTEGER;
                 VAR S    : ARRAY OF LONGREAL;
                 VAR E    : ARRAY OF LONGREAL;
                 VAR U,V  : PMATRIX;
                     job  : INTEGER;
                 VAR info : INTEGER);

          (*----------------------------------------------------------------*)
          (* Modula-2 Version der LinPack SUBROUTINE DSVDC                  *)
          (*                                                                *)
          (* Anpassung an Modula-2 : Michael Riedl                          *)
          (* Lizenz                : GNU public licence                     *)
          (*                                                                *)
          (* 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,       *)
          (*   C version by                                                 *)
          (*     John Burkardt.                                             *)
          (*                                                                *)
          (* 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, double 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, int LDA, the leading dimension of the array A.        *)
          (*   LDA must be at least M.                                      *)
          (*                                                                *)
          (*   Input, int M, the number of rows of the matrix.              *)
          (*                                                                *)
          (*   Input, int N, the number of columns of the matrix A.         *)
          (*                                                                *)
          (*   Output, double 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, double E[MM], where MM = min(M+1,N), ordinarily      *)
          (*   contains zeros. However see the discussion of INFO for       *)
          (*   exceptions.                                                  *)
          (*                                                                *)
          (*   Output, double 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, int LDU, the leading dimension of the array U.        *)
          (*   LDU must be at least M.                                      *)
          (*                                                                *)
          (*   Output, double 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, int LDV, the leading dimension of the array V.        *)
          (*   LDV must be at least N.                                      *)
          (*                                                                *)
          (*   Input, int 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.            *)
          (*                                                                *)
          (*   Output, int *DSVDC, status indicator INFO.                   *)
          (*   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   *)
          (*----------------------------------------------------------------*)

END SVDLib2.
