SVDLib1.mod.m2pp
2350 lines (2185 with data), 78.7 kB
IMPLEMENTATION MODULE SVDLib1;
(*========================================================================*)
(* WICHTIG: BITTE NUR DIE DATEI SVDLib1.mod.m2pp EDITIEREN !!! *)
(*========================================================================*)
(* This file contains three version - one calling Modula-2 LibDBlas *)
(* the other the Fortran dblas level 1 routines (LibDBlasF77) *)
(* and the third inlines the code *)
(* *)
(* m2pp -D __{OPTION}__ < SVDLib1.mod.m2pp > SVDLib1.mod *)
(* or *)
(* m2pp -b -D __{OPTION}__ < SVDLib1.mod.m2pp > SVDLib1.mod *)
(* *)
(* to preserve line information from *.m2pp file where OPTION may be one *)
(* of {BLAS|BLASF77|INLINE} *)
(* *)
(* In addition an option debug can be activated which gives out a lot of *)
(* intermediate results - not for normal use ;-) *)
(*------------------------------------------------------------------------*)
(* Anpassung an Modula-2 : Michael Riedl *)
(* *)
(* Letzte Bearbeitung: *)
(* *)
(* 31.03.16, MRi: Erstellend der ersten Rohversion von dSVDc *)
(* 01.04.16, MRi: Erstellen der ersten Version von dSVD *)
(* 02.04.16, MRi: Kleine Korrekturen, umstellen der Indizes für Matrix V *)
(* in dSVD *)
(* 03.04.16, MRi: Erste uebersetzbare Version von dSVDc *)
(* 03.05.16, MRi: Ermittelt S und V & U korrekt (dSVDc) *)
(* 05.05.16, MRi: Erste Version mit offenen Feldern (dSVDc) *)
(* 12.07.16, MRi: Erstellen der ersten Rohversion von MinFit - mit Pascal *)
(* uebersetzbar. *)
(* 15.07.16, MRi: Anpassen der Kontrollstrukturen - elimination von GOTO *)
(* nach Vorlage von MinFit (Praxis Unterroutine) *)
(* 25.07.16, MRi: Abfragen fuer Cancellation/TestFconvergence in MinFit *)
(* ueberarbeitet - XDS M2 macht hier selsame Dinge ... *)
(* 28.07.16, MRi: Erstellen der ersten Version von GivSVD *)
(* 01.08.16, MRi: Erstellen der ersten Rohversion von JacobiSVD *)
(* 02.08.16, MRi: Erste uebersetzbare Version von JacobiSVD *)
(* Kleine Korrekturen in JacobiSVD - erste Testergebnisse *)
(* sehen gut aus *)
(* 03.08.16, MRi: Erstellen der ersten Version von PowSVD, erste Test- *)
(* ergebnisse sind zufriedenstellend *)
(* 08.08.16, MRi: Einfuegen von CalcResiduals *)
(* Korrekturen an GivSVD, Ergebnisse sehen gut aus *)
(* Anpassen von GivSVD an offene Felder *)
(* 21.08.16, MRi: Ersetzen von ErrOut durch setzen von Fehlerflag *)
(* 22.08.16, MRi: Kleine Korrekturen an JacobiSVD *)
(* Umstellen der Schleifenstruktur in CalcLQsol so dass *)
(* die Reorganisation von JacobiSVD nachgezogen wird. *)
(* Zudem den Parameter "rank" eingefuegt um die Rang- *)
(* reduktion durch JacobiSVD oder PowSVD nutzen zu koennen *)
(* Die alte Ordnung (z.B. von dSVD) kann ueber den Para- *)
(* meter IfT#"{Y|y]" angewaehlt werden. *)
(* 22.08.16, MRi: Kleine Korrekturen bei der Durchfuehrung der Rotation *)
(* in JacobiSVD (r->h, g neu) *)
(* 04.12.17, MRi: "verfeinerte" Rotationsdurchfuehrung in JacobiSVD *)
(* 07.12.17, MRi: Einfuegen von SVDSolv *)
(* 01.07.18, MRi: Zusammenfuehren aller Routinen in einem Modul *)
(* 06.07.18, MRi: Korrektur in CalcLQres1 *)
(* 07.07.18, MRi: Erstellen der ersten Rohversion von MinFitV *)
(* 08.07.18, MRi: Alle GOTO in MinFitV eliminiert *)
(* Parameter A und D auf offene Felder umgestellt *)
(* 09.07.18, MRi: Level1 BLAS routinen eingefuegt *)
(* 16.08.18, MRi: Ersetzen MinFit durch neu ubersetzte Version (MinFitV) *)
(* 18.08.18, MRi: Umstellen der Indizierung von B in GivSVD um kompatibel *)
(* mit MinFit zu sein. *)
(* 21.08.18, MRi: F77func Aufrufe durch vorhandene lokale Prozeduren *)
(* ersetzt *)
(* 24.08.18, MRi: Die rechte Seite der Gleichungssysteme auf einheit- *)
(* liches Speicherformat [1..nB,1..M] gebracht *)
(* 25.08.18, MRi: In BerLsg die Schleifenstrukturen optimiert *)
(* 26.08.18, MRi: In BerLsg und MinFit komplett auf offene Felder umge- *)
(* stellt *)
(* 28.08.18, MRi: In MinFit die Indexberechnung bereinigt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Durchsehen der Fehlermeldungen *)
(* - Durchsehen der Kommentare und Erklaerungen *)
(* - Ermitteln besserer Parameter für tol & delta *)
(* - Austesten von dSVDc *)
(* - Korrektur der Indizes (momentan nicht optimal) *)
(* - Fehlerhafte Ausgabe fuer "Test Nr 0 - 5x4" in dSVD ermitteln - der *)
(* Fehler tritt nur in den Versionen mit BLAS oder BLASF77 auf. *)
(* Auch dSVDc ist nicht betroffen (WICHTIG !!!) *)
(* - Parameter m,n und Form der Speicherung in A durchsehen, einheitlich *)
(* machen und dokumentieren !!! *)
(*------------------------------------------------------------------------*)
(* Code based in parts on *)
(* *)
(* Nash, J. C.; "Compact numerical methods for computers: linear algebra *)
(* and function minimisation, Adam Hilger, Bristol (UK) (1990) *)
(* *)
(* Code put the LGPL with written authorisation of the author *)
(*------------------------------------------------------------------------*)
(* Test routine for SVD routines are given in TstSVDLib1a.mod *)
(* Test routine for MinFit, GivSVD and oterh LQ related routines are *)
(* given in TstSVD1u2LQ.mod *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SVDLib1.mod.m2pp,v 1.5 2018/08/28 12:30:50 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
IMPORT Errors;
FROM Errors IMPORT Fehlerflag;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps,Small,MachEpsR4,Pythag,MaxCard;
FROM F77func IMPORT MIN0,MAX0,DMAX,DSIGN;
FROM RandomLib IMPORT Zufall;
FROM MatLib IMPORT SumVek,NeumaierProdSum;
<* IF (__debug__) THEN *>
IMPORT FMatEA,FileSystem;
IMPORT TIO;
<* END *>
<* IF (__INLINE__) THEN *>
FROM LibDBlas IMPORT drotg;
(*******
<* IF NOT (__debug__) THEN *>
IMPORT TIO;
<* END *>
*******)
<* END *>
<* IF (__BLAS__) THEN *>
FROM LibDBlas IMPORT daxpy,ddot,dnrm2,drot,drotg,dscal,dswap;
(*******
<* IF NOT (__debug__) THEN *>
IMPORT TIO;
<* END *>
*******)
<* END *>
<* IF (__BLASF77__) THEN *>
FROM LibDBlasF77 IMPORT daxpy,ddot,dnrm2,drot,drotg,dscal,dswap;
(*******
<* IF NOT (__debug__) THEN *>
IMPORT TIO;
<* END *>
*******)
<* END *>
<* IF (__BLAS__) THEN *>
CONST tag = "SVDLib1 Modula-2 BLAS";
<* END *>
<* IF (__BLASF77__) THEN *>
CONST tag = "SVDLib1 Fortran BLAS";
<* END *>
<* IF (__INLINE__) THEN *>
CONST tag = "SVDLib1 inlined code";
<* END *>
(*=========================== Interne Routinen =============================*)
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
<* IF (__debug__) THEN *>
PROCEDURE WrVec( str : ARRAY OF CHAR;
VAR X : ARRAY OF LONGREAL;
n : INTEGER);
VAR i : INTEGER;
BEGIN
TIO.WrLn; TIO.WrStr(str); TIO.WrLn;
FOR i:=0 TO n-1 DO TIO.WrLngReal(X[i],24,14); TIO.WrLn; END;
TIO.WrLn;
END WrVec;
<* 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;
(*=========================== Nash Routinen ================================*)
PROCEDURE PowSVD( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR sv : ARRAY OF LONGREAL;
VAR U,V : ARRAY OF ARRAY OF LONGREAL;
klimit : INTEGER;
VAR rank : CARDINAL;
VAR iFehl : INTEGER);
CONST eps = MachEps;
delta = MachEpsR4*MachEpsR4; (* 0.01*MachEpsR4; *)
tol = 10.0*MachEps;
MAXINT = MAX(INTEGER);
VAR i,j,k,nm : INTEGER;
counter,climit : INTEGER;
s,norm,d : LONGREAL;
dPRi,Pj,Ri : LONGREAL;
P,Q,R,W : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
rankdef : BOOLEAN;
BEGIN
nm:=M; IF (N > M) THEN nm:=N; END;
ALLOCATE(P,N*TSIZE(LONGREAL));
ALLOCATE(Q,M*TSIZE(LONGREAL));
ALLOCATE(R,N*TSIZE(LONGREAL));
ALLOCATE(W,N*TSIZE(LONGREAL));
iFehl:=0;
rankdef:=FALSE;
climit := nm*(nm+1)*(nm+2); IF (climit < 64) THEN climit := 64; END;
sv[0]:=0.0;
s:=0.0;
FOR i:=0 TO N-1 DO
Ri:=Zufall();
IF (Ri = 0.0) THEN Ri:=100.0*eps; END;
s:=s + Ri*Ri;
R^[i]:=Ri;
END;
s:=sqrt(s);
FOR i:=0 TO N-1 DO
P^[i] := R^[i] / s;
END;
k:=0;
REPEAT (* step 2 *)
counter:=0;
REPEAT (* step 3 *)
INC(counter);
(* step 4 - compute Q = A x P *)
FOR i:=0 TO M-1 DO Q^[i]:=0.0; END;
FOR j:=0 TO N-1 DO
Pj:=P^[j];
FOR i:=0 TO M-1 DO
Q^[i]:=Q^[i] + A[j,i]*Pj;
END;
END;
norm:=0.0;
FOR i:=0 TO M-1 DO norm:=norm + sqr(Q^[i]); END;
norm:= 1.0 / sqrt(norm);
FOR i:=0 TO M-1 DO Q^[i]:=Q^[i]*norm; END;
<* IF (__debug__) THEN *>
WrVec("Q_c",Q^,M);
<* END *>
(* step 5 - compute R = A^{tr} x Q *)
norm:=0.0;
FOR i:=0 TO N-1 DO
Ri:=0.0;
FOR j:=0 TO M-1 DO
Ri:=Ri + A[i,j]*Q^[j];
END;
R^[i]:=Ri;
norm:=norm + Ri*Ri;
END;
s:=sqrt(norm);
norm := 1.0 / sqrt(norm);
FOR i:=0 TO N-1 DO R^[i]:=R^[i]*norm; END;
<* IF (__debug__) THEN *>
WrVec("R_c",R^,N);
<* END *>
(* step 6 *)
d:=0.0;
FOR i:=0 TO N-1 DO
dPRi:=ABS(P^[i] - R^[i]);
IF (dPRi > d) THEN d:=dPRi; END;
END;
IF (d > delta) THEN
IF (d > 10.0*eps) THEN
FOR i:=0 TO N-1 DO
W^[i] := P^[i] - R^[i];
END;
END;
FOR i:=0 TO N-1 DO
P^[i] := R^[i];
END;
END;
UNTIL (d <= delta) OR (counter > climit); (* goto step 3 *)
IF (counter > climit) THEN
iFehl:=2;
Fehlerflag:="climit exceeded (PowSVD)";
END;
IF ((s / (sv[0] + tol)) <= tol) THEN (* step 7 *)
rankdef:=TRUE;
iFehl:=iFehl + 1;
Fehlerflag:="Matrix rank deficite (PowSVD)";
ELSE
sv[k] := s; (* step 8 - convergence achieved *)
FOR i:=0 TO M-1 DO U[k,i]:=Q^[i]; END;
FOR i:=0 TO N-1 DO V[k,i]:=R^[i]; END;
INC(k);
IF (k < klimit) THEN (* step 10 *)
FOR i:=0 TO N-1 DO
Ri:=R^[i];
FOR j:=0 TO M-1 DO (* deflate matrix A *)
A[i,j]:=A[i,j] - s*Q^[j]*Ri;
END;
END;
(* step 11 *)
norm:=0.0;
FOR i:=0 TO N-1 DO norm:=norm + sqr(W^[i]); END;
norm:=sqrt(norm);
IF (norm < eps) THEN (* goto step 1 *)
(* Only "noise" in current W vector - get new values *)
norm:=0.0;
FOR i:=0 TO N-1 DO
W^[i]:=Zufall(); norm:=norm + sqr(W^[i]);
END;
END;
FOR i:=0 TO N-1 DO P^[i]:=W^[i] / norm; END;
END; (* IF k *)
END; (* IF rankdef *)
(* goto step 2; *)
UNTIL (k >= klimit) OR (rankdef);
rank:=k;
DEALLOCATE(P,N*TSIZE(LONGREAL));
DEALLOCATE(Q,M*TSIZE(LONGREAL));
DEALLOCATE(R,N*TSIZE(LONGREAL));
DEALLOCATE(W,N*TSIZE(LONGREAL));
END PowSVD;
PROCEDURE JacobiSVD( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
job : CARDINAL;
VAR rank : CARDINAL;
VAR iFehl : INTEGER);
CONST eps = MachEps;
tol = 0.1*MachEps;
precise = FALSE;
VAR nt,i,j,k : INTEGER;
scount,rcount,slim : INTEGER;
p,q,r : LONGREAL;
IfV,IfU,rot,IfRank : BOOLEAN;
g,h,c,s,tau,vt,e2 : LONGREAL;
BEGIN
(* step 0 - initialisation *)
IfRank := (rank = VAL(CARDINAL,2*(M+N)));
<* IF (__debug__) THEN *>
TIO.WrStr("IfRank = "); TIO.WrBool(IfRank); TIO.WrLn;
<* END *>
iFehl:=0;
c:=0.0; s:=0.0; (* Wg. Compilerwarnung *)
slim := N DIV 4; IF (slim < 12) THEN slim := 12; END;
nt:=N; (* current estimate of rank *)
scount:=0;
IfV:=FALSE; IfU:=FALSE;
IF (job = 1) THEN
IfV:=TRUE;
ELSIF (job = 2) THEN
IfU:=TRUE;
ELSIF (job = 3) THEN
IfV:=TRUE; IfU:=TRUE;
END;
IF IfV THEN (* step 1 *)
FOR i:=0 TO N-1 DO
FOR j:=0 TO N-1 DO V[i,j]:=0.0; END;
V[i,i]:=1.0;
END;
END;
e2 := 10.0*LFLOAT(M)*eps*eps;
REPEAT (* step 2 *)
rcount:=nt*(nt -1) DIV 2;
INC(scount);
FOR i:=0 TO nt-2 DO (* step 3 *)
FOR k:=i+1 TO nt-1 DO (* step 4 *)
p:=0.0; q:=0.0; r:=0.0; (* step 5 *)
IF NOT precise THEN
FOR j:=0 TO M-1 DO
p:=p + A[i,j]*A[k,j];
q:=q + A[i,j]*A[i,j];
r:=r + A[k,j]*A[k,j];
END;
ELSE (* Summation mit "Gard Digit" *)
p := NeumaierProdSum(A[i],A[k],M); (* ddot *)
q := NeumaierProdSum(A[i],A[i],M); (* dnrm2 *)
r := NeumaierProdSum(A[k],A[k],M); (* dnrm2 *)
END;
W[i]:=q; (* step 6 *)
W[k]:=r;
rot:=FALSE;
(* if (q >= r) then goto step 9 *)
IF (q < r) THEN
q := (q / r) - 1.0; (* step 8 *)
p := (p / r);
vt := Pythag(2.0*p,q);
s := sqrt(0.5*(1.0 - (q / vt)));
IF (p < 0.0) THEN s:=-s; END;
c := p / (vt*s);
rot:=TRUE; (* GOTO step 11; *)
ELSE (* step 9 *)
(* if (q*r <= eps*eps) then goto step 13 *)
(* if ((p / q)*(p / r) > eps) then goto step 13 *)
(* IF (q*r > eps*eps) AND ((p / q)*(p / r) > eps) THEN *)
(* Possible alternative (taken from ref[2] *)
(*
IF NOT ((r = 0.0) OR ((ABS(q) < eps) AND (ABS(r) < eps)) OR
((p/q)*(p/r) < eps))
THEN
*)
IF NOT ((q <= e2*W[0]) OR (ABS(p) <= tol*q)) THEN
(* step 10 *)
r := 1.0 - (r / q);
p := p / q;
vt := Pythag(2.0*p,r);
c := sqrt(0.5*(1.0 + (r / vt)));
s := p / (vt*c);
rot:=TRUE;
END;
END;
IF rot THEN
IF NOT precise THEN
FOR j:=0 TO M-1 DO (* step 11 *)
g := A[i,j];
h := A[k,j];
A[i,j] := c*g + s*h;
A[k,j] := - s*g + c*h;
END;
IF IfV THEN
FOR j:=0 TO N-1 DO (* step 12 *)
g := V[i,j];
h := V[k,j];
V[i,j] := c*g + s*h;
V[k,j] := - s*g + c*h;
END;
(* goto step 14 *)
END;
ELSE (* MRi, 04.12.17 *)
tau := s / (1.0 + c);
FOR j:=0 TO M-1 DO (* step 11 *)
g := A[i,j];
h := A[k,j];
A[i,j] := g + s*(h - g*tau);
A[k,j] := h - s*(g + h*tau);
END;
IF IfV THEN
FOR j:=0 TO N-1 DO (* step 12 *)
g := V[i,j];
h := V[k,j];
V[i,j] := g + s*(h - g*tau);
V[k,j] := h - s*(g + h*tau);
END;
(* goto step 14 *)
END;
END; (* IF precise *)
ELSE (* decrement rot. count since rot. was not necessary *)
DEC(rcount); (* step 13 *)
END; (* IF *) (* step 14 *)
END; (* FOR k *) (* step 15 *)
END; (* FOR i *) (* step 16 *)
<* IF (__debug__) THEN *>
TIO.WrStr("End of sweep "); TIO.WrInt(scount,2);
TIO.WrStr(", number of rotations performed "); TIO.WrInt(rcount,3);
TIO.WrLn;
<* END *>
(* Remove the next "while" statement to force higher precision *)
(* computation in cases of severe multicollinearity. Otherwise *)
(* do NOT use results in cases where Z[N]/Z[1] is very small. *)
IF NOT precise AND IfRank THEN
WHILE (nt >= 3) AND ((W[nt-1] / (W[0] + eps)) <= eps) DO
DEC(nt); (* step 18 *)
iFehl:=1;
Fehlerflag:="Matrix rank deficite (JacobiSVD)";
END;
END;
UNTIL (rcount = 0) OR (scount > slim); (* step 19 & 2 (part) *)
IF (scount > slim) THEN (* step 2 (part) *)
iFehl:=iFehl + 2;
Fehlerflag:="too many sweeps (JacobiSVD)";
END;
rank:=nt;
(* Normalize Vektor U *)
FOR i:=0 TO nt-1 DO (* q is Norm of U[i] if W[i] # 0 *)
IF (ABS(W[i]) > Small) THEN q:=sqrt(ABS(W[i])); ELSE q:=0.0; END;
W[i]:=q;
IF IfU THEN
IF (ABS(q) = 0.0) THEN
s:=0.0;
FOR j:=0 TO M-1 DO s:=A[i,j]*A[i,j]; END;
IF (s # 0.0) THEN q:=1.0 / sqrt(s); ELSE q:=0.0; END;
FOR j:=0 TO M-1 DO A[i,j] := A[i,j] * q; END;
ELSE
FOR j:=0 TO M-1 DO A[i,j] := A[i,j] / q; END;
END;
END;
END;
END JacobiSVD;
PROCEDURE CalcLQsol1( M,N : INTEGER;
rank : INTEGER;
IfT : CHAR;
VAR U : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
tol : LONGREAL);
VAR i,j,k : INTEGER;
BEGIN
IF (tol <= 0.0) THEN tol:=100.0*W[0]*MachEps; END;
tol := tol*tol;
IF (CAP(IfT) = "T") THEN
FOR i:=0 TO N-1 DO X[i]:=0.0; END;
FOR i:=0 TO rank-1 DO (* reduced rank < N ??? *)
FOR j:=0 TO N-1 DO
FOR k:=0 TO M-1 DO
IF (W[i] > tol) THEN (* ignore "small" SVs *)
X[j]:=X[j] + V[i,j]*U[i,k]*B[k] / (W[i]);
END;
END;
END;
END;
ELSE (* Use memory order from dSVD *)
FOR i:=0 TO N-1 DO X[i]:=0.0; END;
FOR i:=0 TO rank-1 DO (* reduced rank < N ??? *)
FOR j:=0 TO N-1 DO
FOR k:=0 TO M-1 DO
IF (W[i] > tol) THEN (* ignore "small" SVs *)
X[j]:=X[j] + V[i,j]*U[k,i]*B[k] / (W[i]);
END;
END;
END;
END;
END;
END CalcLQsol1;
PROCEDURE CalcLQres1( M,N : INTEGER;
IfT : CHAR;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR R : ARRAY OF LONGREAL;
VAR r1 : LONGREAL;
VAR r2 : LONGREAL);
VAR i,j : INTEGER;
Ri,s1,s2 : LONGREAL;
BEGIN
s1:=0.0; s2:=0.0;
IF (CAP(IfT) = "Y") OR (CAP(IfT) = "T") THEN
FOR i:=0 TO M-1 DO R[i] := - B[i]; END;
FOR i:=0 TO N-1 DO
FOR j:=0 TO M-1 DO R[j]:=R[j] + A[i,j]*X[i]; END;
END;
FOR i:=0 TO M-1 DO s1:=s1 + R[i]; s2:=s2 + R[i]*R[i]; END;
ELSE (* Use memory order from dSVD *)
FOR i:=0 TO M-1 DO
Ri := - B[i]; (* note form of residual is R = A * X - B *)
FOR j:=0 TO N-1 DO Ri:=Ri + A[i,j]*X[j]; END;
R[i]:=Ri;
s1:=s1 + Ri;
s2:=s2 + Ri*Ri;
END;
END;
r1 := s1;
r2 := s2 / LFLOAT(M);
END CalcLQres1;
PROCEDURE GivSVD( M,N : INTEGER;
nRHS : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF ARRAY OF LONGREAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR rss : ARRAY OF LONGREAL;
VAR svs : ARRAY OF LONGREAL;
VAR W : ARRAY OF ARRAY OF LONGREAL;
tol : LONGREAL;
VAR iFehl : INTEGER);
CONST eps = MachEps;
(* Perform the rotations, will be inlined by compiler *)
PROCEDURE RotNSub(l,j,k,tcol : INTEGER;
c,s : LONGREAL);
VAR i : INTEGER;
r : LONGREAL;
BEGIN
FOR i:=l TO tcol-1 DO
r := W[j,i];
W[j,i] := r*c + s*W[k,i];
W[k,i] := -r*s + c*W[k,i];
END;
END RotNSub;
VAR count,EstRowRank,slimit,sweep,tcol : INTEGER;
i,j,k : INTEGER;
bb,c,e2,p,q,r,s,vt,tol2,trss : LONGREAL;
BEGIN
iFehl :=0;
tcol := N+nRHS;
k := N; (* Eventuell k durch N ersetzen bis Schleife FOR k ... *)
FOR i:=0 TO N-1 DO
FOR j:=0 TO tcol-1 DO
W[i,j]:=0.0;
END;
END;
FOR i:=0 TO nRHS-1 DO rss[i]:=0.0; END;
tol2 := FLOAT(N*N)*eps*eps;
<* IF (__debug__) THEN *>
TIO.WrLn; TIO.WrStr("Matrix (A|B)"); TIO.WrLn; TIO.WrLn;
<* END *>
FOR i:=0 TO M-1 DO
(* Here we could call a routine get a new row of observations *)
(* and not handing them over per parameter. Then we need a *)
(* procedure parameter which also hands over a paramete which *)
(* indicated end of data (insted of a fixed line count "M") *)
(* GetObsN(N,nRHS,W,k,endflag); *)
FOR j:=0 TO N-1 DO (* copy line of A into work Array W *)
W[k,j]:=A[i,j];
END;
FOR j:=0 TO nRHS-1 DO (* copy RHS into work Array W *)
W[k,j+N]:=B[j,i]; (* !!! *)
END;
<* IF (__debug__) THEN *>
FOR j:=0 TO N+nRHS-1 DO TIO.WrLngReal(W[k,j],12,5); END; TIO.WrLn;
<* END *>
FOR j:=0 TO N-1 DO (* step 3 *)
(* loop over the rows of the work array to move information *)
(* into the triangular part of the Givens reduction *)
(* step 4 *)
s:=W[k,j];
c:=W[j,j];
bb := ABS(c); IF (ABS(s) > bb) THEN bb := ABS(s); END;
IF (bb > 0.0) THEN
c := c / bb;
s := s / bb;
p := sqrt(c*c + s*s); (* step 7 *)
s := s / p;
IF (ABS(s) >= tol2) THEN (* step 8 *)
c := c/p;
RotNSub(j,j,k,tcol,c,s);
END;
END; (* FOR j *)
END;
FOR j:=0 TO nRHS-1 DO
(* step 10 - accumulate the residual sums of squares *)
rss[j]:=rss[j] + sqr(W[k,N+j]);
END;
END; (* FOR i=1,nobs *)
<* IF (__debug__) THEN *>
TIO.WrLn; TIO.WrStr('Uncorrelated residuals'); TIO.WrLn;
FOR i:=0 TO nRHS-1 DO TIO.WrLngReal(rss[i],12,-3); END; TIO.WrLn;
TIO.WrLn; TIO.WrStr('Matrix W[1:n+1,1:n+nRHS]'); TIO.WrLn; TIO.WrLn;
FOR i:=0 TO N+1-1 DO
FOR j:=0 TO N+nRHS-1 DO TIO.WrLngReal(W[i,j],12,-3); END; TIO.WrLn;
END; TIO.WrLn;
<* END *>
slimit := N DIV 4; IF (slimit < 6) THEN slimit := 6;END;
sweep := 0;
e2 := 10.0*FLOAT(N)*eps*eps;
EstRowRank := N;
REPEAT
count := 0;
FOR j:=0 TO EstRowRank-2 DO
FOR k:=j+1 TO EstRowRank-1 DO
p := 0.0; q:=0.0; r:=0.0;
FOR i:=0 TO N-1 DO
p:=p + W[j,i]*W[k,i];
q:=q + sqr(W[j,i]);
r:=r + sqr(W[k,i]);
END;
svs[j] := q; svs[k] := r;
IF (q >= r) THEN
IF NOT ((q <= e2*svs[0]) OR (ABS(p) <= tol*q)) THEN
p := p/q;
r := 1.0 - r/q;
vt := sqrt(4.0*p*p + r*r);
c := sqrt(0.5*(1.0 + r/vt));
s := p / (vt*c);
RotNSub(0,j,k,tcol,c,s);
count:=count+1;
END;
ELSE
p := p/r;
q := q/r - 1.0;
vt := sqrt(4.0*p*p + q*q);
s := sqrt(0.5*(1.0 - q/vt));
IF (p < 0) THEN s := -s;END;
c := p / (vt*s);
RotNSub(0,j,k,tcol,c,s);
INC(count);
END;
END;
END;
INC(sweep);
(* Remove the next "while" statement to force higher precision *)
(* computation in cases of severe multicollinearity. Otherwise do *)
(* NOT use results in cases where svs[n]/svs[1] is very small. *)
WHILE (EstRowRank >= 3) AND (svs[EstRowRank-1] <= svs[0]*tol+tol*tol) DO
EstRowRank:=EstRowRank - 1;
iFehl:=1;
END;
UNTIL (sweep > slimit) OR (count = 0);
FOR j:=0 TO N-1 DO
s := svs[j];
s := sqrt(s);
svs[j] := s; (* Singular value *)
IF (s >= tol) THEN
FOR i:=0 TO N-1 DO W[j,i]:=W[j,i] / s; END;
END;
END;
<* IF (__debug__) THEN *>
TIO.WrLn; TIO.WrStr("Matrix W'[1:n+1,1:n+nRHS]"); TIO.WrLn; TIO.WrLn;
FOR i:=0 TO N+1-1 DO
FOR j:=0 TO N+nRHS-1 DO TIO.WrLngReal(W[i,j],12,-3); END;
TIO.WrLn;
END; TIO.WrLn;
TIO.WrLn; TIO.WrStr("B^tr"); TIO.WrLn; TIO.WrLn;
FOR i:=0 TO nRHS-1 DO
FOR j:=0 TO N-1 DO
TIO.WrStr('W['); TIO.WrInt(j,2); TIO.WrChar(',');
TIO.WrInt((N+i),2);
TIO.WrStr('] = '); TIO.WrLngReal(W[j,N+i],14,6); TIO.WrLn;
END;
END;
TIO.WrLn;
<* END *>
FOR i:=0 TO nRHS-1 DO
trss := rss[i];
FOR j:=0 TO N-1 DO
p:=0.0;
FOR k:=0 TO N-1 DO
IF (svs[k] > Small) THEN p:=p + W[k,j]*W[k,N+i] / svs[k]; END;
END;
X[i,j] := p;
IF (svs[j] <= MachEps) THEN trss:=trss + sqr(W[j,N+i]); END;
END;
rss[i]:=trss;
END;
END GivSVD;
(*=========================== Eispack Routinen =============================*)
PROCEDURE dSVD(VAR A : ARRAY OF ARRAY OF LONGREAL;
m,n : INTEGER;
MatU : CARDINAL;
VAR W : ARRAY OF LONGREAL;
MatV : CARDINAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR IErr : INTEGER);
VAR i,j,k,l,ll,its,jj,nm : INTEGER;
c,f,h,s,x,y,z : LONGREAL;
anorm,g,scale : LONGREAL;
split : BOOLEAN;
RV1 : POINTER TO
ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
WantU,WantV : BOOLEAN;
BEGIN
Errors.Fehler:=FALSE;
IErr:=0;
IF (m < n) THEN
Errors.Fehlerflag:='WARNUNG M < N (dSVD)';
Errors.ErrOut(Errors.Fehlerflag);
Errors.Fehler:=TRUE;
IErr:=-1;
END;
WantU := MatU = 1;
WantV := MatV = 1;
anorm := 0.0;
scale := 0.0;
g := 0.0;
ALLOCATE(RV1,n*TSIZE(LONGREAL));
(* Householder reduction to bidiagonal form *)
l:=n;
FOR i:=0 TO n-1 DO (* left-hand reduction *)
l := i + 1;
RV1^[i] := scale*g;
g := 0.0; s := 0.0; scale := 0.0;
IF (i < m) THEN
FOR k:=i TO m-1 DO
scale := scale + ABS(A[k,i]); (* <= *)
END;
IF (scale # 0.0) THEN
FOR k:=i TO m-1 DO
A[k,i]:=A[k,i] / scale; (* <= *)
s:=s + A[k,i]*A[k,i]; (* <= *)
END;
f := A[i,i];
g := -dsign(sqrt(s), f);
h := f*g - s;
A[i,i] := (f - g);
IF (i # n-1) THEN
FOR j:=l TO n-1 DO
s:=0.0;
FOR k:=i TO m-1 DO;
s:=s + A[k,i]*A[k,j]; (* <= *)
END;
f := s / h;
FOR k:=i TO m-1 DO
A[k,j]:=A[k,j] + f*A[k,i]; (* <= *)
END;
END;
END;
FOR k:=i TO m-1 DO
A[k,i] := A[k,i]*scale; (* <= *)
END;
END;
END;
W[i] := scale*g;
(* right-hand reduction *)
g:=0.0; s:=0.0; scale:=0.0;
IF (i < m) AND (i # n-1) THEN
FOR k:=l TO n-1 DO
scale:=scale + ABS(A[i,k]); (* !! *)
END;
IF (scale # 0.0) THEN
FOR k:=l TO n-1 DO
A[i,k]:=A[i,k] / scale; (* !! *)
s:=s + A[i,k]*A[i,k]; (* !! *)
END;
f := A[i,l];
g := -dsign(sqrt(s),f);
h := f*g - s;
A[i,l] := (f - g);
FOR k:=l TO n-1 DO
RV1^[k] := A[i,k] / h; (* !! *)
END;
IF (i # m-1) THEN (* diff *)
FOR j:=l TO m-1 DO
s := 0.0;
FOR k:=l TO n-1 DO s:=s + (A[j,k]*A[i,k]); END; (* !! O(nm) *)
FOR k:=l TO n-1 DO A[j,k]:=A[j,k] + (s*RV1^[k]); END; (* !! n**2 O(nm)*)
END;
END;
FOR k:=l TO n-1 DO A[i,k]:=A[i,k]*scale; END; (* !! *)
END;
END;
anorm := dmax(anorm, (ABS(W[i]) + ABS(RV1^[i])));
END; (* FOR i *)
(* accumulate the right-hand transformation *)
IF WantV THEN
FOR i:=n-1 TO 0 BY -1 DO
IF (i < n-1) THEN
IF (g # 0.0) THEN
FOR j:=l TO n-1 DO (* double division to avoid underflow *)
V[i,j] := ((A[i,j] / A[i,l]) / g); (* !! *)
END;
FOR j:=l TO n-1 DO
s := 0.0;
FOR k:=l TO n-1 DO s:=s + A[i,k]*V[j,k]; END; (* !! *)
FOR k:=l TO n-1 DO V[j,k]:=V[j,k] + s*V[i,k]; END;
END;
END;
FOR j:=l TO n-1 DO V[i,j]:=0.0; V[j,i]:=0.0; END;
END;
V[i,i] := 1.0;
g := RV1^[i];
l := i;
END; (* FOR i *)
END; (* IF WantV *)
(* accumulate the left-hand transformation *)
IF WantU THEN
FOR i:=n-1 TO 0 BY -1 DO
l := i + 1;
g := W[i];
IF (i < n-1) THEN
FOR j:=l TO n-1 DO
A[i,j] := 0.0; (* !! *)
END;
END;
IF (g # 0.0) THEN
g := 1.0 / g;
IF (i # n-1) THEN (* diff *)
FOR j:=l TO n-1 DO
s := 0.0;
FOR k:=l TO m-1 DO s:=s + (A[k,i]*A[k,j]); END; (* <= *)
f := (s / A[i,i])*g;
FOR k:=i TO m-1 DO A[k,j]:=A[k,j] + f*A[k,i]; END; (* <= *)
END;
END;
FOR j:=i TO m-1 DO A[j,i]:=A[j,i]*g; END; (* <= *)
ELSE
FOR j:=i TO m-1 DO A[j,i]:=0.0; END; (* <= *)
END; (* IF g *)
A[i,i]:=A[i,i] + 1.0; (* .. *)
END; (* FOR i *)
END; (* IF WantU *)
(* diagonalize the bidiagonal form *)
FOR k:=n-1 TO 0 BY -1 DO (* loop over singular values *)
its:=0;
LOOP (* loop over allowed iterations *)
INC(its);
split := FALSE;
LOOP
ll:=0;
nm:=0;
FOR l:=k TO 0 BY -1 DO (* test for splitting *)
nm := l - 1;
IF ((ABS(RV1^[l]) + anorm) = anorm) THEN
ll:=l;
split := TRUE;
EXIT;
END;
IF ((ABS(W[nm])+anorm) = anorm) THEN
ll:=l;
EXIT;
END;
END;
EXIT;
END;
l:=ll;
IF NOT split THEN
s:=1.0; (* c:=0.0; *)
FOR i:=l TO k DO
f := s*RV1^[i];
IF ((ABS(f) + anorm) # anorm) THEN
g := W[i];
h := Pythag(f, g);
W[i] := h;
h := 1.0 / h;
c := g*h;
s := (- f*h);
IF WantU THEN
FOR j:=0 TO m-1 DO
y := A[j,nm]; (* <= *)
z := A[j,i]; (* <= *)
A[j,nm] := (y*c + z*s); (* <= *)
A[j,i] := (z*c - y*s); (* <= *)
END;
END;
END;
END; (* FOR i *)
END; (* IF *)
z := W[k];
IF (l = k) THEN (* convergence *)
IF (z < 0.0) THEN (* make singular value nonnegative *)
W[k] := (-z);
IF WantV THEN
FOR j:=0 TO n-1 DO
V[k,j] := (-V[k,j]);
END;
END;
END;
EXIT; (* !!! LOOP its !!! *)
END;
IF (its >= 30) THEN
IF Errors.Fehler THEN
IErr := -VAL(INTEGER,k+1);
ELSE
IErr := VAL(INTEGER,k+1);
END;
Errors.Fehlerflag:='No convergence after 30 iterations ';
Errors.Fehler:=TRUE;
RETURN;
END;
(* shift from bottom 2 x 2 minor *)
x := W[l];
nm := k - 1;
y := W[nm];
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 + dsign(g, f))) - h)) / x;
(* next QR transformation *)
c:=1.0; s:=1.0;
FOR j:=l TO nm DO
i := j + 1;
g := RV1^[i];
y := W[i];
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 := g*c - x*s;
h := y*s;
y := y*c;
IF WantV THEN
FOR jj:=0 TO n-1 DO
x := V[j,jj];
z := V[i,jj];
V[j,jj] := (x*c + z*s);
V[i,jj] := (z*c - x*s);
END;
END;
z := Pythag(f, h);
W[j] := z;
IF (z # 0.0) THEN
z := 1.0 / z;
c := f*z;
s := h*z;
END;
f := (c*g) + (s*y);
x := (c*y) - (s*g);
IF WantU THEN
FOR jj:=0 TO m-1 DO
y := A[jj,j]; (* <= *)
z := A[jj,i]; (* <= *)
A[jj,j] := (y*c + z*s); (* <= *)
A[jj,i] := (z*c - y*s); (* <= *)
END;
END;
END; (* FOR j *)
RV1^[l] := 0.0;
RV1^[k] := f;
W[k] := x;
END; (* LOOP its *)
END; (* FOR k *)
DEALLOCATE(RV1,n*TSIZE(LONGREAL));
END dSVD;
PROCEDURE OrderSVD(VAR U : ARRAY OF ARRAY OF LONGREAL;
VAR W : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
m,n : CARDINAL);
VAR i,j,MaxI : CARDINAL;
zw,MaxW : LONGREAL;
BEGIN
FOR i:=0 TO n-1 DO
MaxW := W[i]; MaxI := i;
FOR j:=i+1 TO n-1 DO
IF (W[j] > MaxW) THEN MaxW:=W[j]; MaxI:=j; END;
END;
IF (i # MaxI) THEN
zw := W[i]; W[i]:=MaxW; W[MaxI]:= zw;
FOR j:=0 TO n-1 DO
zw := V[i,j];
V[i,j] := V[MaxI,j];
V[MaxI,j] := zw;
END;
FOR j:=0 TO m-1 DO
zw := U[j,i];
U[j,i] := U[j,MaxI];
U[j,MaxI] := zw;
END;
END; (* IF *)
END;
END OrderSVD;
PROCEDURE SVDSolv(VAR U : ARRAY OF ARRAY OF LONGREAL; (* m,n *)
Utr : CHAR;
VAR W : ARRAY OF LONGREAL; (* n *)
VAR V : ARRAY OF ARRAY OF LONGREAL; (* 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 [0..8190] OF LONGREAL; (* n *)
Svek : POINTER TO ARRAY [0..8190] 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
Errors.FatalError("Kein Freispeicher vorhanden (SVDLib0.SVDSolv) !");
END;
FOR j:=0 TO N-1 DO
s := 0.0;
IF (W[j] # 0.0) THEN
IF (CAP(Utr) # "T") THEN
FOR i:=0 TO M-1 DO Svek^[i] := U[i,j]*C[i]; END;
ELSE
FOR i:=0 TO M-1 DO Svek^[i] := U[j,i]*C[i]; END;
END;
s := SumVek(Svek^,1,M);
s := s / W[j];
END;
Zw^[j] := s;
END;
FOR j:=0 TO N-1 DO
FOR jj:=0 TO N-1 DO Svek^[jj] := V[jj,j]*Zw^[jj]; END;
s := SumVek(Svek^,1,N);
X[j] := s;
END;
DEALLOCATE(Svek,MaxCard(N,M)*TSIZE(LONGREAL));
DEALLOCATE(Zw,N*TSIZE(LONGREAL));
END SVDSolv;
PROCEDURE MinFit( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
Nb : INTEGER;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* Adaption nach Modula-2 von Michael Riedl, July 2016 *)
(*----------------------------------------------------------------*)
PROCEDURE CacheLen(Ld : INTEGER) : INTEGER;
(*-----------------------------------------------------*)
(* Estimated number of columns of a longreal matrix *)
(* fitting in the cache. Ld is first dimension of the *)
(* matrix array *)
(*-----------------------------------------------------*)
BEGIN
RETURN 2048 / (8*MAX0(Ld,1));
END CacheLen;
PROCEDURE DSC16( S : LONGREAL;
VAR S16 : LONGREAL;
VAR R16 : LONGREAL);
(*-----------------------------------------------------*)
(* Integer power of 16 "near" a positive scalar and *)
(* reciprocal *)
(* *)
(* s : positive scalar *)
(* s16 : integer power of 16 "near" s *)
(* r16 : reciprocal of s16 *)
(*-----------------------------------------------------*)
BEGIN
S16:=1.0;
WHILE (S16 < ABS(S)) DO S16:=S16*16.0; END;
R16:=1.0 / S16;
END DSC16;
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST maxit = 30;
VAR i,j,k,kk,l,ll,mn,it,kache : INTEGER;
kp1,ip1 :INTEGER;
c,s,x,y,z,f,g,h,bnorm,s16,r16 : LONGREAL;
cancellation,exit : BOOLEAN;
E : POINTER TO ARRAY
[0..MAX(INTEGER)-1] OF LONGREAL;
BEGIN
ALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
iErr := 0;
E^[0] := 0.0;
mn := MIN0(M,N);
kache := CacheLen(HIGH(A));
kk := N - kache;
(* reduction to upper-bidiagonal form by householder transformations *)
kp1:=1;
LOOP (* FOR k:=1 TO mn DO *)
k:=kp1-1;
(* IF (k > mn) THEN EXIT; END; *)
(* L1 norm of the subdiagonal part of column k *)
x:=0.0;
FOR i:=kp1 TO M-1 DO
x:=x + ABS(A[k,i]);
END;
IF (x # 0.0) THEN
x:=x + ABS(A[k,k]);
(* scaled householder vector *)
DSC16(x,s16,r16);
s:=0.0;
FOR i:=k TO M-1 DO
A[k,i] := r16*A[k,i];
s:=s + sqr(A[k,i]);
END;
s := -DSIGN(sqrt(s),A[k,k]);
x := 1.0 / s;
<* IF (__BLAS__) THEN *>
dscal(M-kp1+1,x,A[k,k],1);
<* ELSE *>
FOR i:=k TO M-1 DO
A[k,i] := x*A[k,i];
END;
<* END *>
A[k,k] := A[k,k] - 1.0;
(* diagonal element *)
D[k] := s16*s;
(* left transformation *)
x := 1.0 / A[k,k];
IF (kp1 < kk) THEN
(* path for small arrays *)
FOR j:=k+1 TO N-1 DO
s:=0.0;
FOR i:=k TO M-1 DO
s:=s + A[j,i]*A[k,i];
END;
s := s*x;
FOR i:=k TO M-1 DO
A[j,i]:=A[j,i] + s*A[k,i];
END;
E^[j] := A[j,k];
END; (* FOR *)
ELSE
(* path for large arrays *)
FOR j:=k+1 TO N-1 DO
<* IF (__BLAS__) THEN *>
s := x*ddot(M-kp1+1,A[j,k],1,A[k,k],1);
daxpy(M-kp1+1,s,A[k,k],1,A[j,k],1);
<* ELSE *>
s:=0.0;
FOR i:=k TO M-1 DO
s:=s + A[j,i]*A[k,i];
END;
s := s*x;
FOR i:=k TO M-1 DO
A[j,i] := A[j,i] + s*A[k,i];
END;
<* END *>
E^[j] := A[j,k];
END; (* FOR j *)
END; (* IF k < kk *)
(* transformation of the right-hand side *)
FOR j:=0 TO Nb-1 DO
s:=0.0;
FOR i:=k TO M-1 DO
s:=s + A[k,i]*B[i,j]; (* !!! *)
END;
s := s*x;
FOR i:=k TO M-1 DO
B[i,j]:=B[i,j] + s*A[k,i]; (* !!! *)
END;
END; (* FOR j *)
ELSE
(* the transformation is bypassed if the norm is zero *)
D[k] := A[k,k];
FOR j:=k+1 TO N-1 DO
E^[j] := A[j,k];
END;
END; (* IF x # 0 *)
IF (kp1 = N) THEN
EXIT; (* LOOP k *)
END;
(* L1 norm of the supercodiagonal part of row k *)
x:=0.0;
FOR j:=k+2 TO N-1 DO
x:=x + ABS(E^[j]);
END;
IF (x # 0.0) THEN
(* scaled householder vector *)
x:=x + ABS(E^[kp1]);
DSC16(x,s16,r16);
s:=0.0;
FOR j:=k+1 TO N-1 DO
E^[j] := r16*E^[j];
s:=s + sqr(E^[j]);
END;
y := -DSIGN(sqrt(s),E^[kp1]);
E^[kp1] := E^[kp1] - y;
x := 1.0 / y;
FOR j:=k+1 TO N-1 DO
E^[j] := x*E^[j];
A[k,j] := E^[j];
END;
x := 1.0 / E^[kp1];
(* right transformation *)
FOR i:=kp1 TO M-1 DO
s:=0.0;
FOR j:=k+1 TO N-1 DO
s:=s + A[j,i]*E^[j];
END;
s := s*x;
FOR j:=k+1 TO N-1 DO
A[j,i]:=A[j,i] + E^[j]*s;
END;
END; (* FOR *)
(* codiagonal element *)
E^[kp1] := s16*y;
ELSE (* x = 0 *)
(* the transformation is bypassed if the norm is zero *)
FOR i:=kp1 TO N-1 DO
A[k,i]:=0.0;
END;
END; (* IF x *)
INC(kp1);
END; (* LOOP k *)
l := mn;
IF (M < N) THEN
INC(l);
END; (* IF *)
(* right matrix of the bidiagonal decomposition *)
FOR j:=l TO N-1 DO
FOR i:=0 TO N-1 DO
A[j,i]:=0.0;
END;
A[j,j]:=1.0;
END; (* FOR j *)
kk := N - kache;
FOR kp1:=l TO 2 BY -1 DO
k:=kp1-1;
FOR i:=0 TO k-1 DO
A[k,i] := 0.0;
END;
FOR i:=k TO N-1 DO
A[k,i] := A[k-1,i];
END;
IF (A[k,k] # 0.0) THEN
x := 1.0 / A[k,k];
IF (kp1 < kk) THEN
(* path for small arrays *)
FOR j:=k+1 TO N-1 DO
s:=0.0;
FOR i:=k TO N-1 DO
s:=s + A[k,i]*A[j,i];
END;
s := s*x;
FOR i:=k TO N-1 DO
A[j,i]:=A[j,i] + s*A[k,i];
END;
END; (* FOR *)
ELSE
(* path for large arrays *)
FOR j:=k+1 TO N-1 DO
<* IF (__BLAS__) THEN *>
s := x*ddot(N-kp1+1,A[k,k],1,A[j,k],1);
daxpy(N-kp1+1,s,A[k,k],1,A[j,k],1);
<* ELSE *>
s := 0.0;
FOR i:=k TO N-1 DO
s:=s + A[k,i]*A[j,i];
END;
s := s*x;
FOR i:=k TO N-1 DO
A[j,i]:=A[j,i] + s*A[k,i];
END;
<* END *>
END; (* FOR j *)
END; (* IF k < kk *)
END; (* IF A[k-1,k-1] # 0 *)
A[k,k]:=A[k,k] + 1.0;
END; (* FOR k *)
A[0,0]:=1.0;
FOR i:=1 TO N-1 DO
A[0,i]:=0.0;
END;
(* annihilation of the last superdiagonal element when m < n *)
IF (M < N) THEN
c := 1.0;
s := -1.0;
kp1:=mn; exit:=FALSE;
WHILE (kp1 > 1) AND NOT exit DO (* FOR k:=mn TO 1 BY -1 DO *)
x := -s*E^[kp1];
E^[kp1] := c*E^[kp1];
IF (x = 0.0) THEN
exit:=TRUE;
ELSE
IF (ABS(D[k]) < ABS(x)) THEN
c := D[k] / x;
y := sqrt(1.0 + c*c);
s := 1.0 / y;
c := c*s;
D[k] := y*x;
ELSE
s := x / D[k];
y := sqrt(1.0 + s*s);
c := 1.0 / y;
s := s*c;
D[k] := y*D[k];
END; (* IF *)
FOR i:=0 TO N-1 DO
x := A[k,i];
y := A[l-1,i];
A[k,i] := c*x + s*y;
A[l-1,i] := c*y - s*x;
END; (* FOR *)
END; (* IF *)
DEC(kp1);
END; (* WHILE k *)
END; (* IF M < N *)
FOR i:=mn TO N-1 DO
FOR j:=0 TO Nb-1 DO
B[i,j]:=0.0;
END;
D[i]:=0.0;
END;
(* singular-value decomposition of the upper- *)
(* bidiagonal matrix of order min(m,n) *)
(* L1 norm of the bidiagonal matrix *)
bnorm := 0.0;
FOR i:=0 TO mn-1 DO
bnorm := DMAX(ABS(D[i])+ABS(E^[i]),bnorm);
END;
FOR k:=mn TO 1 BY -1 DO (* diagonalization *)
kp1:=k+1;
it:=0;
REPEAT (* L100 *)
(* tests for splitting *)
cancellation:=TRUE;
ll:=kp1; (* k=mn,1,-1 *)
LOOP (* FOR l:=k TO 1 BY -1 DO *)
x := bnorm + ABS(E^[ll-1]);
IF (x = bnorm) THEN
cancellation:=FALSE;
EXIT;
END; (* IF *)
IF (ll = 0) THEN EXIT; END;
x := bnorm + ABS(D[ll-2]);
IF (x = bnorm) THEN
EXIT;
END; (* IF *)
DEC(ll);
END; (* LOOP FOR l *)
l:=ll;
IF cancellation THEN
(* cancellation of E[l] *)
c := 0.0;
s := 1.0;
ip1 := l; exit := FALSE;
WHILE (l <= kp1) AND NOT exit DO (* FOR i:=l TO k DO *)
i:=ip1-1;
f := s*E^[i];
E^[i] := c*E^[i];
x := bnorm + ABS(f);
IF (x = bnorm) THEN
exit := TRUE;
ELSE
g := D[i];
IF (ABS(g) < ABS(f)) THEN
c := g / f;
h := sqrt(1.0 + c*c);
s := -1.0 / h;
c := -c*s;
D[i] := f*h;
ELSE
s := f / g;
h := sqrt(1.0 + s*s);
c := 1.0 / h;
s := -s*c;
D[i] := g*h;
END; (* IF *)
(* update of the right-hand side *)
FOR j:=0 TO Nb-1 DO
y := B[l-2,j];
z := B[i,j];
B[l-2,j] := y*c + z*s;
B[i,j] := -y*s + z*c;
END;
INC(ip1);
END; (* IF *)
END; (* WHILE i *)
END; (* IF cancellation *)
(* test for convergence *)
z := D[k];
IF (l # kp1) THEN
(* shift from bottom principal matrix of order 2 *)
INC(it);
IF (it >= maxit) THEN
(* no convergence to a singular value after maxit iterations *)
iErr := kp1;
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
Errors.Pause(999);
RETURN;
END; (* IF *)
x := D[l-1];
y := D[kp1-2];
g := E^[kp1-2];
h := E^[k];
f := 0.5*(((g + z) / h)*((g - z) / y) + y/h - h/y);
IF (ABS(f) >= 1.0) THEN
f := x - (z/x)*z +
(h/x)*(y / (f*(1.0 + sqrt(1.0 + sqr(1.0/f)))) - h);
ELSE
f := x - (z/x)*z +
(h/x)*(y / (f + DSIGN(sqrt(1.0 + f*f),f)) - h);
END; (* IF *)
(* QR sweep *)
c := 1.0;
s := 1.0;
FOR i:=l-1 TO k-1 DO
ip1:=i+1;
g := E^[ip1];
y := D[ip1];
h := s*g;
g := c*g;
IF (ABS(f) < ABS(h)) THEN
c := f / h;
z := sqrt(1.0 + c*c);
s := 1.0 / z;
c := c*s;
E^[i] := h*z;
ELSE
s := h / f;
z := sqrt(1.0 + s*s);
c := 1.0 / z;
s := s*c;
E^[i] := f*z;
END; (* IF |f| < |h| *)
f := x*c + g*s;
g := -x*s + g*c;
h := y*s;
y := y*c;
(* update of matrix V *)
FOR j:=0 TO N-1 DO
x := A[i,j];
z := A[ip1,j];
A[i,j] := x*c + z*s;
A[ip1,j] := -x*s + z*c;
END; (* FOR j *)
z := DMAX(ABS(f),ABS(h));
z := z*sqrt(sqr(f/z) + sqr(h/z));
D[i] := z;
IF (z#0.0) THEN
c := f / z;
s := h / z;
END; (* IF *)
f := c*g + s*y;
x := -s*g + c*y;
(* update of the right-hand side *)
FOR j:=0 TO Nb-1 DO
y := B[i,j];
z := B[ip1,j];
B[i,j] := y*c + z*s;
B[ip1,j] := -y*s + z*c;
END; (* FOR *)
END; (* FOR i *)
E^[l-1] := 0.0;
E^[k] := f;
D[k] := x;
END; (* IF l # k *)
UNTIL (l = kp1); (* LOOP L100 *)
IF (z < 0.0) THEN
(* k-th singular value *)
D[k] := -z;
FOR j:=0 TO N-1 DO
A[k,j] := -A[k,j];
END;
END;
END; (* FOR k *)
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
END MinFit;
PROCEDURE BerLsg(VAR M,N : INTEGER;
IP : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL; (* [1..N,1..N] *)
VAR A0 : ARRAY OF ARRAY OF LONGREAL; (* [1..M,1..N] *)
VAR B,B0 : ARRAY OF ARRAY OF LONGREAL; (* [1..M,IP] *)
VAR X : ARRAY OF LONGREAL; (* [1..N] *)
VAR W : ARRAY OF LONGREAL; (* [1..N] *)
VAR R : ARRAY OF LONGREAL; (* [1..M] *)
VAR ifR : BOOLEAN);
VAR i,j : INTEGER;
BEGIN
IP:=IP - 1; (* Offene Felder *)
FOR i:=0 TO N-1 DO
(* Multiply B by reciprocals of singular values *)
IF (W[i] # 0.0) THEN
B[i,IP]:=B[i,IP] / W[i];
ELSE
B[i,IP]:=B[i,IP];
END;
END;
FOR j:=0 TO N-1 DO X[j]:=0.0; END;
FOR i:=0 TO N-1 DO
(* Multiply B by right singular vectors of A *)
FOR j:=0 TO N-1 DO
X[j]:=X[j] + A[i,j]*B[i,IP];
END;
END;
IF ifR THEN
(* Compute the residual out from the computed solution *)
(* vector X and the original B vector and A matrix *)
FOR i:=0 TO M-1 DO R[i] := - B0[i,IP]; END;
FOR i:=0 TO N-1 DO
FOR j:=0 TO M-1 DO
R[j]:=R[j] + A0[i,j]*X[i];
END;
END;
END; (* IF *)
END BerLsg;
(*=========================== LinPack Routinen =============================*)
PROCEDURE dSVDc(VAR A : ARRAY OF ARRAY OF LONGREAL;
m : INTEGER;
n : INTEGER;
VAR S : ARRAY OF LONGREAL;
VAR E : ARRAY OF LONGREAL;
VAR U,V : ARRAY OF ARRAY OF LONGREAL;
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 *>
<* IF (__INLINE__) THEN *>
(* Version mit "inline" code *)
<* END *>
(*----------------------------------------------------------------*)
<* IF (__INLINE__) THEN *>
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
<* END *>
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,ncu,nrt,nrtp1 : INTEGER;
scale,shift,sl,sm,smm1,sn,t,t1,test : LONGREAL;
wantu,wantv : BOOLEAN;
ztest : LONGREAL;
work : WorkSpace;
<* IF (__INLINE__) THEN *>
ii : INTEGER;
s,tmp,x1,y1 : LONGREAL;
<* END *>
<* IF (__BLASF77__) THEN *>
ione,num : INTEGER;
tmp : LONGREAL;
<* END *>
BEGIN
<* IF (__debug__) THEN *>
TIO.WrStr(tag); TIO.WrLn;
<* END *>
(* 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 := n;
END;
IF (jobu # 0) THEN
wantu := TRUE;
END;
IF ((job MOD 10) # 0) THEN
wantv := TRUE;
END;
<* IF (__debug__) THEN *>
TIO.WrStr("n = "); TIO.WrInt(n ,1); TIO.WrLn;
TIO.WrStr("m = "); TIO.WrInt(m ,1); TIO.WrLn;
TIO.WrStr("ncu = "); TIO.WrInt(ncu ,1); TIO.WrLn;
TIO.WrStr("jobu = "); TIO.WrInt(jobu,1); TIO.WrLn;
TIO.WrStr("job MOD 10 = "); TIO.WrInt(job MOD 10,1); TIO.WrLn;
TIO.WrStr("wantu = "); TIO.WrBool(wantu); TIO.WrLn;
TIO.WrStr("wantv = "); TIO.WrBool(wantv); TIO.WrLn;
<* 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);
<* IF (__debug__) THEN *>
TIO.WrStr("nct = "); TIO.WrInt(nct ,1); TIO.WrLn;
TIO.WrStr("nrt = "); TIO.WrInt(nrt ,1); TIO.WrLn;
TIO.WrStr("lu = "); TIO.WrInt(lu ,1); TIO.WrLn;
<* END *>
ALLOCATE(work,m*TSIZE(LONGREAL));
IF (work = NIL) THEN
Errors.Fehlerflag:="Kein Freispeicher vorhanden (SVDLib3.dSVDc)";
Errors.Fehler:=TRUE;
Errors.ErrOut(Errors.Fehlerflag);
info:=-1;
RETURN;
END;
<* IF (__BLASF77__) THEN *>
ione:=1; (* Fuer F77 Blas routinen *)
<* END *>
FOR l:=0 TO lu-1 DO
(* Compute the transformation for the L-th column and *)
(* place the L-th diagonal in S(L). *)
IF (l < nct) THEN
<* IF (__INLINE__) THEN *>
s:=0.0;
FOR ii:=l TO m-1 DO s:=s + sqr(A[l,ii]); END;
S[l]:=sqrt(s);
<* END *>
<* IF (__BLAS__) THEN *>
S[l] := dnrm2(m-l,A[l,l],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l;
S[l] := dnrm2(num,A[l,l],ione);
<* END *>
IF (S[l] # 0.0) THEN
IF (A[l,l] # 0.0) THEN
S[l] := dsign(S[l],A[l,l]);
END;
<* IF (__INLINE__) THEN *>
tmp:= 1.0 / S[l];
FOR ii:=l TO m-1 DO A[l,ii]:=A[l,ii]*tmp; END;
<* END *>
<* IF (__BLAS__) THEN *>
dscal(m-l,(1.0/S[l]),A[l,l],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l;
tmp:= 1.0 / S[l];
dscal(num,tmp,A[l,l],ione);
<* END *>
A[l,l]:=A[l,l] + 1.0;
END;
S[l] := -S[l];
END;
FOR j:=l+1 TO n-1 DO
(* Apply the transformation. *)
IF (l < nct) AND (S[l] # 0.0) THEN (* redundant ??? *)
<* IF (__INLINE__) THEN *>
t:=0.0;
FOR ii:=l TO m-1 DO t:=t + A[l,ii]*A[j,ii]; END;
t := - t / A[l,l];
FOR ii:=l TO m-1 DO
A[j,ii]:=A[j,ii] + t*A[l,ii];
END;
<* END *>
<* IF (__BLAS__) THEN *>
t := - ddot(m-l,A[l,l],1,A[j,l],1) / A[l,l];
daxpy(m-l,t,A[l,l],1,A[j,l],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l;
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] := 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-1 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 (__INLINE__) THEN *>
s:=0.0;
FOR ii:=l+1 TO n-1 DO s:=s + sqr(E[ii]); END;
E[l]:=sqrt(s);
<* END *>
<* IF (__BLAS__) THEN *>
E[l] := dnrm2(n-l,E[l+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-l;
E[l] := dnrm2(num,E[l+1],ione);
<* END *>
IF (E[l] # 0.0) THEN
IF (E[l+1] # 0.0) THEN
E[l] := dsign(E[l],E[l+1]);
END;
<* IF (__INLINE__) THEN *>
tmp:=1.0 / E[l];
FOR ii:=l+1 TO n-1 DO E[ii]:=E[ii]*tmp; END;
<* END *>
<* IF (__BLAS__) THEN *>
dscal(n-l-1,(1.0 / E[l]),E[l+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-l-1;
tmp:=1.0 / E[l];
dscal(num,tmp,E[l+1],ione);
<* END *>
E[l+1]:=E[l+1] + 1.0;
END;
E[l] := -E[l];
IF (l+1 < m) AND (E[l] # 0.0) THEN
(* Apply the transformation *)
FOR i:=l+1 TO m-1 DO
work^[i] := 0.0;
END;
FOR j:=l+1 TO n-1 DO
<* IF (__INLINE__) THEN *>
FOR ii:=l+1 TO m-1 DO
work^[ii]:=work^[ii] + E[j]*A[j,ii];
END;
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(m-l-1,E[j],A[j,l+1],1,work^[l+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l-1;
daxpy(num,E[j],A[j,l+1],ione,work^[l+1],ione);
<* END *>
END;
FOR j:=l+1 TO n-1 DO
<* IF (__INLINE__) THEN *>
tmp:= E[j] / E[l+1];
FOR ii:=l+1 TO m-1 DO
A[j,ii]:=A[j,ii] - tmp*work^[ii];
END;
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(m-l-1,-E[j]/E[l+1],work^[l+1],1,A[j,l+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:= m-l-1;
tmp:= -E[j] / E[l+1];
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-1 DO
V[l,i] := E[i];
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[nct] := A[nct,nct];
END;
IF (m < mn) THEN
S[mn-1] := 0.0;
END;
IF (nrtp1 < mn) THEN
E[nrt] := A[mn-1,nrt];
END;
E[mn-1] := 0.0;
IF (wantu) THEN (* If required, generate U. *)
FOR j:=nct TO ncu-1 DO
FOR i:=0 TO m-1 DO
U[j,i]:=0.0;
END;
U[j,j]:=1.0;
END;
(*
FOR l:=nct-1 TO 0 BY -1 DO
*)
FOR ll:=1 TO nct DO
l := nct - ll;
IF (S[l] # 0.0) THEN
FOR j:=l+1 TO ncu-1 DO
<* IF (__INLINE__) THEN *>
t:=0.0;
FOR ii:=l TO m-1 DO t:=t + U[l,ii]*U[j,ii]; END;
t:= - t / U[l,l];
FOR ii:=l TO m-1 DO U[j,ii]:=U[j,ii] + t*U[l,ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
t := - ddot(m-l,U[l,l],1,U[j,l],1) / U[l,l];
daxpy(m-l,t,U[l,l],1,U[j,l],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l;
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 (__INLINE__) THEN *>
FOR ii:=l TO m-1 DO U[l,ii]:= - U[l,ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
dscal(m-l,-1.0,U[l,l],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=m-l;
tmp:=-1.0;
dscal(num,tmp,U[l,l],ione);
<* END *>
U[l,l]:=U[l,l] + 1.0;
FOR i:=0 TO l-1 DO
U[l,i]:=0.0;
END;
ELSE
FOR i:=0 TO m-1 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-1 TO 0 BY -1 DO
*)
FOR ll:=1 TO n DO
l := n - ll;
lp1 := l + 1;
IF (l < nrt) AND (E[l] # 0.0) THEN
FOR j:=lp1 TO n-1 DO
<* IF (__INLINE__) THEN *>
t:=0.0;
FOR ii:=lp1 TO n-1 DO t:=t + V[l,ii]*V[j,ii]; END;
t := - t / V[l,lp1];
FOR ii:=lp1 TO n-1 DO
V[j,ii]:=V[j,ii] + t*V[l,ii];
END;
<* END *>
<* IF (__BLAS__) THEN *>
t := - ddot(n-l-1,V[l,lp1],1,V[j,lp1],1) / V[l,lp1];
daxpy(n-l-1,t,V[l,lp1],1,V[j,lp1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-l-1;
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:=0 TO n-1 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 offenen Feldern"); TIO.WrLn; TIO.WrLn;
TIO.WrStr("A +"); FMatEA.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:=0 TO n-1 DO
TIO.WrInt(l+1,3);
TIO.WrLngReal(E[l],12,4); TIO.WrLngReal(S[l],12,4);
TIO.WrLn;
END; TIO.WrLn;
mm:=imax(m,n)+1;
TIO.WrStr("U"); FMatEA.MATaus(FileSystem.Output,mm,mm,U,80,8,3);
TIO.WrStr("V"); FMatEA.MATaus(FileSystem.Output,mm,mm,V,80,8,3);
<* END *>
(* Main iteration loop for the singular values. *)
sn:=1.0; cs:=0.0; (* Wg. Compilerwarnung *)
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-1 TO l-1 BY -1 DO
*)
FOR kk:=l TO mm1 DO
k := mm1 - kk + l - 1;
t1 := S[k];
<* IF (__INLINE__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
<* IF (__BLAS__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
S[k] := t1;
IF (k+1 # l) THEN
f := -sn*E[k-1];
E[k-1] := cs*E[k-1];
END;
IF (wantv) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO (* rot V[k],V[mn-1,cs,sn *)
x1 := V[k ,ii];
y1 := V[mn-1,ii];
V[k ,ii] := cs*x1 + sn*y1;
V[mn-1,ii] := -sn*x1 + cs*y1;
END;
<* END *>
<* IF (__BLAS__) THEN *>
drot(n,V[k,0],1,V[mn-1,0],1,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drot(n,V[k,0],ione,V[mn-1,0],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-1 TO mn-1 DO
t1 := S[k];
<* IF (__INLINE__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
<* IF (__BLAS__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drotg(t1,f,cs,sn);
<* END *>
S[k] := t1;
f := - sn*E[k];
E[k] := cs*E[k];
IF (wantu) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO m-1 DO (* rot U[k],U[l-2,cs,sn *)
x1 := U[k ,ii];
y1 := U[l-2,ii];
U[k ,ii] := cs*x1 + sn*y1;
U[l-2,ii] := -sn*x1 + cs*y1;
END;
<* END *>
<* IF (__BLAS__) THEN *>
drot(m,U[k,0],1,U[l-2,0],1,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drot(m,U[k,0],ione,U[l-2,0],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-1 TO mm1-1 DO
<* IF (__INLINE__) THEN *>
drotg(f,g,cs,sn);
<* END *>
<* IF (__BLAS__) THEN *>
drotg(f,g,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drotg(f,g,cs,sn);
<* END *>
IF (k+1 # l) THEN
E[k-1] := f;
END;
f := cs*S[k] + sn*E[k];
E[k] := cs*E[k] - sn*S[k];
g := sn*S[k+1];
S[k+1] := cs*S[k+1];
IF (wantv) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO (* rot V[k],V[k+1,cs,sn *)
x1 := V[k ,ii];
y1 := V[k+1,ii];
V[k ,ii] := cs*x1 + sn*y1;
V[k+1,ii] := -sn*x1 + cs*y1;
END;
<* END *>
<* IF (__BLAS__) THEN *>
drot(n,V[k,0],1,V[k+1,0],1,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drot(n,V[k,0],ione,V[k+1,0],ione,cs,sn);
<* END *>
END;
<* IF (__INLINE__) THEN *>
drotg(f,g,cs,sn);
<* END *>
<* IF (__BLAS__) THEN *>
drotg(f,g,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drotg(f,g,cs,sn);
<* END *>
S[k] := f;
f := cs*E[k] + sn*S[k+1];
S[k+1] := -sn*E[k] + cs*S[k+1];
g := sn*E[k+1];
E[k+1] := cs*E[k+1];
IF (wantu) AND (k < m) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO m-1 DO (* rot U[k],U[k+1,cs,sn *)
x1 := U[k ,ii];
y1 := U[k+1,ii];
U[k ,ii] := cs*x1 + sn*y1;
U[k+1,ii] := -sn*x1 + cs*y1;
END;
<* END *>
<* IF (__BLAS__) THEN *>
drot(m,U[k,0],1,U[k+1,0],1,cs,sn);
<* END *>
<* IF (__BLASF77__) THEN *>
drot(m,U[k,0],ione,U[k+1,0],ione,cs,sn);
<* END *>
END;
END; (* FOR k *)
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 (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO V[l-1,ii]:= - V[l-1,ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
dscal(n,-1.0,V[l-1,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
tmp:=-1.0;
dscal(n,tmp,V[l-1,0],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 (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO (* swap V[l],V[l-1] *)
tmp :=V[l-1,ii];
V[l-1,ii]:=V[l ,ii];
V[l ,ii]:=tmp;
END;
<* END *>
<* IF (__BLAS__) THEN *>
dswap(n,V[l-1,0],1,V[l,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
dswap(n,V[l-1,0],ione,V[l,0],ione);
<* END *>
END;
IF (wantu) AND (l < m) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO m-1 DO (* swap U[l],U[l-1] *)
tmp :=U[l-1,ii];
U[l-1,ii]:=U[l ,ii];
U[l ,ii]:=tmp;
END;
<* END *>
<* IF (__BLAS__) THEN *>
dswap(m,U[l-1,0],1,U[l,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
dswap(m,U[l-1,0],ione,U[l,0],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 dSVDc;
BEGIN
Errors.WriteString(tag); Errors.WriteLn;
END SVDLib1.