IMPLEMENTATION MODULE SVDLib2;
(*------------------------------------------------------------------------*)
(* Modul zur Singul"arwertzerlegung. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 10.01.97, MRi: Durchsicht *)
(* 29.09.17, MRi: Herausloesen aus LinLib *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Prozedur OrderSVD auf schnelleren Sortieralgorithmus umstellen *)
(*------------------------------------------------------------------------*)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SVDLib2.mod,v 1.1 2018/01/16 09:20:42 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE,ADR;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Deklera IMPORT PMATRIX,PVEKTOR;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MaxCard,Pythag;
FROM MatLib IMPORT AbsSumVek,SumVek;
PROCEDURE SVD(VAR A : PMATRIX; (* Dynamische Matrix A[1..n,1..m] *)
m,n : INTEGER;
VAR W : ARRAY OF LONGREAL;
VAR V : PMATRIX);
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 SVD;
PROCEDURE SVDSolv(VAR U : PMATRIX; (* m,n *)
Utr : CHAR;
VAR W : ARRAY OF LONGREAL; (* n *)
VAR V : PMATRIX; (* n,n *)
VAR X : ARRAY OF LONGREAL; (* n *)
VAR C : ARRAY OF LONGREAL;
m,n : CARDINAL);
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 (SVDLib2.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 SVDSolv;
PROCEDURE OrderSVD(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 OrderSVD;
END SVDLib2.