Switch to side-by-side view

--- a
+++ b/SVDLib3.def
@@ -0,0 +1,138 @@
+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.