Open Source Projects Europe

--- a/SVDLib2.mod
+++ b/SVDLib2.mod.m2pp
@@ -1,35 +1,86 @@
 IMPLEMENTATION MODULE SVDLib2;
 
+  (*========================================================================*)
+  (* WICHTIG: BITTE NUR DIE DATEI SVDLib2.mod.m2pp EDITIEREN !!!            *)
+  (*========================================================================*)
+  (* This file contains two versions - one calling Modula-2 LibDBlas        *)
+  (* the other the Fortran dblas level 1 routines (LibDBlasF77) in          *)
+  (* procedure PdSVDc                                                       *)
+  (*                                                                        *)
+  (*    m2pp    -D __{BLAS|BLASF77}__ < SVDLib2.mod.m2pp > SVDLib2.mod      *)
+  (*                                                                        *)
+  (*  or                                                                    *)
+  (*                                                                        *)
+  (*    m2pp -b -D __{BLAS|BLASF77}__ < SVDLib2.mod.m2pp > SVDLib2.mod      *)
+  (*                                                                        *)
+  (* to preserve line information from *.m2pp file                          *)
+  (*                                                                        *)
+  (* There is also a debug options included  which prints out some          *)
+  (* additonal but anoying info for the check of some procedures            *)
   (*------------------------------------------------------------------------*)
-  (* Modul zur Singul"arwertzerlegung.                                      *)
+  (* 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.                      *)
   (*------------------------------------------------------------------------*)
   (* Letzte Bearbeitung:                                                    *)
   (*                                                                        *)
+  (*       95, MRi: Erstellen der ersten Version                            *)
   (* 10.01.97, MRi: Durchsicht                                              *)
-  (* 29.09.17, MRi: Herausloesen aus LinLib                                 *)
+  (* 29.09.17, MRi: Herausloesen aus LinLib (PSVD,PSVDSolv,POrderSVD)       *)
+  (* 07.12.17, MRi: Erstellen der ersten Verion von PdSVDc aus dSVDc        *)
+  (* 01.07.18, MRi: Zusammenfuehren aller Routinen mit dynamischer PMATRIX  *)
+  (*                in ein Modul, E und S in PdSVDc auf offene Felder       *)
+  (*                umgestellt.                                             *)
   (*------------------------------------------------------------------------*)
   (* Offene Punkte:                                                         *)
   (*                                                                        *)
-  (* - Prozedur OrderSVD auf schnelleren Sortieralgorithmus umstellen       *)
+  (* - Weiteres testen von PdSVDc                                           *)
+  (* - Parameter m,n und Form der Speicherung in A durchsehen, einheitlich  *)
+  (*   machen und dokumentieren !!!                                         *)
   (*------------------------------------------------------------------------*)
-  (* Licence : GNU Lesser General Public License (LGPL)                     *)
+  (* Tests for PdSVD & PdSVDc are given in TstSVD1u2a                       *)
   (*------------------------------------------------------------------------*)
-
-  (* $Id: SVDLib2.mod,v 1.1 2018/01/16 09:20:42 mriedl Exp mriedl $ *)
-
-FROM SYSTEM       IMPORT TSIZE,ADR;
-FROM Storage      IMPORT ALLOCATE,DEALLOCATE;
-FROM Deklera      IMPORT PMATRIX,PVEKTOR;
-FROM Errors       IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
-FROM LongMath     IMPORT sqrt;
-FROM LMathLib     IMPORT MaxCard,Pythag;
-FROM MatLib       IMPORT AbsSumVek,SumVek;
-
-
-PROCEDURE SVD(VAR A   : PMATRIX; (* Dynamische Matrix A[1..n,1..m] *)
-                  m,n : INTEGER;
-              VAR W   : ARRAY OF LONGREAL;
-              VAR V   : PMATRIX);
+  (* Implementation : Michael Riedl                                         *)
+  (* Licence        : GNU Lesser General Public License (LGPL)              *)
+  (*------------------------------------------------------------------------*)
+
+  (* $Id: SVDLib2.mod.m2pp,v 1.3 2018/08/25 15:04:22 mriedl Exp mriedl $ *)
+
+
+FROM SYSTEM      IMPORT TSIZE;
+FROM Storage     IMPORT ALLOCATE,DEALLOCATE;
+FROM Deklera     IMPORT PMATRIX,PVEKTOR;
+FROM Errors      IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
+                 IMPORT Errors;
+FROM LongMath    IMPORT sqrt;
+FROM LMathLib    IMPORT MaxCard,Pythag;
+FROM MatLib      IMPORT AbsSumVek,SumVek;
+<* IF (__debug__) THEN *>
+                 IMPORT TIO;
+                 IMPORT MatEA,FileSystem;
+<* END *> 
+<* IF (__BLAS__) THEN *>
+FROM LibDBlas    IMPORT daxpy,ddot,dnrm2,drot,drotg,dscal,dswap;
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+FROM LibDBlasF77 IMPORT daxpy,ddot,dnrm2,drot,drotg,dscal,dswap;
+<* END *> 
+
+<* IF (__BLAS__) THEN *>
+CONST tag = "SVDLib2 Modula-2 BLAS";
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+CONST tag = "SVDLib2 Fortran BLAS";
+<* END *> 
+
+PROCEDURE PdSVD(VAR A   : PMATRIX; (* Dynamische Matrix A[1..n,1..m] *)
+                    m,n : INTEGER;
+                VAR W   : ARRAY OF LONGREAL;
+                VAR V   : PMATRIX);
 
   PROCEDURE sign(x,y: LONGREAL): LONGREAL;
   
@@ -284,15 +335,15 @@
         END;
       END;
  *)
