--- a/SVDLib1.mod.m2pp
+++ b/SVDLib1.mod.m2pp
@@ -73,6 +73,9 @@
(* 26.08.18, MRi: In BerLsg und MinFit komplett auf offene Felder umge- *)
(* stellt *)
(* 28.08.18, MRi: In MinFit die Indexberechnung bereinigt *)
+ (* 27.09.18, MRi: In PowSVD Initialisierung "restlicher" Singulaerwerte *)
+ (* und -vektoren mit 0 eingefuegt falls der Rang der *)
+ (* Matrix < N ist. *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
@@ -102,7 +105,7 @@
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
- (* $Id: SVDLib1.mod.m2pp,v 1.5 2018/08/28 12:30:50 mriedl Exp mriedl $ *)
+ (* $Id: SVDLib1.mod.m2pp,v 1.7 2019/02/01 22:24:03 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
@@ -209,7 +212,7 @@
tol = 10.0*MachEps;
MAXINT = MAX(INTEGER);
- VAR i,j,k,nm : INTEGER;
+ VAR i,j,k,kk,nm : INTEGER;
counter,climit : INTEGER;
s,norm,d : LONGREAL;
dPRi,Pj,Ri : LONGREAL;
@@ -294,12 +297,20 @@
IF (counter > climit) THEN
iFehl:=2;
- Fehlerflag:="climit exceeded (PowSVD)";
+ Fehlerflag:="climit exceeded (SVDLib1.PowSVD)";
END;
IF ((s / (sv[0] + tol)) <= tol) THEN (* step 7 *)
rankdef:=TRUE;
iFehl:=iFehl + 1;
- Fehlerflag:="Matrix rank deficite (PowSVD)";
+ Fehlerflag:="Matrix rank deficite (SVDLib1.PowSVD)";
+(********)
+ FOR kk:=k TO N-1 DO (* Added 27.09.18 *)
+ sv[kk] := 0.0;
+ FOR i:=0 TO M-1 DO U[kk,i]:=0.0; END;
+ FOR i:=0 TO N-1 DO V[kk,i]:=0.0; END;
+ END;
+ DEC(k);
+(********)
ELSE
sv[k] := s; (* step 8 - convergence achieved *)
FOR i:=0 TO M-1 DO U[k,i]:=Q^[i]; END;
@@ -635,7 +646,7 @@
END;
FOR i:=0 TO nRHS-1 DO rss[i]:=0.0; END;
- tol2 := FLOAT(N*N)*eps*eps;
+ tol2 := LFLOAT(N*N)*eps*eps;
<* IF (__debug__) THEN *>
TIO.WrLn; TIO.WrStr("Matrix (A|B)"); TIO.WrLn; TIO.WrLn;
<* END *>
@@ -691,7 +702,7 @@
slimit := N DIV 4; IF (slimit < 6) THEN slimit := 6;END;
sweep := 0;
- e2 := 10.0*FLOAT(N)*eps*eps;
+ e2 := 10.0*LFLOAT(N)*eps*eps;
EstRowRank := N;
REPEAT
count := 0;
@@ -719,7 +730,7 @@
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;
+ IF (p < 0.0) THEN s := -s;END;
c := p / (vt*s);
RotNSub(0,j,k,tcol,c,s);
INC(count);
@@ -1141,7 +1152,7 @@
PROCEDURE MinFit( M,N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
- Nb : INTEGER;
+ nB : INTEGER;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
@@ -1266,7 +1277,7 @@
END; (* IF k < kk *)
(* transformation of the right-hand side *)
- FOR j:=0 TO Nb-1 DO
+ 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]; (* !!! *)
@@ -1430,7 +1441,7 @@
END; (* IF M < N *)
FOR i:=mn TO N-1 DO
- FOR j:=0 TO Nb-1 DO
+ FOR j:=0 TO nB-1 DO
B[i,j]:=0.0;
END;
D[i]:=0.0;
@@ -1494,7 +1505,7 @@
D[i] := g*h;
END; (* IF *)
(* update of the right-hand side *)
- FOR j:=0 TO Nb-1 DO
+ FOR j:=0 TO nB-1 DO
y := B[l-2,j];
z := B[i,j];
B[l-2,j] := y*c + z*s;
@@ -1572,7 +1583,7 @@
f := c*g + s*y;
x := -s*g + c*y;
(* update of the right-hand side *)
- FOR j:=0 TO Nb-1 DO
+ FOR j:=0 TO nB-1 DO
y := B[i,j];
z := B[ip1,j];
B[i,j] := y*c + z*s;