SVDLib2.mod.m2pp
989 lines (926 with data), 32.9 kB
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 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 (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: *)
(* *)
(* - Weiteres testen von PdSVDc *)
(* - Parameter m,n und Form der Speicherung in A durchsehen, einheitlich *)
(* machen und dokumentieren !!! *)
(*------------------------------------------------------------------------*)
(* Tests for PdSVD & PdSVDc are given in TstSVD1u2a *)
(*------------------------------------------------------------------------*)
(* 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;
BEGIN
IF (y >= 0.0) THEN RETURN ABS(x); ELSE RETURN -ABS(x); END;
END sign;
PROCEDURE max(x,y: LONGREAL): LONGREAL;
BEGIN
IF (x > y) THEN RETURN x; ELSE RETURN y; END;
END max;
CONST MaxIter = 30;
precise = TRUE;
VAR nm,l,k,j,jj,its,i,mnmin : INTEGER;
z,y,x,scale,s,h,g,f,c,anorm : LONGREAL;
Split,Konv : BOOLEAN;
rv1,tmpvek : PVEKTOR;
BEGIN
ALLOCATE(rv1,n*TSIZE(LONGREAL));
IF precise THEN ALLOCATE(tmpvek,m*TSIZE(LONGREAL)); END;
g:=0.0;
scale:=0.0;
anorm:=0.0;
FOR i:=1 TO n DO (* Housholder-Transformation *)
l := i + 1;
rv1^[i] := scale*g;
g:=0.0;
s:=0.0;
IF (i <= m) THEN
IF precise THEN
FOR k:=i TO m DO tmpvek^[k] := A^[k]^[i]; END;
scale := AbsSumVek(tmpvek^,i,m);
ELSE
scale := 0.0;
FOR k:=i TO m DO scale:=scale + ABS(A^[k]^[i]); END;
END;
IF (scale # 0.0) THEN
FOR k:=i TO m DO
A^[k]^[i] := A^[k]^[i] / scale;
s:=s + A^[k]^[i]*A^[k]^[i];
END;
f := A^[i]^[i];
g := -sign(sqrt(s),f);
h := f*g-s;
A^[i]^[i] := f - g;
FOR j:=l TO n DO
s:=0.0; FOR k:=i TO m DO s:=s + A^[k]^[i]*A^[k]^[j]; END;
f := s / h;
FOR k:=i TO m DO A^[k]^[j]:=A^[k]^[j] + f*A^[k]^[i]; END;
END;
FOR k:=i TO m DO A^[k]^[i] := scale*A^[k]^[i]; END;
END
END; (* IF i *)
W[i-1] := scale*g;
g:=0.0;
s:=0.0;
scale := 0.0;
IF (i <= m) AND (i # n) THEN
FOR k:=l TO n DO scale:=scale + ABS(A^[i]^[k]); END;
IF (scale # 0.0) THEN
FOR k:=l TO n DO
A^[i]^[k] := A^[i]^[k] / scale;
s:=s + A^[i]^[k]*A^[i]^[k];
END;
f := A^[i]^[l];
g := -sign(sqrt(s),f);
h := f*g - s;
A^[i]^[l] := f-g;
FOR k:=l TO n DO rv1^[k] := A^[i]^[k] / h; END;
FOR j:=l TO m DO
s:=0.0;
FOR k:=l TO n DO s:=s + A^[j]^[k]*A^[i]^[k]; END;
FOR k:=l TO n DO A^[j]^[k]:=A^[j]^[k] + s*rv1^[k]; END;
END;
FOR k:=l TO n DO A^[i]^[k] := scale*A^[i]^[k]; END;
END
END;
anorm := max(anorm,(ABS(W[i-1])+ABS(rv1^[i])));
END; (* FOR i *)
(* Akkumulation der Housholdertransformationen *)
(* Hier muesste l=n+1 sein ... *)
V^[n]^[n] := 1.0; g := rv1^[n]; l := n;
FOR i:=n-1 TO 1 BY -1 DO (* Rechtsseitige Transformation *)
IF (g # 0.0) THEN
FOR j:=l TO n DO V^[i]^[j] := (A^[i]^[j] / A^[i]^[l]) / g; END;
FOR j:=l TO n DO
s:=0.0;
FOR k:=l TO n DO s:=s + A^[i]^[k]*V^[j]^[k]; END;
FOR k:=l TO n DO V^[j]^[k]:=V^[j]^[k] + s*V^[i]^[k]; END;
END;
END;
FOR j:=l TO n DO V^[j]^[i]:=0.0; V^[i]^[j]:=0.0; END;
V^[i]^[i] := 1.0;
g := rv1^[i];
l := i;
END; (* FOR i *)
(*
l:=n+1; (* Wg. Compilerwarung *)
FOR i:=n TO 1 BY -1 DO (* Rechtsseitige Transformation *)
IF (i < n) THEN
IF (g # 0.0) THEN
FOR j:=l TO n DO V^[i]^[j] := (A^[i]^[j] / A^[i]^[l]) / g; END;
FOR j:=l TO n DO
s:=0.0;
FOR k:=l TO n DO s:=s + A^[i]^[k]*V^[j]^[k]; END;
FOR k:=l TO n DO V^[j]^[k]:=V^[j]^[k] + s*V^[i]^[k]; END;
END;
END;
FOR j:=l TO n DO V^[j]^[i]:=0.0; V^[i]^[j]:=0.0; END;
END;
V^[i]^[i] := 1.0;
g := rv1^[i];
l := i;
END; (* FOR i *)
*)
IF (m < n) THEN mnmin:=m; ELSE mnmin:=n; END;
FOR i:=mnmin TO 1 BY -1 DO (* Linkssseitige Transformation *)
l := i+1;
g := W[i-1];
FOR j:=l TO n DO A^[i]^[j] := 0.0; END;
IF (g # 0.0) THEN
g := 1.0 / g;
FOR j:=l TO n DO
s := 0.0;
FOR k:=l TO m DO s:=s + A^[k]^[i]*A^[k]^[j]; END;
f := (s/A^[i]^[i])*g;
FOR k:=i TO m DO A^[k]^[j]:=A^[k]^[j] + f*A^[k]^[i]; END;
END;
FOR j:=i TO m DO A^[j]^[i] := A^[j]^[i]*g; END;
ELSE
FOR j:=i TO m DO A^[j]^[i] := 0.0; END;
END;
A^[i]^[i]:=A^[i]^[i] + 1.0
END; (* FOR i *)
FOR k:=n TO 1 BY -1 DO (* Diagonalisierung der Bidiagonalmatrix *)
Konv:=FALSE; its:=0;
REPEAT
INC(its);
l:=k; Split:=FALSE;
LOOP (* FOR l:=k TO 1 BY -1 DO *) (* Teste auf Blockungen *)
nm := l-1;
IF ((ABS(rv1^[l]) + anorm) = anorm) THEN Split:=TRUE; EXIT END;
IF (nm > 0) THEN
IF ((ABS(W[nm-1]) + anorm) = anorm) THEN EXIT END;
END;
IF (l = 1) THEN EXIT END;
DEC(l);
END; (* LOOP *)
IF NOT Split THEN
c := 0.0; s := 1.0;
i:=l;
WHILE (i <= k) AND NOT Split DO (* FOR i:=l TO k DO *)
f := s*rv1^[i];
rv1^[i] := c*rv1^[i];
IF ((ABS(f) + anorm) = anorm) THEN
Split:=TRUE; (* !!! *)
ELSE
g := W[i-1];
h := Pythag(f,g);
W[i-1] := h;
h := 1.0 / h;
c := (g*h);
s := -(f*h);
FOR j:=1 TO m DO
y := A^[j]^[nm];
z := A^[j]^[i];
A^[j]^[nm] := (y*c) + (z*s);
A^[j]^[i] := -(y*s) + (z*c);
END;
END; (* IF *)
INC(i);
END; (* WHILE *)
END; (* IF NOT Split *)
z := W[k-1];
IF (l = k) THEN (* Konvergenz ! *)
Konv:= TRUE;
IF (z < 0.0) THEN (* Mache den Singul"arwert positiv. *)
W[k-1] := -z;
FOR j:=1 TO n DO V^[k]^[j] := -V^[k]^[j]; END;
END;
ELSE (* l # k *)
IF (its = 30) THEN
Fehler:=TRUE;
Fehlerflag:= 'MaxIter ueberschritten (SVD)';
ErrOut(Fehlerflag); RETURN;
END;
x := W[l-1];
nm := k-1;
y := W[nm-1];
g := rv1^[nm];
h := rv1^[k];
f := ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
g := Pythag(f,1.0);
f := ((x-z)*(x+z) + h*((y / (f + sign(g,f))) - h)) / x;
c := 1.0; s:=1.0;
FOR j:=l TO nm DO
i := j+1;
g := rv1^[i];
y := W[i-1];
h := s*g;
g := c*g;
z := Pythag(f,h);
rv1^[j] := z;
c := f / z;
s := h / z;
f := (x*c) + (g*s);
g := -(x*s) + (g*c);
h := y*s;
y := y*c;
FOR jj:=1 TO n DO
x := V^[j]^[jj];
z := V^[i]^[jj];
V^[j]^[jj] := (x*c) + (z*s);
V^[i]^[jj] := -(x*s) + (z*c)
END;
z := Pythag(f,h);
W[j-1] := z;
IF (z # 0.0) THEN
z := 1.0 / z;
c := f*z;
s := h*z;
END;
f := (c*g) + (s*y);
x := -(s*g) + (c*y);
FOR jj:=1 TO m DO
y := A^[jj]^[j];
z := A^[jj]^[i];
A^[jj]^[j] := (y*c) + (z*s);
A^[jj]^[i] := -(y*s) + (z*c);
END;
END;
rv1^[l] := 0.0;
rv1^[k] := f;
W[k-1] := x;
END; (* IF Konvergenz ! *)
UNTIL Konv OR (its > MaxIter);
END; (* FOR k *)
IF precise THEN DEALLOCATE(tmpvek,m*TSIZE(LONGREAL)); END;
DEALLOCATE(rv1,n*TSIZE(LONGREAL));
(*
FOR i:=1 TO n DO
FOR j:=1 TO n DO
x:=V^[i]^[j]; V^[i]^[j]:=V^[j]^[i]; V^[j]^[i]:=x;
END;
END;
*)
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;
Zw : POINTER TO ARRAY [1..8191] OF LONGREAL; (* n *)
Svek : POINTER TO ARRAY [1..8191] OF LONGREAL; (* MAX(n,m) *)
BEGIN
ALLOCATE(Zw,n*TSIZE(LONGREAL));
ALLOCATE(Svek,MaxCard(n,m)*TSIZE(LONGREAL));
IF (Zw = NIL) OR (Svek = NIL) THEN
FatalError("Kein Freispeicher vorhanden (SVDLib1.SVDSolv) !");
END;
FOR j:=1 TO n DO
s := 0.0;
IF (W[j-1] # 0.0) THEN
IF (CAP(Utr) # "T") THEN
FOR i:=1 TO m DO Svek^[i] := U^[i]^[j]*C[i-1]; END;
ELSE
FOR i:=1 TO m DO Svek^[i] := U^[j]^[i]*C[i-1]; END;
END;
s := SumVek(Svek^,1,m);
s := s / W[j-1];
END;
Zw^[j] := s;
END;
FOR j:=1 TO n DO
FOR jj:=1 TO n DO Svek^[jj] := V^[jj]^[j]*Zw^[jj]; END;
s := SumVek(Svek^,1,n);
X[j-1] := s;
END;
DEALLOCATE(Svek,MaxCard(n,m)*TSIZE(LONGREAL));
DEALLOCATE(Zw,n*TSIZE(LONGREAL));
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;
BEGIN
FOR i:=1 TO n DO
MaxW := W[i-1]; MaxI := i;
FOR j:=i+1 TO n DO
IF (W[j-1] > MaxW) THEN MaxW:=W[j-1]; MaxI:=j; END;
END;
IF (i # MaxI) THEN
zw := W[i-1]; W[i-1]:=MaxW; W[MaxI-1]:= zw;
FOR j:=1 TO n DO
zw := V^[i]^[j];
V^[i]^[j] := V^[MaxI]^[j];
V^[MaxI]^[j] := zw;
END;
FOR j:=1 TO m DO
zw := U^[j]^[i];
U^[j]^[i] := U^[j]^[MaxI];
U^[j]^[MaxI] := zw;
END;
END; (* IF *)
END;
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.