--- a/SVDLib1.mod.m2pp
+++ b/SVDLib1.mod.m2pp
@@ -70,6 +70,8 @@
(* 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 *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
@@ -99,11 +101,10 @@
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
- (* $Id: SVDLib1.mod.m2pp,v 1.5 2018/08/25 15:04:22 mriedl Exp mriedl $ *)
+ (* $Id: SVDLib1.mod.m2pp,v 1.6 2018/08/26 14:25:59 mriedl Exp $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
-FROM Deklera IMPORT VEKTOR,MATRIX;
IMPORT Errors;
FROM Errors IMPORT Fehlerflag;
FROM LongMath IMPORT sqrt;
@@ -1140,7 +1141,7 @@
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
Nb : INTEGER;
- VAR B : MATRIX;
+ VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
@@ -1194,8 +1195,6 @@
mn := MIN0(M,N);
kache := CacheLen(HIGH(A));
kk := N - kache;
-
-kk:=2;
(* reduction to upper-bidiagonal form by householder transformations *)
@@ -1263,11 +1262,11 @@
FOR j:=1 TO Nb DO
s:=0.0;
FOR i:=k TO M DO
- s:=s + A[k-1,i-1]*B[i,j]; (* !!! *)
+ s:=s + A[k-1,i-1]*B[i-1,j-1]; (* !!! *)
END;
s := s*x;
FOR i:=k TO M DO
- B[i,j]:=B[i,j] + s*A[k-1,i-1]; (* !!! *)
+ B[i-1,j-1]:=B[i-1,j-1] + s*A[k-1,i-1]; (* !!! *)
END;
END; (* FOR j *)
ELSE
@@ -1424,7 +1423,7 @@
FOR i:=mn+1 TO N DO
FOR j:=1 TO Nb DO
- B[i,j]:=0.0;
+ B[i-1,j-1]:=0.0;
END;
D[i-1]:=0.0;
END;
@@ -1488,10 +1487,10 @@
END; (* IF *)
(* update of the right-hand side *)
FOR j:=1 TO Nb DO
- y := B[l-1,j];
- z := B[i,j];
- B[l-1,j] := y*c + z*s;
- B[i,j] := -y*s + z*c;
+ 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;
END;
INC(i);
END; (* IF *)
@@ -1565,10 +1564,10 @@
x := -s*g + c*y;
(* update of the right-hand side *)
FOR j:=1 TO Nb DO
- y := B[i,j];
- z := B[i+1,j];
- B[i,j] := y*c + z*s;
- B[i+1,j] := -y*s + z*c;
+ 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;
END; (* FOR *)
END; (* FOR *)
E^[l-1] := 0.0;
@@ -1591,18 +1590,18 @@
PROCEDURE BerLsg(VAR M,N : INTEGER;
IP : INTEGER;
- VAR A : MATRIX; (* [1..N,1..N] *)
- VAR A0 : MATRIX; (* [1..M,1..N] *)
- VAR B,B0 : MATRIX; (* [1..M,IP] *)
- VAR X : VEKTOR; (* [1..N] *)
- VAR W : VEKTOR; (* [1..N] *)
- VAR R : VEKTOR; (* [1..M] *)
+ 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;
- Xi : LONGREAL;
BEGIN
- FOR i:=1 TO N DO
+ 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];
@@ -1611,10 +1610,10 @@
END;
END;
- FOR j:=1 TO N DO X[j]:=0.0; END;
- FOR i:=1 TO N DO
+ 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:=1 TO N DO
+ FOR j:=0 TO N-1 DO
X[j]:=X[j] + A[i,j]*B[i,IP];
END;
END;
@@ -1622,9 +1621,9 @@
IF ifR THEN
(* Compute the residual out from the computed solution *)
(* vector X and the original B vector and A matrix *)
- FOR i:=1 TO M DO R[i] := - B0[i,IP]; END;
- FOR i:=1 TO N DO
- FOR j:=1 TO M DO
+ 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;