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.
