--- a
+++ b/SVDLib3.mod
@@ -0,0 +1,675 @@
+IMPLEMENTATION MODULE SVDLib3;
+
+  (*------------------------------------------------------------------------*)
+  (* Module provides routines for Takagi factorisation and a complex value  *)
+  (* SVD routine.                                                           *)
+  (*                                                                        *)
+  (* 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              *)
+  (*------------------------------------------------------------------------*)
+  (* Letzte Veraenderung                                                    *)
+  (*                                                                        *)
+  (* 09.08.11, ThH: last modified of Fortran orgiginal                      *)
+  (* 14.09.16, MRi: Erstellen der ersten uebersetzbaren Takagi Version      *)
+  (* 15.09.16, MRi: Fehlerkorrektur - Takagi scheint zu funktionieren       *)
+  (* 28.09.16, MRi: Umstellen der Sortierung in SVD auf externe Routine     *)
+  (* 20.06.18, MRi: Erstellen der erstern uebersetzbaren Version von        *)
+  (*                zSVDC und der benoetigten BLAS Routinen                 *)
+  (* 21.06.18, MRi: Korrekturen in dznrm2 und zdotc                         *)
+  (*                Testmatrix wird mit zSVDc korrekt berechnet             *)
+  (*------------------------------------------------------------------------*)
+  (* Testroutinen fuer zSVDc in TstSVDLib3a, fuer Takagi in TstSVDLib3b     *)
+  (*------------------------------------------------------------------------*)
+  (* Offene Punkte                                                          *)
+  (*                                                                        *)
+  (* - Weitere Verbesserung der Indizierung in zSVDc ([i-1] Problem)        *)
+  (*------------------------------------------------------------------------*)
+  (* Implementation : Michael Riedl                                         *)
+  (* Licence        : GNU Lesser General Public License (LGPL)              *)
+  (*------------------------------------------------------------------------*)
+
+  (* $Id: SVDLib3.mod,v 1.3 2018/07/28 07:06:02 mriedl Exp mriedl $ *)
+
+FROM Storage  IMPORT ALLOCATE,DEALLOCATE;
+FROM Deklera  IMPORT SIZELONGCMPLX;
+              IMPORT Errors;
+FROM LongMath IMPORT sqrt;
+              IMPORT LongComplexMath;
+FROM LMathLib IMPORT MachEps,CardPot,sign2;
+FROM LibDBlas IMPORT drotg,zswap,zdotc,dznrm2,zscal,zaxpy,zdrot;
+
+CONST CABS       = LongComplexMath.abs;
+      conj       = LongComplexMath.conj;
+      scalarMult = LongComplexMath.scalarMult;
+
+      zero       = LongComplexMath.zero;
+      one        = LongComplexMath.one;
+
+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);
+
+          (* X = U * s * conj(tran(V))                                      *)
+          (* X[P,N], V[P,P], U[N,N]                                         *)
+          (* Benutzt BLAS zaxpy,zdotc,zscal,zswap,dznrm2,drotg              *)
+
+          PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
+          
+          PROCEDURE MAX0(a,b : CARDINAL) : CARDINAL;
+          
+          BEGIN
+                IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
+          END MAX0;
+          
+          PROCEDURE MIN0(a,b : CARDINAL) : CARDINAL;
+          
+          BEGIN
+                IF (a >= b) THEN RETURN b; ELSE RETURN a; END;
+          END MIN0;
+          
+          PROCEDURE DMAX1(a,b : LONGREAL) : LONGREAL;
+          
+          BEGIN
+                IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
+          END DMAX1;
+          
+          PROCEDURE CABS1(x : LONGCOMPLEX) : LONGREAL;
+          
+          BEGIN
+                RETURN ABS(RE(x)) + ABS(IM(x));
+          END CABS1;
+          
+          PROCEDURE CSIGN(a,b : LONGCOMPLEX) : LONGCOMPLEX;
+          
+          BEGIN
+                RETURN scalarMult(1.0/CABS(b), scalarMult(CABS(a),b));
+          END CSIGN;
+          CONST maxit = 30;
+
+          VAR i,j,k,l,m,kk,ll,mm                    : INTEGER;
+              iter,jobu,kase,nctp1,nrtp1            : INTEGER;
+              lls,lm1,lp1,ls,lu,mm1,mp1,nct,ncu,nrt : INTEGER;
+              t,r                                   : LONGCOMPLEX;
+              b,c,cs,sn,el,emm1,f,g                 : LONGREAL;
+              scale,shift,test,ztest                : LONGREAL;
+              sl,sm,smm1,t1                         : LONGREAL;
+              wantu,wantv                           : BOOLEAN;
+BEGIN
+      (* determine what is to be computed. *)
+
+      wantu := FALSE;
+      wantv := FALSE;
+      jobu := (Job MOD 100) DIV 10;
+      ncu := N;
+      IF (jobu > 1) THEN
+        IF (N <= P) THEN ncu:=N; ELSE ncu:=P; END;
+      END;
+      IF (jobu # 0) THEN
+        wantu := TRUE;
+      END;
+      IF ((Job MOD 10) # 0) THEN
+        wantv := TRUE;
+      END;
+
+      (* reduce x to bidiagonal form, storing the diagonal   *)
+      (* elements in s and the super-diagonal elements in e. *)
+
+      Info := 0;
+      nct  := MIN0(N-1,P);
+      nrt  := MAX0(0,MIN0(P-2,N));
+      lu   := MAX0(nct,nrt);
+      IF (lu >= 1) THEN
+       FOR l:=1 TO lu DO
+          lp1 := l + 1;
+          IF (l <= nct) THEN
+            (* compute the transformation for the l-th column *)
+            (* and place the l-th diagonal in s(l).           *)
+            S[l-1] := CMPLX(dznrm2(N-l+1,X[l-1,l-1],1),0.0);
+            IF (CABS1(S[l-1]) # 0.0) THEN
+              IF (CABS1(X[l-1,l-1]) # 0.0) THEN
+                S[l-1] := CSIGN(S[l-1],X[l-1,l-1]);
+              END;
+              zscal(N-l+1,1.0/S[l-1],X[l-1,l-1],1);
+              X[l-1,l-1] := one + X[l-1,l-1];
+            END; (* IF *)
+            S[l-1] := -S[l-1];
+          END; (* IF *)
+          IF (P >= lp1) THEN
+           FOR j:=lp1 TO P DO
+              IF (l <= nct) THEN
+                IF (CABS1(S[l-1]) # 0.0) THEN
+                  (* apply the transformation. *)
+                  t := -zdotc(N-l+1,X[l-1,l-1],1,X[j-1,l-1],1)/X[l-1,l-1];
+                  zaxpy(N-l+1,t,X[l-1,l-1],1,X[j-1,l-1],1);
+                END; (* IF *)
+              END; (* IF *)
+              (* place the l-th row of x into e for the            *)
+              (* subsequent calculation of the row transformation. *)
+              E[j-1] := conj(X[j-1,l-1]);
+            END; (* FOR *)
+          END; (* IF *)
+          IF NOT ((NOT wantu) OR (l > nct)) THEN
+            (* place the transformation in u for *)
+            (* subsequent back multiplication.   *)
+            FOR i:=l TO N DO
+              U[l-1,i-1] := X[l-1,i-1];
+            END; (* FOR *)
+          END; (* IF *)
+          IF (l <= nrt) THEN
+           (* compute the l-th row transformation and *)
+           (* place the l-th super-diagonal in e(l).  *)
+            E[l-1] := CMPLX(dznrm2(P-l,E[lp1-1],1),0.0);
+            IF (CABS1(E[l-1]) # 0.0) THEN
+              IF (CABS1(E[lp1-1]) # 0.0) THEN
+                E[l-1] := CSIGN(E[l-1],E[lp1-1]);
+              END; (* IF *)
+              zscal(P-l,1.0/E[l-1],E[lp1-1],1);
+              E[lp1-1] := one + E[lp1-1];
+            END; (* IF *)
+            E[l-1] := -conj(E[l-1]);
+            IF (lp1 <= N) AND (CABS1(E[l-1]) # 0.0) THEN
+             (* apply the transformation. *)
+             FOR i:=lp1 TO N DO Work[i-1] := zero; END;
+             FOR j:=lp1 TO P DO
+               zaxpy(N-l,E[j-1],X[j-1,lp1-1],1,Work[lp1-1],1);
+             END;
+             FOR j:=lp1 TO P DO
+               zaxpy(N-l,conj(-E[j-1]/E[lp1-1]),Work[lp1-1],1,X[j-1,lp1-1],1);
+              END; (* FOR *)
+            END; (* IF *)
+            IF wantv THEN
+              (* place the transformation in V for *)
+              (* subsequent back multiplication.   *)
+              FOR i:=lp1 TO P DO
+                V[l-1,i-1] := E[i-1];
+              END;
+            END; (* IF *)
+          END; (* IF *)
+        END; (* FOR *)
+      END; (* IF *)
+
+      (* set up the final bidiagonal matrix or order m *)
+
+      m := MIN0(P,N+1);
+      nctp1 := nct + 1;
+      nrtp1 := nrt + 1;
+      IF (nct < P) THEN
+        S[nctp1-1] := X[nctp1-1,nctp1-1];
+      END; (* IF *)
+      IF (N < m) THEN
+        S[m-1] := zero;
+      END; (* IF *)
+      IF (nrtp1 < m) THEN
+        E[nrtp1-1] := X[m-1,nrtp1-1];
+      END; (* IF *)
+      E[m-1] := zero;
+
+      IF (wantu) THEN (* if required, generate U *)
+        IF (ncu >= nctp1) THEN
+          FOR j:=nctp1 TO ncu DO
+            FOR i:=1 TO N DO U[j-1,i-1] := zero; END;
+            U[j-1,j-1] := one;
+          END; (* FOR *)
+        END; (* IF *)
+        IF (nct >= 1) THEN
+          FOR ll:=1 TO nct DO
+            l := nct - ll + 1;
+            IF (CABS1(S[l-1]) = 0.0) THEN
+              FOR i:=1 TO N DO
+                U[l-1,i-1] := zero;
+              END; (* FOR *)
+              U[l-1,l-1] := one;
+            ELSE
+              lp1 := l + 1;
+              IF (ncu >= lp1) THEN
+               FOR j:=lp1 TO ncu DO
+                  t := -zdotc(N-l+1,U[l-1,l-1],1,U[j-1,l-1],1) /
+                        U[l-1,l-1];
+                  zaxpy(N-l+1,t,U[l-1,l-1],1,U[j-1,l-1],1);
+                END; (* FOR *)
+              END; (* IF *)
+              zscal(N-l+1,-one,U[l-1,l-1],1);
+              U[l-1,l-1] := one + U[l-1,l-1];
+              lm1 := l - 1;
+              IF (lm1 >= 1) THEN
+                FOR i := 1 TO lm1 DO U[l-1,i-1] := zero; END;
+              END; (* IF *)
+            END; (* IF *)
+          END; (* FOR *)
+        END; (* IF *)
+      END; (* IF *)
+
+      IF (wantv) THEN (* if it is required, generate V *)
+        FOR ll:=1 TO P DO
+          l := P - ll + 1;
+          lp1 := l + 1;
+          IF (l <= nrt) THEN
+            IF (CABS1(E[l-1]) # 0.0) THEN
+              FOR j:=lp1 TO P DO
+                t := -zdotc(P-l,V[l-1,lp1-1],1,V[j-1,lp1-1],1)/V[l-1,lp1-1];
+                zaxpy(P-l,t,V[l-1,lp1-1],1,V[j-1,lp1-1],1);
+              END;
+            END; (* IF *)
+          END; (* IF *)
+          FOR i:=0 TO P-1 DO V[l-1,i] := zero; END;
+          V[l-1,l-1] := one;
+        END; (* FOR *)
+      END; (* IF *)
+
+
+      (* transform s and e so that they are double precision. *)
+
+      FOR i:=1 TO m DO
+        IF (CABS1(S[i-1]) # 0.0) THEN
+          t := CMPLX(CABS(S[i-1]),0.0);
+          r := S[i-1] / t;
+          S[i-1] := t;
+          IF (i < m) THEN
+            E[i-1] := E[i-1] / r;
+          END;
+          IF (wantu) THEN
+            zscal(N,r,U[i-1,1 -1],1);
+          END; (* IF *)
+        END; (* IF *)
+
+        IF (i # m) THEN
+          IF (CABS1(E[i-1]) # 0.0) THEN
+            t := CMPLX(CABS(E[i-1]),0.0);
+            r := t / E[i-1];
+            E[i-1] := t;
+            S[i+1-1] := S[i+1-1]*r;
+            IF (wantv) THEN
+              zscal(P,r,V[i+1-1,1-1],1);
+            END;
+          END; (* IF *)
+        END; (* IF *)
+      END; (* FOR *)
+
+      (* main iteration loop for the singular values. *)
+
+      mm := m;
+      iter := 0;
+
+      (* quit if all the singular values have been found. *)
+
+      (* exit ??? *)
+
+      WHILE (m # 0) DO
+
+        (* if too many iterations have been performed, set *)
+        (* flag and return. *)
+
+        IF (iter < maxit) THEN
+
+          (* 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(m) and e(l-1) are negligible and l.lt.m    *)
+          (* kase = 2 if s(l) is negligible and l.lt.m                *)
+          (* kase = 3 if e(l-1) is negligible, l.lt.m, and            *)
+          (* s(l), ..., s(m) are not negligible (qr step).            *)
+          (* kase = 4 if e(m-1) is negligible (convergence).          *)
+
+          ll:=1;
+          LOOP
+            l := m - ll;
+            IF (l = 0) THEN EXIT; END;
+            test := CABS(S[l-1]) + CABS(S[l+1-1]);
+            ztest := test + CABS(E[l-1]);
+            IF (ztest = test) THEN
+              E[l-1] := zero;
+              EXIT;
+            END;
+            INC(ll);
+          END;
+
+          IF (l # m-1) THEN
+            lp1 := l + 1;
+            mp1 := m + 1;
+
+            ls := l; (* Wg. Compilerwarnung *)
+            lls:=lp1;
+            LOOP
+              IF (lls > mp1) THEN EXIT; END;
+              ls := m - lls + lp1;
+              IF (ls = l) THEN
+                EXIT;
+              END;
+              test := 0.0;
+              IF (ls # m) THEN
+                test:= CABS(E[ls-1]); (* + test *)
+              END;
+              IF (ls # l+1) THEN
+                test := test + CABS(E[ls-1-1]);
+              END;
+              ztest := test + CABS(S[ls-1]);
+              IF (ztest = test) THEN
+                S[ls-1-1] := zero;
+                EXIT;
+              END;
+              INC(lls);
+            END;
+
+            IF (ls = l) THEN
+              kase := 3;
+            ELSIF (ls # m) THEN
+              kase := 2;
+              l := ls;
+            ELSE
+              kase := 1;
+            END;
+          ELSE
+            kase := 4;
+          END; (* IF *)
+          INC(l);
+
+          (* perform the task indicated by kase. *)
+
+          IF (kase = 1) THEN (* deflate negligible s(m). *)
+
+            mm1 := m - 1;
+            f := RE(E[m-1-1]);
+            E[m-1-1] := zero;
+            FOR kk:=l TO mm1 DO
+              k := mm1 - kk + l;
+              t1 := RE(S[k-1]);
+              drotg(t1,f,cs,sn);
+              S[k-1] := CMPLX(t1,0.0);
+              IF (k # l) THEN
+                f := -sn*RE(E[k-1-1]);
+                E[k-1-1] := scalarMult(cs,E[k-1-1]);
+              END;
+              IF (wantv) THEN
+                zdrot(P,V[k-1,1 -1],1,V[m-1,1 -1],1,cs,sn);
+              END; (* IF *)
+            END; (* FOR *)
+
+          ELSIF (kase = 2) THEN (* split at negligible s(l). *)
+
+            f := RE(E[l-1-1]);
+            E[l-1-1] := zero;
+            FOR k:=l TO m DO
+              t1 := RE(S[k-1]);
+              drotg(t1,f,cs,sn);
+              S[k-1] := CMPLX(t1,0.0);
+              f := -sn*RE(E[k-1]);
+              E[k-1] := scalarMult(cs,E[k-1]);
+              IF (wantu) THEN
+                zdrot(N,U[k-1,1 -1],1,U[l-1 -1,1 -1],1,cs,sn);
+              END;
+            END; (* FOR *)
+
+          ELSIF (kase = 3) THEN (* perform one qr step *)
+
+            (* calculate the shift. *)
+            scale := DMAX1(DMAX1(CABS(S[m -1]),CABS(S[m-1-1])),
+                           DMAX1(CABS(E[m-1-1]),CABS(S[l -1])));
+            scale := DMAX1(scale,CABS(E[l-1]));
+            sm := RE(S[m -1]) / scale;
+            smm1 := RE(S[m-1-1]) / scale;
+            emm1 := RE(E[m-1-1]) / scale;
+            sl := RE(S[l -1]) / scale;
+            el := RE(E[l -1]) / scale;
+            b := ((smm1 + sm)*(smm1 - sm) + emm1*emm1) / 2.0;
+            c := sqr(sm*emm1);
+            shift := 0.0;
+            IF (b # 0.0) OR (c # 0.0) THEN
+              shift := sqrt(b*b + c);
+              IF (b < 0.0) THEN
+                shift := -shift;
+              END;
+              shift := c/(b+shift);
+            END; (* IF *)
+            f := (sl+sm)*(sl-sm) + shift;
+            g := sl*el;
+
+            (* chase zeros. *)
+
+            mm1 := m - 1;
+            FOR k:=l TO mm1 DO
+              drotg(f,g,cs,sn);
+              IF (k # l) THEN
+                E[k-1-1] := CMPLX(f,0.0);
+              END;
+              f := cs*RE(S[k-1]) + sn*RE(E[k-1]);
+              E[k-1] := scalarMult(cs,E[k-1]) - scalarMult(sn,S[k-1]);
+              g := sn*RE(S[k+1-1]);
+              S[k+1-1] := scalarMult(cs,S[k+1-1]);
+              IF wantv THEN
+                zdrot(P,V[k-1,1 -1],1,V[k+1 -1,1 -1],1,cs,sn);
+              END;
+              drotg(f,g,cs,sn);
+              S[k-1] := CMPLX(f,0.0);
+              f := cs*RE(E[k-1]) + sn*RE(S[k+1-1]);
+              S[k+1-1] := -scalarMult(sn,E[k-1]) + scalarMult(cs,S[k+1-1]);
+              g := sn*RE(E[k+1-1]);
+              E[k+1-1] := scalarMult(cs,E[k+1-1]);
+              IF wantu AND (k < N) THEN
+                zdrot(N,U[k-1,1 -1],1,U[k+1 -1,1 -1],1,cs,sn);
+              END;
+            END; (* FOR *)
+            E[m-1-1] := CMPLX(f,0.0);
+            INC(iter);
+
+          ELSIF (kase = 4) THEN
+
+            (* convergence, make the singular value positive *)
+
+            IF (RE(S[l-1]) < 0.0) THEN
+              S[l-1] := -S[l-1];
+              IF (wantv) THEN
+                zscal(P,CMPLX(-1.0,0.0),V[l-1,1 -1],1)
+              END;
+            END;
+            (* order the singular value. *)
+            WHILE (l # mm) AND (RE(S[l-1]) < RE(S[l+1-1])) DO
+              t := S[l-1];
+              S[l-1] := S[l+1-1];
+              S[l+1-1] := t;
+              IF wantv AND (l < P) THEN
+                zswap(P,V[l-1,1 -1],1,V[l+1 -1,1 -1],1)
+              END; (* IF *)
+              IF wantu AND (l < N) THEN
+                zswap(N,U[l-1,1 -1],1,U[l+1 -1,1 -1],1)
+              END; (* IF *)
+              INC(l);
+            END; (* WHILE *)
+            iter := 0;
+            DEC(m);
+
+          END; (* IF kase *)
+
+        ELSE
+          Info := m;
+          RETURN;
+        END; (* IF *)
+      END; (* FOR *)
+END zSVDc;
+
+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);
+
+          PROCEDURE SQ(c : LONGCOMPLEX) : LONGREAL;
+
+          BEGIN
+                 RETURN RE(c*conj(c));
+          END SQ;
+
+          CONST eps  = 1.0E+01*MachEps; 
+                MaxI = MAX(INTEGER);
+
+          VAR   p,q,j          : INTEGER;
+                red,off,thresh : LONGREAL;
+                sqp,sqq,t,invc : LONGREAL;
+                f,x,y          : LONGCOMPLEX;
+                ev1,ev2        : POINTER TO ARRAY [0..MaxI-1] OF LONGCOMPLEX;
+                sweep          : INTEGER;
+BEGIN
+      ALLOCATE(ev1,N*SIZELONGCMPLX);
+      ALLOCATE(ev2,N*SIZELONGCMPLX);
+ 
+      IF (ev1 = NIL) OR (ev2 = NIL) THEN
+        Errors.ErrOut("Kein Freispeicher vorhanden (Tagati)");
+        IF (ev1 # NIL) THEN DEALLOCATE(ev1,N*SIZELONGCMPLX); END;
+        RETURN;
+      END;
+ 
+      FOR p:=0 TO N-1 DO
+        ev1^[p] := zero;
+        ev2^[p] := A[p,p];
+      END;
+ 
+      FOR p:=0 TO N-1 DO
+        FOR q:=0 TO N-1 DO
+          U[p,q] := zero;
+        END;
+        U[p,p] := one;
+      END;
+         
+      red := 0.04 / CardPot(VAL(LONGREAL,N),4);
+ 
+      sweep:=0;
+      LOOP
+        INC(sweep);
+ 
+        IF (sweep > 50) THEN
+          Errors.ErrOut("Bad convergence in TakagiFactor");
+          EXIT;
+        END;
+ 
+        off := 0.0;
+        FOR q:=1 TO N-1 DO
+          FOR p:=0 TO q-1 DO
+            off:=off + SQ(A[q,p]);
+          END;
+        END;
+        IF (off <=  eps*eps) THEN
+          EXIT;
+        END;
+        thresh := 0.0;
+        IF (sweep < 4) THEN
+          thresh := off*red
+        END;
+ 
+        FOR q:=1 TO N-1 DO
+          FOR p:=0 TO q-1 DO
+            off := SQ(A[q,p]);
+            sqp := SQ(ev2^[p]);
+            sqq := SQ(ev2^[q]);
+            IF (sweep > 4) AND (off < eps*(sqp+sqq)) THEN
+              A[q,p] := zero;
+            ELSIF (off > thresh) THEN
+              t := 0.5*ABS(sqp - sqq);
+              IF (t > 0.0) THEN
+                f := scalarMult(sign2(1.0,sqp-sqq),(ev2^[q]*conj(A[q,p]) + 
+                                                    conj(ev2^[p])*A[q,p]));
+              ELSE
+                f := one;
+                IF (sqp # 0) THEN
+                  f := LongComplexMath.sqrt(ev2^[q] / ev2^[p])
+                END;
+              END;
+              t:=t + sqrt(t*t + SQ(f));
+              f:=f / CMPLX(t,0.0);
+ 
+              ev1^[p] := ev1^[p] + A[q,p]*conj(f);
+              ev2^[p] := A[p,p]  + ev1^[p];
+              ev1^[q] := ev1^[q] - A[q,p]*f;
+              ev2^[q] := A[q,q]  + ev1^[q];
+ 
+              t := SQ(f);
+              invc := sqrt(t +1.0);
+              f:=f / CMPLX(invc,0.0);
+              t:=t / (invc*(invc + 1.0));
+ 
+              FOR j:=0 TO p-1 DO
+                x := A[p,j];
+                y := A[q,j];
+                A[p,j] := x + (conj(f)*y - CMPLX(t,0.0)*x);
+                A[q,j] := y -      (f *x + CMPLX(t,0.0)*y);
+              END;
+ 
+              FOR j:=p+1 TO q-1 DO
+                x := A[j,p];
+                y := A[q,j];
+                A[j,p] := x + (conj(f)*y - CMPLX(t,0.0)*x);
+                A[q,j] := y -      (f *x + CMPLX(t,0.0)*y);
+              END;
+ 
+              FOR j:=q+1 TO N-1 DO
+                x := A[j,p];
+                y := A[j,q];
+                A[j,p] := x + (conj(f)*y - CMPLX(t,0.0)*x);
+                A[j,q] := y -      (f *x + CMPLX(t,0.0)*y);
+              END;
+ 
+              A[q,p] := zero;
+ 
+              FOR j:=0 TO N-1 DO
+                x := U[j,p];
+                y := U[j,q];
+                U[j,p] := x +      (f *y - CMPLX(t,0.0)*x);
+                U[j,q] := y - (conj(f)*x + CMPLX(t,0.0)*y);
+              END;
+ 
+            END; (* IF *)
+          END;
+        END;
+ 
+        FOR p:=0 TO N-1 DO
+          ev1^[p] := zero;
+          A[p,p]  := ev2^[p];
+        END;
+ 
+      END; (* LOOP *)
+ 
+      (* make the diagonal elements nonnegative *)
+ 
+      FOR p:=0 TO N-1 DO
+        D[p] := CABS(A[p,p]);
+        IF (D[p] > eps) AND (D[p] # RE(A[p,p])) THEN
+          f := LongComplexMath.sqrt(A[p,p] / CMPLX(D[p],0.0));
+          FOR q:=0 TO N-1 DO
+            U[q,p] := U[q,p]*f;
+          END;
+        END; (* IF *)
+      END;
+ 
+      IF (sort # 0) THEN (* sort the eigenvalues *)
+        FOR p:=0 TO N-2 DO
+          j := p;
+          t := D[p];
+          FOR q := p+1 TO N-1 DO
+            IF (VAL(LONGREAL,sort)*(t - D[q]) > 0.0) THEN
+              j := q;
+              t := D[q];
+            END; (* IF *)
+          END;
+          IF (j # p) THEN
+            D[j] := D[p];
+            D[p] := t;
+            FOR q:=0 TO N-1 DO
+              x      := U[q,p];
+              U[q,p] := U[q,j];
+              U[q,j] := x;
+            END;
+          END; (* IF *)
+        END;
+      END; (* IF sort *)
+ 
+      DEALLOCATE(ev2,N*SIZELONGCMPLX);
+      DEALLOCATE(ev1,N*SIZELONGCMPLX);
+END Takagi;
+
+END SVDLib3.