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.