-END SVD;
-
-PROCEDURE SVDSolv(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;
-                      m,n : CARDINAL);
+END PdSVD;
+
+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;
+                       m,n : CARDINAL);
 
           VAR  jj,j,i : CARDINAL;
                s      : LONGREAL;
@@ -302,7 +353,7 @@
       ALLOCATE(Zw,n*TSIZE(LONGREAL));
       ALLOCATE(Svek,MaxCard(n,m)*TSIZE(LONGREAL));
       IF (Zw = NIL) OR (Svek = NIL) THEN
-        FatalError("Kein Freispeicher vorhanden (SVDLib2.SVDSolv) !");
+        FatalError("Kein Freispeicher vorhanden (SVDLib1.SVDSolv) !");
       END;
       FOR j:=1 TO n DO
         s := 0.0;
@@ -324,12 +375,12 @@
       END;
       DEALLOCATE(Svek,MaxCard(n,m)*TSIZE(LONGREAL));
       DEALLOCATE(Zw,n*TSIZE(LONGREAL));
-END SVDSolv;
-
-PROCEDURE OrderSVD(VAR U    : PMATRIX;
-                   VAR W    : ARRAY OF LONGREAL;
-                   VAR V    : PMATRIX;
-                       m,n  : CARDINAL);
+END PSVDSolv;
+
+PROCEDURE POrderSVD(VAR U    : PMATRIX;
+                    VAR W    : ARRAY OF LONGREAL;
+                    VAR V    : PMATRIX;
+                        m,n  : CARDINAL);
 
           VAR i,j,MaxI : CARDINAL;
               zw,MaxW  : LONGREAL;
@@ -353,6 +404,585 @@
           END;
         END; (* IF *)
       END;
