--- a/SVDLib1.mod.m2pp
+++ b/SVDLib1.mod.m2pp
@@ -72,6 +72,7 @@
(* 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 *)
(* *)
@@ -101,7 +102,7 @@
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
- (* $Id: SVDLib1.mod.m2pp,v 1.6 2018/08/26 14:25:59 mriedl Exp $ *)
+ (* $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;
@@ -1183,6 +1184,7 @@
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
@@ -1198,129 +1200,134 @@
(* reduction to upper-bidiagonal form by householder transformations *)
- k:=1;
+ 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:=k+1 TO M DO
- x:=x + ABS(A[k-1,i-1]);
+ FOR i:=kp1 TO M-1 DO
+ x:=x + ABS(A[k,i]);
END;
IF (x # 0.0) THEN
- x:=x + ABS(A[k-1,k-1]);
+ x:=x + ABS(A[k,k]);
(* scaled householder vector *)
DSC16(x,s16,r16);
s:=0.0;
- FOR i:=k TO M DO
- A[k-1,i-1] := r16*A[k-1,i-1];
- s:=s + sqr(A[k-1,i-1]);
- END;
- s := -DSIGN(sqrt(s),A[k-1,k-1]);
+ 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-k+1,x,A[k-1,k-1],1);
+ dscal(M-kp1+1,x,A[k,k],1);
<* ELSE *>
- FOR i:=k TO M DO A[k-1,i-1] := x*A[k-1,i-1]; END;
+ FOR i:=k TO M-1 DO
+ A[k,i] := x*A[k,i];
+ END;
<* END *>
- A[k-1,k-1] := A[k-1,k-1] - 1.0;
+ A[k,k] := A[k,k] - 1.0;
(* diagonal element *)
- D[k-1] := s16*s;
+ D[k] := s16*s;
(* left transformation *)
- x := 1.0 / A[k-1,k-1];
- IF (k < kk) THEN
+ x := 1.0 / A[k,k];
+ IF (kp1 < kk) THEN
(* path for small arrays *)
- FOR j:=k+1 TO N DO
+ FOR j:=k+1 TO N-1 DO
s:=0.0;
- FOR i:=k TO M DO
- s:=s + A[j-1,i-1]*A[k-1,i-1];
+ FOR i:=k TO M-1 DO
+ s:=s + A[j,i]*A[k,i];
END;
s := s*x;
- FOR i:=k TO M DO
- A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
- END;
- E^[j-1] := A[j-1,k-1];
+ 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 DO
+ FOR j:=k+1 TO N-1 DO
<* IF (__BLAS__) THEN *>
- s := x*ddot(M-k+1,A[j-1,k-1],1,A[k-1,k-1],1);
- daxpy(M-k+1,s,A[k-1,k-1],1,A[j-1,k-1],1);
+ 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 DO s:=s + A[j-1,i-1]*A[k-1,i-1]; END;
+ FOR i:=k TO M-1 DO
+ s:=s + A[j,i]*A[k,i];
+ END;
s := s*x;
- FOR i:=k TO M DO
- A[j-1,i-1] := A[j-1,i-1] + s*A[k-1,i-1];
+ FOR i:=k TO M-1 DO
+ A[j,i] := A[j,i] + s*A[k,i];
END;
<* END *>
- E^[j-1] := A[j-1,k-1];
+ E^[j] := A[j,k];
END; (* FOR j *)
END; (* IF k < kk *)
(* transformation of the right-hand side *)
- FOR j:=1 TO Nb DO
+ FOR j:=0 TO Nb-1 DO
s:=0.0;
- FOR i:=k TO M DO
- s:=s + A[k-1,i-1]*B[i-1,j-1]; (* !!! *)
+ FOR i:=k TO M-1 DO
+ s:=s + A[k,i]*B[i,j]; (* !!! *)
END;
s := s*x;
- FOR i:=k TO M DO
- B[i-1,j-1]:=B[i-1,j-1] + s*A[k-1,i-1]; (* !!! *)
+ 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-1] := A[k-1,k-1];
- FOR j:=k+1 TO N DO
- E^[j-1] := A[j-1,k-1];
+ D[k] := A[k,k];
+ FOR j:=k+1 TO N-1 DO
+ E^[j] := A[j,k];
END;
END; (* IF x # 0 *)
- IF (k = N) THEN
+ 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 DO
- x:=x + ABS(E^[j-1]);
+ 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^[k]);
+ x:=x + ABS(E^[kp1]);
DSC16(x,s16,r16);
s:=0.0;
- FOR j:=k+1 TO N DO
- E^[j-1] := r16*E^[j-1];
- s:=s + sqr(E^[j-1]);
- END;
- y := -DSIGN(sqrt(s),E^[k]);
- E^[k] := E^[k] - y;
+ 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 DO
- E^[j-1] := x*E^[j-1];
- A[k-1,j-1] := E^[j-1];
- END;
- x := 1.0 / E^[k];
+ 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:=k+1 TO M DO
+ FOR i:=kp1 TO M-1 DO
s:=0.0;
- FOR j:=k+1 TO N DO
- s:=s + A[j-1,i-1]*E^[j-1];
+ 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 DO
- A[j-1,i-1]:=A[j-1,i-1] + E^[j-1]*s;
+ FOR j:=k+1 TO N-1 DO
+ A[j,i]:=A[j,i] + E^[j]*s;
END;
END; (* FOR *)
(* codiagonal element *)
- E^[k] := s16*y;
+ E^[kp1] := s16*y;
ELSE (* x = 0 *)
(* the transformation is bypassed if the norm is zero *)
- FOR i:=k+1 TO N DO
- A[k-1,i-1]:=0.0;
+ FOR i:=kp1 TO N-1 DO
+ A[k,i]:=0.0;
END;
END; (* IF x *)
- INC(k);
+ INC(kp1);
END; (* LOOP k *)
l := mn;
@@ -1329,121 +1336,121 @@
END; (* IF *)
(* right matrix of the bidiagonal decomposition *)
- FOR j:=l+1 TO N DO
- FOR i:=1 TO N DO
- A[j-1,i-1]:=0.0;
- END;
- A[j-1,j-1]:=1.0;
+ 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 k:=l TO 2 BY -1 DO
- FOR i:=1 TO k-1 DO
- A[k-1,i-1] := 0.0;
- END;
- FOR i:=k TO N DO
- A[k-1,i-1] := A[k-1-1,i-1];
- END;
- IF (A[k-1,k-1] # 0.0) THEN
- x := 1.0 / A[k-1,k-1];
- IF (k < kk) THEN
+ 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 DO
+ FOR j:=k+1 TO N-1 DO
s:=0.0;
- FOR i:=k TO N DO
- s:=s + A[k-1,i-1]*A[j-1,i-1];
+ FOR i:=k TO N-1 DO
+ s:=s + A[k,i]*A[j,i];
END;
s := s*x;
- FOR i:=k TO N DO
- A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
+ 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 DO
+ FOR j:=k+1 TO N-1 DO
<* IF (__BLAS__) THEN *>
- s := x*ddot(N-k+1,A[k-1,k-1],1,A[j-1,k-1],1);
- daxpy(N-k+1,s,A[k-1,k-1],1,A[j-1,k-1],1);
+ 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 DO
- s:=s + A[k-1,i-1]*A[j-1,i-1];
+ FOR i:=k TO N-1 DO
+ s:=s + A[k,i]*A[j,i];
END;
s := s*x;
- FOR i:=k TO N DO
- A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
+ 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-1,k-1]:=A[k-1,k-1] + 1.0;
+ A[k,k]:=A[k,k] + 1.0;
END; (* FOR k *)
- A[1-1,1-1]:=1.0;
- FOR i:=2 TO N DO
- A[1-1,i-1]:=0.0;
+ 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;
- k:=mn; exit:=FALSE;
- WHILE (k > 1) AND NOT exit DO (* FOR k:=mn TO 1 BY -1 DO *)
- x := -s*E^[k];
- E^[k] := c*E^[k];
+ 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-1]) < ABS(x)) THEN
- c := D[k-1] / x;
+ 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-1] := y*x;
+ D[k] := y*x;
ELSE
- s := x / D[k-1];
+ s := x / D[k];
y := sqrt(1.0 + s*s);
c := 1.0 / y;
s := s*c;
- D[k-1] := y*D[k-1];
+ D[k] := y*D[k];
END; (* IF *)
- FOR i:=1 TO N DO
- x := A[k-1,i-1];
- y := A[l-1,i-1];
- A[k-1,i-1] := c*x + s*y;
- A[l-1,i-1] := c*y - s*x;
+ 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(k);
+ DEC(kp1);
END; (* WHILE k *)
END; (* IF M < N *)
- FOR i:=mn+1 TO N DO
- FOR j:=1 TO Nb DO
- B[i-1,j-1]:=0.0;
- END;
- D[i-1]:=0.0;
+ 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:=1 TO mn DO
- bnorm := DMAX(ABS(D[i-1])+ABS(E^[i-1]),bnorm);
+ 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:=k; (* k=mn,1,-1 *)
+ 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
@@ -1463,56 +1470,57 @@
(* cancellation of E[l] *)
c := 0.0;
s := 1.0;
- i := l; exit := FALSE;
- WHILE (l <= k) AND NOT exit DO (* FOR i:=l TO k DO *)
- f := s*E^[i-1];
- E^[i-1] := c*E^[i-1];
+ 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-1];
+ 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-1] := f*h;
+ D[i] := f*h;
ELSE
s := f / g;
h := sqrt(1.0 + s*s);
c := 1.0 / h;
s := -s*c;
- D[i-1] := g*h;
+ D[i] := g*h;
END; (* IF *)
(* update of the right-hand side *)
- FOR j:=1 TO Nb DO
- y := B[l-2,j-1];
- z := B[i-1,j-1];
- B[l-2,j-1] := y*c + z*s;
- B[i-1,j-1] := -y*s + z*c;
+ 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(i);
+ INC(ip1);
END; (* IF *)
- END; (* FOR i *)
+ END; (* WHILE i *)
END; (* IF cancellation *)
(* test for convergence *)
- z := D[k-1];
- IF (l # k) THEN
+ 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 := k;
+ iErr := kp1;
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
Errors.Pause(999);
RETURN;
END; (* IF *)
x := D[l-1];
- y := D[k-2];
- g := E^[k-2];
- h := E^[k-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 +
@@ -1524,9 +1532,10 @@
(* QR sweep *)
c := 1.0;
s := 1.0;
- FOR i:=l TO k-1 DO
- g := E^[i];
- y := D[i];
+ 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
@@ -1534,28 +1543,28 @@
z := sqrt(1.0 + c*c);
s := 1.0 / z;
c := c*s;
- E^[i-1] := h*z;
+ E^[i] := h*z;
ELSE
s := h / f;
z := sqrt(1.0 + s*s);
c := 1.0 / z;
s := s*c;
- E^[i-1] := f*z;
+ 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:=1 TO N DO
- x := A[i-1,j-1];
- z := A[i+1-1,j-1];
- A[i-1,j-1] := x*c + z*s;
- A[i+1-1,j-1] := -x*s + z*c;
+ 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-1] := z;
+ D[i] := z;
IF (z#0.0) THEN
c := f / z;
s := h / z;
@@ -1563,24 +1572,24 @@
f := c*g + s*y;
x := -s*g + c*y;
(* update of the right-hand side *)
- FOR j:=1 TO Nb DO
- y := B[i-1,j-1];
- z := B[i+1-1,j-1];
- B[i-1,j-1] := y*c + z*s;
- B[i+1-1,j-1] := -y*s + z*c;
+ 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 *)
+ END; (* FOR i *)
E^[l-1] := 0.0;
- E^[k-1] := f;
- D[k-1] := x;
+ E^[k] := f;
+ D[k] := x;
END; (* IF l # k *)
- UNTIL (l = k); (* LOOP L100 *)
+ UNTIL (l = kp1); (* LOOP L100 *)
IF (z < 0.0) THEN
(* k-th singular value *)
- D[k-1] := -z;
- FOR j:=1 TO N DO
- A[k-1,j-1] := -A[k-1,j-1];
+ D[k] := -z;
+ FOR j:=0 TO N-1 DO
+ A[k,j] := -A[k,j];
END;
END;