-END OrderSVD;
-
+END POrderSVD;
+
+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                    *)
+          (*                                                               *)
+<* IF (__BLAS__) THEN *>
+          (* Version basierend auf Modula-2 BLAS routinen                  *)
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+          (* Version basierend auf Fortran BLAS routinen                   *)
+<* END *> 
+          (*---------------------------------------------------------------*)
+
+          PROCEDURE dsign(a,b : LONGREAL) : LONGREAL;
+          BEGIN
+            IF (b >= 0.0) THEN RETURN ABS(a); ELSE RETURN -ABS(a); END;
+          END dsign;
+
+          PROCEDURE dmax(a,b : LONGREAL) : LONGREAL;
+          BEGIN
+            IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
+          END dmax;
+
+          PROCEDURE imax(a,b : INTEGER) : INTEGER;
+          BEGIN
+            IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
+          END imax;
+
+          PROCEDURE imin(a,b : INTEGER) : INTEGER;
+          BEGIN
+            IF (a <= b) THEN RETURN a; ELSE RETURN b; END;
+          END imin;
+
+          CONST maxit     = 30;
+
+          TYPE  WorkSpace =  POINTER TO ARRAY [1..MAX(INTEGER)-1]  OF LONGREAL;
+                PVEKTOR   =  WorkSpace;
+
+          VAR   b,c,cs,el,emm1,f,g                    : LONGREAL;
+                i,j,k,l,kk,ll,lls,ls,lu               : INTEGER;
+                lp1                                   : INTEGER;
+                iter,jobu,kase                        : INTEGER;
+                mm,mm1,mn,mp1,nct,nctp1,ncu,nrt,nrtp1 : INTEGER;
+                scale,shift,sl,sm,smm1,sn,t,t1,test   : LONGREAL;
+                wantu,wantv                           : BOOLEAN;
+                ztest                                 : LONGREAL;
+                work                                  : WorkSpace;
+<* IF (__BLASF77__) THEN *>
+                ione,num                              : INTEGER;
+                tmp                                   : LONGREAL;
+<* END *> 
+BEGIN
+      (* Determine what is to be computed.  *)
+
+      wantu := FALSE;
+      wantv := FALSE;
+      jobu := (job MOD 100) DIV 10;
+      IF (jobu > 1) THEN
+        ncu := imin(m,n);
+      ELSE
+        ncu := m;
+      END;
+      IF (jobu # 0) THEN
+        wantu := TRUE;
+      END;
+      IF ((job MOD 10) # 0) THEN
+        wantv := TRUE;
+      END;
+
+      (* Reduce A to bidiagonal form, storing the diagonal elements *)
+      (* in S and the super-diagonal elements in E.                 *)
+
+      info := 0;
+      nct  := imin(m-1,n);
+      nrt  := imax(0,imin(n-2,m));
+      lu   := imax(nct, nrt);
+
+      ALLOCATE(work,m*TSIZE(LONGREAL));
+
+      IF (work = NIL) THEN
+        Errors.Fehlerflag:="Kein Freispeicher vorhanden (SVDLib4.dSVDc)";
+        Errors.Fehler:=TRUE;
+        Errors.ErrOut(Errors.Fehlerflag);
+        info:=-1;
+        RETURN;
+      END;
+
+<* IF (__BLASF77__) THEN *>
+      ione:=1; (* Fuer F77 Blas routinen *)
+<* END *> 
+
+      FOR l:=1 TO lu DO
+        (* Compute the transformation for the L-th column and *)
+        (* place the L-th diagonal in S(L).                   *)
+        IF (l <= nct) THEN
+<* IF (__BLAS__) THEN *>
+          S[l-1] := dnrm2(m-l+1,A^[l]^[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+          num:=m-l+1;
+          S[l-1] := dnrm2(num,A^[l]^[l],ione);
+<* END *> 
+          IF (S[l-1] # 0.0) THEN
+            IF (A^[l]^[l] # 0.0) THEN
+              S[l-1] := dsign(S[l-1],A^[l]^[l]);
+            END;
+<* IF (__BLAS__) THEN *>
+            dscal(m-l+1,(1.0/S[l-1]),A^[l]^[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            num:=m-l+1;
+            tmp:= 1.0 / S[l-1];
+            dscal(num,tmp,A^[l]^[l],ione);
+<* END *> 
+            A^[l]^[l]:=A^[l]^[l] + 1.0;
+          END;
+          S[l-1] := -S[l-1];
+        END;
+        FOR j:=l+1 TO n DO
+          (* Apply the transformation.  *)
+          IF (l <= nct) AND (S[l-1] # 0.0) THEN
+<* IF (__BLAS__) THEN *>
+            t := - ddot(m-l+1,A^[l]^[l],1,A^[j]^[l],1) / A^[l]^[l];
+            daxpy(m-l+1,t,A^[l]^[l],1,A^[j]^[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            num:=m-l+1;
+(*          Korrektur (j-1 gegen j) MRi, 21.08.16
+            t := - ddot(num,A^[l]^[l],ione,A^[j-1]^[l],ione) / A^[l]^[l];
+            daxpy(num,t,A^[l]^[l],ione,A^[j-1]^[l],ione);
+ *)
+            t := - ddot(num,A^[l]^[l],ione,A^[j]^[l],ione) / A^[l]^[l];
+            daxpy(num,t,A^[l]^[l],ione,A^[j]^[l],ione);
+<* END *> 
+          END;
+          (* Palace the L-th row of A into E for the           *)
+          (* subsequent calculation of the row transformation. *)
+          E[j-1] := A^[j]^[l];
+        END; (* FOR j *)
+
+        (* Place the transformation in U for *)
+        (* subsequent back multiplication.   *)
+
+        IF (wantu) AND (l <= nct) THEN
+          FOR i:=l TO m DO
+            U^[l]^[i] := A^[l]^[i];
+          END;
+        END;
+        IF (l <= nrt) THEN
+          (* Compute the l-th row transformation and *)
+          (* place the l-th superdiagonal in E(l).   *)
+<* IF (__BLAS__) THEN *>
+          E[l-1] := dnrm2(n-l,E[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+          num:=n-l;
+          E[l-1] := dnrm2(num,E[l],ione);
+<* END *> 
+          IF (E[l-1] # 0.0) THEN
+            IF (E[l] # 0.0) THEN
+              E[l-1] := dsign(E[l-1],E[l]);
+            END;
+<* IF (__BLAS__) THEN *>
+            dscal(n-l,(1.0 / E[l-1]),E[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            num:=n-l;
+            tmp:=1.0 / E[l-1];
+            dscal(num,tmp,E[l],ione);
+<* END *> 
+            E[l]:=E[l] + 1.0;
+          END;
+          E[l-1] := -E[l-1];
+
+          IF (l+1 <= m) AND (E[l-1] # 0.0) THEN
+            (* Apply the transformation *)
+            FOR i:=l+1 TO m DO
+              work^[i] := 0.0;
+            END;
+            FOR j:=l+1 TO n DO
+<* IF (__BLAS__) THEN *>
+              daxpy(m-l,E[j-1],A^[j]^[l+1],1,work^[l+1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              num:=m-l;
+              daxpy(num,E[j-1],A^[j]^[l+1],ione,work^[l+1],ione);
+<* END *> 
+            END;
+            FOR j:=l+1 TO n DO
+<* IF (__BLAS__) THEN *>
+              daxpy(m-l,-E[j-1]/E[l],work^[l+1],1,A^[j]^[l+1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              num:=m-l;
+              tmp:=-E[j-1]/E[l];
+              daxpy(num,tmp,work^[l+1],ione,A^[j]^[l+1],ione);
+<* END *> 
+            END;
+          END;
+          (* Place the transformation in V for *)
+          (* subsequent back multiplication.   *)
+          IF (wantv) THEN
+            FOR i:=l+1 TO n DO
+              V^[l]^[i] := E[i-1];
+            END;
+          END;
+        END; (* IF l *)
+      END; (* FOR l *)
+      DEALLOCATE(work,m*TSIZE(LONGREAL));
+  
+      (* Set up the final bidiagonal matrix of order MN.  *)
+
+      mn := imin(m+1,n);
+      nctp1 := nct + 1;
+      nrtp1 := nrt + 1;
+      IF (nct < n) THEN
+        S[nctp1-1] := A^[nctp1]^[nctp1];
+      END;
+      IF (m < mn) THEN
+        S[mn-1] := 0.0;
+      END;
+      IF (nrtp1 < mn) THEN
+        E[nrtp1-1] := A^[mn]^[nrtp1];
+      END;
+      E[mn-1] := 0.0;
+
+      IF (wantu) THEN (* If required, generate U. *)
+        FOR j:=nctp1 TO ncu DO
+          FOR i:=1 TO m DO
+            U^[j]^[i]:=0.0;
+          END;
+          U^[j]^[j]:=1.0;
+        END;
+(*
+        FOR l:=nct TO 1 BY -1 DO
+ *)
+        FOR ll:=1 TO nct DO
+          l := nct - ll + 1;
+          IF (S[l-1] # 0.0) THEN
+            FOR j:=l+1 TO ncu DO
+<* IF (__BLAS__) THEN *>
+              t := - ddot(m-l+1,U^[l]^[l],1,U^[j]^[l],1) / U^[l]^[l];
+              daxpy(m-l+1,t,U^[l]^[l],1,U^[j]^[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              num:=m-l+1;
+              t := - ddot(num,U^[l]^[l],ione,U^[j]^[l],ione) / U^[l]^[l];
+              daxpy(num,t,U^[l]^[l],ione,U^[j]^[l],ione);
+<* END *> 
+            END;
+<* IF (__BLAS__) THEN *>
+            dscal(m-l+1,-1.0,U^[l]^[l],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            num:=m-l+1;
+            tmp:=-1.0;
+            dscal(num,tmp,U^[l]^[l],ione);
+<* END *> 
+            U^[l]^[l]:=U^[l]^[l] + 1.0;
+            FOR i:=1 TO l-1 DO
+              U^[l]^[i]:=0.0;
+            END;
+          ELSE
+            FOR i:=1 TO m DO
+              U^[l]^[i]:=0.0;
+            END;
+            U^[l]^[l]:=1.0;
+          END;
+        END; (* FOR ll *)
+      END; (* IF wantu *)
+
+      IF (wantv) THEN (* If it is required, generate V.  *)
+(*
+        FOR l:=n TO 1 BY -1 DO
+ *)
+        FOR ll:=1 TO n DO
+          l   := n - ll + 1;
+          lp1 := l + 1;
+          IF (l <= nrt) AND (E[l-1] # 0.0) THEN
+            FOR j:=lp1 TO n DO
+<* IF (__BLAS__) THEN *>
+              t := - ddot(n-l,V^[l]^[lp1],1,V^[j]^[lp1],1) / V^[l]^[lp1];
+              daxpy(n-l,t,V^[l]^[lp1],1,V^[j]^[lp1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              num:=n-l;
+              t := - ddot(num,V^[l]^[lp1],ione,V^[j]^[lp1],ione) / V^[l]^[lp1];
+              daxpy(num,t,V^[l]^[lp1],ione,V^[j]^[lp1],ione);
+<* END *> 
+            END;
+          END;
+          FOR i:=1 TO n DO
+            V^[l]^[i] := 0.0;
+          END;
+          V^[l]^[l] := 1.0;
+        END;
+      END; (* IF wantu *)
+
+<* IF (__debug__  ) THEN *>
+  TIO.WrLn; TIO.WrStr("dSVDc mit MATRIX/VEKTOR"); TIO.WrLn; TIO.WrLn;
+
+  TIO.WrStr("A +"); MatEA.MATaus(FileSystem.Output,mn,mn,A,80,8,3);
+
+  TIO.WrLn; TIO.WrStr("E_i,S_i  i = 1,n");
+  TIO.WrLn; TIO.WrLn;
+  FOR l:=1 TO n DO
+    TIO.WrInt(l,3);
+    TIO.WrLngReal(E[l-1],12,4); TIO.WrLngReal(S[l-1],12,4); TIO.WrLn;
+  END;
+  TIO.WrLn;
+
+  mm:=imax(m,n)+1;
+  TIO.WrStr("U"); MatEA.MATaus(FileSystem.Output,mm,mm,U,80,8,3);
+  TIO.WrStr("V"); MatEA.MATaus(FileSystem.Output,mm,mm,V,80,8,3);
+<* END *> 
+
+      (* Main iteration loop for the singular values. *)
+
+      mm := mn;
+      iter := 0;
+      WHILE (mn # 0) AND (iter <= maxit) DO
+
+        (* This section of the program inspects for negligible elements  *)
+        (* in the S and E arrays. On completion the variables KASE and L *)
+        (* are set as follows:                                           *)
+        (*                                                               *)
+        (* KASE = 1     if S(MN) and E(L-1) are negligible and L < MN    *)
+        (* KASE = 2     if S(L) is negligible and L < MN                 *)
+        (* KASE = 3     if E(L-1) is negligible, L < MN, and             *)
+        (*              S(L), ..., S(MN) are not negligible (QR step).   *)
+        (* KASE = 4     if E(MN-1) is negligible (convergence).          *)
+
+        ll:=1;
+        LOOP
+          l := mn - ll;
+          IF (l = 0) THEN EXIT; END;
+          test := ABS(S[l-1]) + ABS(S[l]);
+          ztest := test + ABS(E[l-1]);
+          IF (ztest = test) THEN
+            E[l-1] := 0.0;
+            EXIT;
+          END;
+          INC(ll);
+        END;
+
+        IF (l = mn-1) THEN
+          kase := 4;
+        ELSE
+          lp1 := l  + 1;
+          mp1 := mn + 1;
+
+          ls := l; (* Wg. Compilerwarnung *)
+          lls:=lp1;
+          LOOP
+            IF (lls > mp1) THEN EXIT; END;
+            ls := mn - lls + lp1;
+            IF (ls = l) THEN
+              EXIT;
+            END;
+            test := 0.0;
+            IF (ls # mn) THEN
+              test:= ABS(E[ls-1]); (* + test *)
+            END;
+            IF (ls # l+1) THEN
+              test := test + ABS(E[ls-2]);
+            END;
+            ztest := test + ABS(S[ls-1]);
+            IF (ztest = test) THEN
+              S[ls-1] := 0.0;
+              EXIT;
+            END;
+            INC(lls);
+          END;
+          IF (ls = l) THEN
+            kase := 3;
+          ELSIF (ls = mn) THEN
+            kase := 1;
+          ELSE
+            kase := 2;
+            l := ls;
+          END;
+        END; (* IF l *)
+        l := l + 1;
+
+        IF (kase = 1) THEN
+          (* Deflate negligible S(MN). *) 
+          mm1     := mn - 1;
+          f       := E[mn-2];
+          E[mn-2] := 0.0;
+(*
+          FOR k:=mm1 TO l BY -1 DO
+ *)
+          FOR kk:=l TO mm1 DO
+            k  := mm1 - kk + l;
+            t1 := S[k-1];
+<* IF (__BLAS__) THEN *>
+            drotg(t1,f,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            drotg(t1,f,cs,sn);
+<* END *> 
+            S[k-1] := t1;
+            IF (k # l) THEN
+              f      := -sn*E[k-2];
+              E[k-2] :=  cs*E[k-2];
+            END;
+            IF (wantv) THEN
+<* IF (__BLAS__) THEN *>
+              drot(n,V^[k]^[1],1,V^[mn]^[1],1,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              drot(n,V^[k]^[1],ione,V^[mn]^[1],ione,cs,sn);
+<* END *> 
+            END;
+          END;
+        ELSIF (kase = 2) THEN
+          (* Split at negligible S(L).  *)
+          f      := E[l-2];
+          E[l-2] := 0.0;
+          FOR k:=l TO mn DO
+            t1   := S[k-1];
+<* IF (__BLAS__) THEN *>
+            drotg(t1,f,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            drotg(t1,f,cs,sn);
+<* END *> 
+            S[k-1] := t1;
+            f      := - sn*E[k-1];
+            E[k-1] :=   cs*E[k-1];
+            IF (wantu) THEN
+<* IF (__BLAS__) THEN *>
+              drot(m,U^[k]^[1],1,U^[l-1]^[1],1,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              drot(m,U^[k]^[1],ione,U^[l-1]^[1],ione,cs,sn);
+<* END *> 
+            END;
+          END;
+        ELSIF (kase = 3) THEN (* Perform one QR step.  *)
+
+          scale := dmax(ABS(S[mn-1]), (* Calculate the shift.  *)
+                     dmax (ABS(S[mn-2]),
+                       dmax (ABS(E[mn-2]),
+                         dmax(ABS(S[l-1]),ABS(E[l-1])
+                   ))));
+
+          sm   := S[mn-1]   / scale;
+          smm1 := S[mn-2] / scale;
+          emm1 := E[mn-2] / scale;
+          sl   := S[l-1]    / scale;
+          el   := E[l-1]    / scale;
+          b  := ((smm1 + sm)*(smm1 - sm) + emm1*emm1) / 2.0;
+          cs := (sm*emm1); c:=cs*cs;
+          shift := 0.0;
+          IF (b # 0.0) OR (c # 0.0) THEN
+(*          shift := sqrt(b*b + c); *)
+            shift := Pythag(b,cs);
+            IF (b < 0.0) THEN
+              shift := -shift;
+            END;
+            shift := c / (b + shift);
+          END;
+          f := (sl + sm)*(sl - sm) - shift;
+          g := sl*el;
+
+          (* Chase zeros *)
+
+          mm1 := mn - 1;
+          FOR k:=l TO mm1 DO
+<* IF (__BLAS__) THEN *>
+            drotg(f,g,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            drotg(f,g,cs,sn);
+<* END *> 
+            IF (k # l) THEN
+              E[k-2] := f;
+            END;
+            f      := cs*S[k-1] + sn*E[k-1];
+            E[k-1] := cs*E[k-1] - sn*S[k-1];
+            g      := sn*S[k];
+            S[k]   := cs*S[k];
+            IF (wantv) THEN
+<* IF (__BLAS__) THEN *>
+              drot(n,V^[k]^[1],1,V^[k+1]^[1],1,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              drot(n,V^[k]^[1],ione,V^[k+1]^[1],ione,cs,sn);
+<* END *> 
+            END;
+<* IF (__BLAS__) THEN *>
+            drotg(f,g,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+            drotg(f,g,cs,sn);
+<* END *> 
+            S[k-1] :=  f;
+            f      :=  cs*E[k-1] + sn*S[k];
+            S[k]   := -sn*E[k-1] + cs*S[k];
+            g      :=  sn*E[k];
+            E[k]   :=  cs*E[k];
+            IF (wantu) AND (k < m) THEN
+<* IF (__BLAS__) THEN *>
+              drot(m,U^[k]^[1],1,U^[k+1]^[1],1,cs,sn);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              drot(m,U^[k]^[1],ione,U^[k+1]^[1],ione,cs,sn);
+<* END *> 
+            END;
+          END;
+          E[mn-2] := f;
+          INC(iter);
+        ELSIF (kase = 4) THEN
+          (* Convergence.  *)
+          IF (S[l-1] < 0.0) THEN
+            (* Make the singular value positive.  *)
+            S[l-1] := -S[l-1];
+            IF (wantv) THEN
+<* IF (__BLAS__) THEN *>
+              dscal(n,-1.0,V^[l]^[1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              tmp:=-1.0;
+              dscal(n,tmp,V^[l]^[1],ione);
+<* END *> 
+            END;
+          END;
+
+          (* Order the singular value.  *)
+
+          WHILE (l # mm) AND (S[l-1] < S[l]) DO
+            t      := S[l-1]; 
+            S[l-1] := S[l]; 
+            S[l]   := t;
+            IF (wantv) AND (l < n) THEN
+<* IF (__BLAS__) THEN *>
+              dswap(n,V^[l]^[1],1,V^[l+1]^[1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              dswap(n,V^[l]^[1],ione,V^[l+1]^[1],ione);
+<* END *> 
+            END;
+            IF (wantu) AND (l < m) THEN
+<* IF (__BLAS__) THEN *>
+              dswap(m,U^[l]^[1],1,U^[l+1]^[1],1);
+<* END *> 
+<* IF (__BLASF77__) THEN *>
+              dswap(m,U^[l]^[1],ione,U^[l+1]^[1],ione);
+<* END *> 
+            END;
+            INC(l);
+          END; (* WHILE*)
+          iter := 0;
+          DEC(mn);
+        END; (* IF kase = 4 *)
+      END; (* WHILE *)
+      IF (iter > maxit) THEN
+        (* If too many iterations have been performed set flag *) 
+        Errors.Fehlerflag:=
+        "Maximale Anzahl der Iterationen ueberschritten (LibSVD.dSVDc)";
+        Errors.Fehler:=TRUE;
+        info := mn;
+      END;
+END PdSVDc;
+
+BEGIN
+      Errors.WriteString(tag);
+      Errors.WriteLn;
 END SVDLib2.