|
a/SVDLib1.mod.m2pp |
|
b/SVDLib1.mod.m2pp |
|
... |
|
... |
| 70 |
(* 24.08.18, MRi: Die rechte Seite der Gleichungssysteme auf einheit- *)
|
70 |
(* 24.08.18, MRi: Die rechte Seite der Gleichungssysteme auf einheit- *)
|
| 71 |
(* liches Speicherformat [1..nB,1..M] gebracht *)
|
71 |
(* liches Speicherformat [1..nB,1..M] gebracht *)
|
| 72 |
(* 25.08.18, MRi: In BerLsg die Schleifenstrukturen optimiert *)
|
72 |
(* 25.08.18, MRi: In BerLsg die Schleifenstrukturen optimiert *)
|
| 73 |
(* 26.08.18, MRi: In BerLsg und MinFit komplett auf offene Felder umge- *)
|
73 |
(* 26.08.18, MRi: In BerLsg und MinFit komplett auf offene Felder umge- *)
|
| 74 |
(* stellt *)
|
74 |
(* stellt *)
|
|
|
75 |
(* 28.08.18, MRi: In MinFit die Indexberechnung bereinigt *)
|
| 75 |
(*------------------------------------------------------------------------*)
|
76 |
(*------------------------------------------------------------------------*)
|
| 76 |
(* Offene Punkte *)
|
77 |
(* Offene Punkte *)
|
| 77 |
(* *)
|
78 |
(* *)
|
| 78 |
(* - Durchsehen der Fehlermeldungen *)
|
79 |
(* - Durchsehen der Fehlermeldungen *)
|
| 79 |
(* - Durchsehen der Kommentare und Erklaerungen *)
|
80 |
(* - Durchsehen der Kommentare und Erklaerungen *)
|
|
... |
|
... |
| 99 |
(*------------------------------------------------------------------------*)
|
100 |
(*------------------------------------------------------------------------*)
|
| 100 |
(* Implementation : Michael Riedl *)
|
101 |
(* Implementation : Michael Riedl *)
|
| 101 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
102 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
| 102 |
(*------------------------------------------------------------------------*)
|
103 |
(*------------------------------------------------------------------------*)
|
| 103 |
|
104 |
|
| 104 |
(* $Id: SVDLib1.mod.m2pp,v 1.6 2018/08/26 14:25:59 mriedl Exp $ *)
|
105 |
(* $Id: SVDLib1.mod.m2pp,v 1.5 2018/08/28 12:30:50 mriedl Exp mriedl $ *)
|
| 105 |
|
106 |
|
| 106 |
FROM SYSTEM IMPORT TSIZE;
|
107 |
FROM SYSTEM IMPORT TSIZE;
|
| 107 |
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
|
108 |
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
|
| 108 |
IMPORT Errors;
|
109 |
IMPORT Errors;
|
| 109 |
FROM Errors IMPORT Fehlerflag;
|
110 |
FROM Errors IMPORT Fehlerflag;
|
|
... |
|
... |
| 1181 |
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
|
1182 |
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
|
| 1182 |
|
1183 |
|
| 1183 |
CONST maxit = 30;
|
1184 |
CONST maxit = 30;
|
| 1184 |
|
1185 |
|
| 1185 |
VAR i,j,k,kk,l,ll,mn,it,kache : INTEGER;
|
1186 |
VAR i,j,k,kk,l,ll,mn,it,kache : INTEGER;
|
|
|
1187 |
kp1,ip1 :INTEGER;
|
| 1186 |
c,s,x,y,z,f,g,h,bnorm,s16,r16 : LONGREAL;
|
1188 |
c,s,x,y,z,f,g,h,bnorm,s16,r16 : LONGREAL;
|
| 1187 |
cancellation,exit : BOOLEAN;
|
1189 |
cancellation,exit : BOOLEAN;
|
| 1188 |
E : POINTER TO ARRAY
|
1190 |
E : POINTER TO ARRAY
|
| 1189 |
[0..MAX(INTEGER)-1] OF LONGREAL;
|
1191 |
[0..MAX(INTEGER)-1] OF LONGREAL;
|
| 1190 |
BEGIN
|
1192 |
BEGIN
|
|
... |
|
... |
| 1196 |
kache := CacheLen(HIGH(A));
|
1198 |
kache := CacheLen(HIGH(A));
|
| 1197 |
kk := N - kache;
|
1199 |
kk := N - kache;
|
| 1198 |
|
1200 |
|
| 1199 |
(* reduction to upper-bidiagonal form by householder transformations *)
|
1201 |
(* reduction to upper-bidiagonal form by householder transformations *)
|
| 1200 |
|
1202 |
|
| 1201 |
k:=1;
|
1203 |
kp1:=1;
|
| 1202 |
LOOP (* FOR k:=1 TO mn DO *)
|
1204 |
LOOP (* FOR k:=1 TO mn DO *)
|
|
|
1205 |
k:=kp1-1;
|
| 1203 |
(* IF (k > mn) THEN EXIT; END; *)
|
1206 |
(* IF (k > mn) THEN EXIT; END; *)
|
| 1204 |
(* L1 norm of the subdiagonal part of column k *)
|
1207 |
(* L1 norm of the subdiagonal part of column k *)
|
| 1205 |
x:=0.0;
|
1208 |
x:=0.0;
|
| 1206 |
FOR i:=k+1 TO M DO
|
1209 |
FOR i:=kp1 TO M-1 DO
|
| 1207 |
x:=x + ABS(A[k-1,i-1]);
|
1210 |
x:=x + ABS(A[k,i]);
|
| 1208 |
END;
|
1211 |
END;
|
| 1209 |
IF (x # 0.0) THEN
|
1212 |
IF (x # 0.0) THEN
|
| 1210 |
x:=x + ABS(A[k-1,k-1]);
|
1213 |
x:=x + ABS(A[k,k]);
|
| 1211 |
(* scaled householder vector *)
|
1214 |
(* scaled householder vector *)
|
| 1212 |
DSC16(x,s16,r16);
|
1215 |
DSC16(x,s16,r16);
|
| 1213 |
s:=0.0;
|
1216 |
s:=0.0;
|
| 1214 |
FOR i:=k TO M DO
|
1217 |
FOR i:=k TO M-1 DO
|
| 1215 |
A[k-1,i-1] := r16*A[k-1,i-1];
|
1218 |
A[k,i] := r16*A[k,i];
|
| 1216 |
s:=s + sqr(A[k-1,i-1]);
|
1219 |
s:=s + sqr(A[k,i]);
|
| 1217 |
END;
|
1220 |
END;
|
| 1218 |
s := -DSIGN(sqrt(s),A[k-1,k-1]);
|
1221 |
s := -DSIGN(sqrt(s),A[k,k]);
|
| 1219 |
x := 1.0 / s;
|
1222 |
x := 1.0 / s;
|
| 1220 |
<* IF (__BLAS__) THEN *>
|
1223 |
<* IF (__BLAS__) THEN *>
|
| 1221 |
dscal(M-k+1,x,A[k-1,k-1],1);
|
1224 |
dscal(M-kp1+1,x,A[k,k],1);
|
| 1222 |
<* ELSE *>
|
1225 |
<* ELSE *>
|
| 1223 |
FOR i:=k TO M DO A[k-1,i-1] := x*A[k-1,i-1]; END;
|
1226 |
FOR i:=k TO M-1 DO
|
|
|
1227 |
A[k,i] := x*A[k,i];
|
|
|
1228 |
END;
|
| 1224 |
<* END *>
|
1229 |
<* END *>
|
| 1225 |
A[k-1,k-1] := A[k-1,k-1] - 1.0;
|
1230 |
A[k,k] := A[k,k] - 1.0;
|
| 1226 |
(* diagonal element *)
|
1231 |
(* diagonal element *)
|
| 1227 |
D[k-1] := s16*s;
|
1232 |
D[k] := s16*s;
|
| 1228 |
(* left transformation *)
|
1233 |
(* left transformation *)
|
| 1229 |
x := 1.0 / A[k-1,k-1];
|
1234 |
x := 1.0 / A[k,k];
|
| 1230 |
IF (k < kk) THEN
|
1235 |
IF (kp1 < kk) THEN
|
| 1231 |
(* path for small arrays *)
|
1236 |
(* path for small arrays *)
|
| 1232 |
FOR j:=k+1 TO N DO
|
1237 |
FOR j:=k+1 TO N-1 DO
|
| 1233 |
s:=0.0;
|
1238 |
s:=0.0;
|
| 1234 |
FOR i:=k TO M DO
|
1239 |
FOR i:=k TO M-1 DO
|
| 1235 |
s:=s + A[j-1,i-1]*A[k-1,i-1];
|
1240 |
s:=s + A[j,i]*A[k,i];
|
| 1236 |
END;
|
1241 |
END;
|
| 1237 |
s := s*x;
|
1242 |
s := s*x;
|
| 1238 |
FOR i:=k TO M DO
|
1243 |
FOR i:=k TO M-1 DO
|
| 1239 |
A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
|
1244 |
A[j,i]:=A[j,i] + s*A[k,i];
|
| 1240 |
END;
|
1245 |
END;
|
| 1241 |
E^[j-1] := A[j-1,k-1];
|
1246 |
E^[j] := A[j,k];
|
| 1242 |
END; (* FOR *)
|
1247 |
END; (* FOR *)
|
| 1243 |
ELSE
|
1248 |
ELSE
|
| 1244 |
(* path for large arrays *)
|
1249 |
(* path for large arrays *)
|
| 1245 |
FOR j:=k+1 TO N DO
|
1250 |
FOR j:=k+1 TO N-1 DO
|
| 1246 |
<* IF (__BLAS__) THEN *>
|
1251 |
<* IF (__BLAS__) THEN *>
|
| 1247 |
s := x*ddot(M-k+1,A[j-1,k-1],1,A[k-1,k-1],1);
|
1252 |
s := x*ddot(M-kp1+1,A[j,k],1,A[k,k],1);
|
| 1248 |
daxpy(M-k+1,s,A[k-1,k-1],1,A[j-1,k-1],1);
|
1253 |
daxpy(M-kp1+1,s,A[k,k],1,A[j,k],1);
|
| 1249 |
<* ELSE *>
|
1254 |
<* ELSE *>
|
| 1250 |
s:=0.0;
|
1255 |
s:=0.0;
|
| 1251 |
FOR i:=k TO M DO s:=s + A[j-1,i-1]*A[k-1,i-1]; END;
|
1256 |
FOR i:=k TO M-1 DO
|
|
|
1257 |
s:=s + A[j,i]*A[k,i];
|
|
|
1258 |
END;
|
| 1252 |
s := s*x;
|
1259 |
s := s*x;
|
| 1253 |
FOR i:=k TO M DO
|
1260 |
FOR i:=k TO M-1 DO
|
| 1254 |
A[j-1,i-1] := A[j-1,i-1] + s*A[k-1,i-1];
|
1261 |
A[j,i] := A[j,i] + s*A[k,i];
|
| 1255 |
END;
|
1262 |
END;
|
| 1256 |
<* END *>
|
1263 |
<* END *>
|
| 1257 |
E^[j-1] := A[j-1,k-1];
|
1264 |
E^[j] := A[j,k];
|
| 1258 |
END; (* FOR j *)
|
1265 |
END; (* FOR j *)
|
| 1259 |
|
1266 |
|
| 1260 |
END; (* IF k < kk *)
|
1267 |
END; (* IF k < kk *)
|
| 1261 |
(* transformation of the right-hand side *)
|
1268 |
(* transformation of the right-hand side *)
|
| 1262 |
FOR j:=1 TO Nb DO
|
1269 |
FOR j:=0 TO Nb-1 DO
|
| 1263 |
s:=0.0;
|
1270 |
s:=0.0;
|
| 1264 |
FOR i:=k TO M DO
|
1271 |
FOR i:=k TO M-1 DO
|
| 1265 |
s:=s + A[k-1,i-1]*B[i-1,j-1]; (* !!! *)
|
1272 |
s:=s + A[k,i]*B[i,j]; (* !!! *)
|
| 1266 |
END;
|
1273 |
END;
|
| 1267 |
s := s*x;
|
1274 |
s := s*x;
|
| 1268 |
FOR i:=k TO M DO
|
1275 |
FOR i:=k TO M-1 DO
|
| 1269 |
B[i-1,j-1]:=B[i-1,j-1] + s*A[k-1,i-1]; (* !!! *)
|
1276 |
B[i,j]:=B[i,j] + s*A[k,i]; (* !!! *)
|
| 1270 |
END;
|
1277 |
END;
|
| 1271 |
END; (* FOR j *)
|
1278 |
END; (* FOR j *)
|
| 1272 |
ELSE
|
1279 |
ELSE
|
| 1273 |
(* the transformation is bypassed if the norm is zero *)
|
1280 |
(* the transformation is bypassed if the norm is zero *)
|
| 1274 |
D[k-1] := A[k-1,k-1];
|
1281 |
D[k] := A[k,k];
|
| 1275 |
FOR j:=k+1 TO N DO
|
1282 |
FOR j:=k+1 TO N-1 DO
|
| 1276 |
E^[j-1] := A[j-1,k-1];
|
1283 |
E^[j] := A[j,k];
|
| 1277 |
END;
|
1284 |
END;
|
| 1278 |
END; (* IF x # 0 *)
|
1285 |
END; (* IF x # 0 *)
|
| 1279 |
IF (k = N) THEN
|
1286 |
IF (kp1 = N) THEN
|
| 1280 |
EXIT; (* LOOP k *)
|
1287 |
EXIT; (* LOOP k *)
|
| 1281 |
END;
|
1288 |
END;
|
| 1282 |
(* L1 norm of the supercodiagonal part of row k *)
|
1289 |
(* L1 norm of the supercodiagonal part of row k *)
|
| 1283 |
x:=0.0;
|
1290 |
x:=0.0;
|
| 1284 |
FOR j:=k+2 TO N DO
|
1291 |
FOR j:=k+2 TO N-1 DO
|
| 1285 |
x:=x + ABS(E^[j-1]);
|
1292 |
x:=x + ABS(E^[j]);
|
| 1286 |
END;
|
1293 |
END;
|
| 1287 |
IF (x # 0.0) THEN
|
1294 |
IF (x # 0.0) THEN
|
| 1288 |
(* scaled householder vector *)
|
1295 |
(* scaled householder vector *)
|
| 1289 |
x:=x + ABS(E^[k]);
|
1296 |
x:=x + ABS(E^[kp1]);
|
| 1290 |
DSC16(x,s16,r16);
|
1297 |
DSC16(x,s16,r16);
|
| 1291 |
s:=0.0;
|
1298 |
s:=0.0;
|
| 1292 |
FOR j:=k+1 TO N DO
|
1299 |
FOR j:=k+1 TO N-1 DO
|
| 1293 |
E^[j-1] := r16*E^[j-1];
|
1300 |
E^[j] := r16*E^[j];
|
| 1294 |
s:=s + sqr(E^[j-1]);
|
1301 |
s:=s + sqr(E^[j]);
|
| 1295 |
END;
|
1302 |
END;
|
| 1296 |
y := -DSIGN(sqrt(s),E^[k]);
|
1303 |
y := -DSIGN(sqrt(s),E^[kp1]);
|
| 1297 |
E^[k] := E^[k] - y;
|
1304 |
E^[kp1] := E^[kp1] - y;
|
| 1298 |
x := 1.0 / y;
|
1305 |
x := 1.0 / y;
|
| 1299 |
FOR j:=k+1 TO N DO
|
1306 |
FOR j:=k+1 TO N-1 DO
|
| 1300 |
E^[j-1] := x*E^[j-1];
|
1307 |
E^[j] := x*E^[j];
|
| 1301 |
A[k-1,j-1] := E^[j-1];
|
1308 |
A[k,j] := E^[j];
|
| 1302 |
END;
|
1309 |
END;
|
| 1303 |
x := 1.0 / E^[k];
|
1310 |
x := 1.0 / E^[kp1];
|
| 1304 |
(* right transformation *)
|
1311 |
(* right transformation *)
|
| 1305 |
FOR i:=k+1 TO M DO
|
1312 |
FOR i:=kp1 TO M-1 DO
|
| 1306 |
s:=0.0;
|
1313 |
s:=0.0;
|
| 1307 |
FOR j:=k+1 TO N DO
|
1314 |
FOR j:=k+1 TO N-1 DO
|
| 1308 |
s:=s + A[j-1,i-1]*E^[j-1];
|
1315 |
s:=s + A[j,i]*E^[j];
|
| 1309 |
END;
|
1316 |
END;
|
| 1310 |
s := s*x;
|
1317 |
s := s*x;
|
| 1311 |
FOR j:=k+1 TO N DO
|
1318 |
FOR j:=k+1 TO N-1 DO
|
| 1312 |
A[j-1,i-1]:=A[j-1,i-1] + E^[j-1]*s;
|
1319 |
A[j,i]:=A[j,i] + E^[j]*s;
|
| 1313 |
END;
|
1320 |
END;
|
| 1314 |
END; (* FOR *)
|
1321 |
END; (* FOR *)
|
| 1315 |
(* codiagonal element *)
|
1322 |
(* codiagonal element *)
|
| 1316 |
E^[k] := s16*y;
|
1323 |
E^[kp1] := s16*y;
|
| 1317 |
ELSE (* x = 0 *)
|
1324 |
ELSE (* x = 0 *)
|
| 1318 |
(* the transformation is bypassed if the norm is zero *)
|
1325 |
(* the transformation is bypassed if the norm is zero *)
|
| 1319 |
FOR i:=k+1 TO N DO
|
1326 |
FOR i:=kp1 TO N-1 DO
|
| 1320 |
A[k-1,i-1]:=0.0;
|
1327 |
A[k,i]:=0.0;
|
| 1321 |
END;
|
1328 |
END;
|
| 1322 |
END; (* IF x *)
|
1329 |
END; (* IF x *)
|
| 1323 |
INC(k);
|
1330 |
INC(kp1);
|
| 1324 |
END; (* LOOP k *)
|
1331 |
END; (* LOOP k *)
|
| 1325 |
|
1332 |
|
| 1326 |
l := mn;
|
1333 |
l := mn;
|
| 1327 |
IF (M < N) THEN
|
1334 |
IF (M < N) THEN
|
| 1328 |
INC(l);
|
1335 |
INC(l);
|
| 1329 |
END; (* IF *)
|
1336 |
END; (* IF *)
|
| 1330 |
|
1337 |
|
| 1331 |
(* right matrix of the bidiagonal decomposition *)
|
1338 |
(* right matrix of the bidiagonal decomposition *)
|
| 1332 |
FOR j:=l+1 TO N DO
|
1339 |
FOR j:=l TO N-1 DO
|
| 1333 |
FOR i:=1 TO N DO
|
1340 |
FOR i:=0 TO N-1 DO
|
| 1334 |
A[j-1,i-1]:=0.0;
|
1341 |
A[j,i]:=0.0;
|
| 1335 |
END;
|
1342 |
END;
|
| 1336 |
A[j-1,j-1]:=1.0;
|
1343 |
A[j,j]:=1.0;
|
| 1337 |
END; (* FOR j *)
|
1344 |
END; (* FOR j *)
|
| 1338 |
|
1345 |
|
| 1339 |
kk := N - kache;
|
1346 |
kk := N - kache;
|
| 1340 |
FOR k:=l TO 2 BY -1 DO
|
1347 |
FOR kp1:=l TO 2 BY -1 DO
|
|
|
1348 |
k:=kp1-1;
|
| 1341 |
FOR i:=1 TO k-1 DO
|
1349 |
FOR i:=0 TO k-1 DO
|
| 1342 |
A[k-1,i-1] := 0.0;
|
1350 |
A[k,i] := 0.0;
|
| 1343 |
END;
|
1351 |
END;
|
| 1344 |
FOR i:=k TO N DO
|
1352 |
FOR i:=k TO N-1 DO
|
| 1345 |
A[k-1,i-1] := A[k-1-1,i-1];
|
1353 |
A[k,i] := A[k-1,i];
|
| 1346 |
END;
|
1354 |
END;
|
| 1347 |
IF (A[k-1,k-1] # 0.0) THEN
|
1355 |
IF (A[k,k] # 0.0) THEN
|
| 1348 |
x := 1.0 / A[k-1,k-1];
|
1356 |
x := 1.0 / A[k,k];
|
| 1349 |
IF (k < kk) THEN
|
1357 |
IF (kp1 < kk) THEN
|
| 1350 |
(* path for small arrays *)
|
1358 |
(* path for small arrays *)
|
| 1351 |
FOR j:=k+1 TO N DO
|
1359 |
FOR j:=k+1 TO N-1 DO
|
| 1352 |
s:=0.0;
|
1360 |
s:=0.0;
|
| 1353 |
FOR i:=k TO N DO
|
1361 |
FOR i:=k TO N-1 DO
|
| 1354 |
s:=s + A[k-1,i-1]*A[j-1,i-1];
|
1362 |
s:=s + A[k,i]*A[j,i];
|
| 1355 |
END;
|
1363 |
END;
|
| 1356 |
s := s*x;
|
1364 |
s := s*x;
|
| 1357 |
FOR i:=k TO N DO
|
1365 |
FOR i:=k TO N-1 DO
|
| 1358 |
A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
|
1366 |
A[j,i]:=A[j,i] + s*A[k,i];
|
| 1359 |
END;
|
1367 |
END;
|
| 1360 |
END; (* FOR *)
|
1368 |
END; (* FOR *)
|
| 1361 |
ELSE
|
1369 |
ELSE
|
| 1362 |
(* path for large arrays *)
|
1370 |
(* path for large arrays *)
|
| 1363 |
FOR j:=k+1 TO N DO
|
1371 |
FOR j:=k+1 TO N-1 DO
|
| 1364 |
<* IF (__BLAS__) THEN *>
|
1372 |
<* IF (__BLAS__) THEN *>
|
| 1365 |
s := x*ddot(N-k+1,A[k-1,k-1],1,A[j-1,k-1],1);
|
1373 |
s := x*ddot(N-kp1+1,A[k,k],1,A[j,k],1);
|
| 1366 |
daxpy(N-k+1,s,A[k-1,k-1],1,A[j-1,k-1],1);
|
1374 |
daxpy(N-kp1+1,s,A[k,k],1,A[j,k],1);
|
| 1367 |
<* ELSE *>
|
1375 |
<* ELSE *>
|
| 1368 |
s := 0.0;
|
1376 |
s := 0.0;
|
| 1369 |
FOR i:=k TO N DO
|
1377 |
FOR i:=k TO N-1 DO
|
| 1370 |
s:=s + A[k-1,i-1]*A[j-1,i-1];
|
1378 |
s:=s + A[k,i]*A[j,i];
|
| 1371 |
END;
|
1379 |
END;
|
| 1372 |
s := s*x;
|
1380 |
s := s*x;
|
| 1373 |
FOR i:=k TO N DO
|
1381 |
FOR i:=k TO N-1 DO
|
| 1374 |
A[j-1,i-1]:=A[j-1,i-1] + s*A[k-1,i-1];
|
1382 |
A[j,i]:=A[j,i] + s*A[k,i];
|
| 1375 |
END;
|
1383 |
END;
|
| 1376 |
<* END *>
|
1384 |
<* END *>
|
| 1377 |
END; (* FOR j *)
|
1385 |
END; (* FOR j *)
|
| 1378 |
|
1386 |
|
| 1379 |
END; (* IF k < kk *)
|
1387 |
END; (* IF k < kk *)
|
| 1380 |
END; (* IF A[k-1,k-1] # 0 *)
|
1388 |
END; (* IF A[k-1,k-1] # 0 *)
|
| 1381 |
A[k-1,k-1]:=A[k-1,k-1] + 1.0;
|
1389 |
A[k,k]:=A[k,k] + 1.0;
|
| 1382 |
END; (* FOR k *)
|
1390 |
END; (* FOR k *)
|
| 1383 |
|
1391 |
|
| 1384 |
A[1-1,1-1]:=1.0;
|
1392 |
A[0,0]:=1.0;
|
| 1385 |
FOR i:=2 TO N DO
|
1393 |
FOR i:=1 TO N-1 DO
|
| 1386 |
A[1-1,i-1]:=0.0;
|
1394 |
A[0,i]:=0.0;
|
| 1387 |
END;
|
1395 |
END;
|
| 1388 |
|
1396 |
|
| 1389 |
(* annihilation of the last superdiagonal element when m < n *)
|
1397 |
(* annihilation of the last superdiagonal element when m < n *)
|
| 1390 |
IF (M < N) THEN
|
1398 |
IF (M < N) THEN
|
| 1391 |
c := 1.0;
|
1399 |
c := 1.0;
|
| 1392 |
s := -1.0;
|
1400 |
s := -1.0;
|
| 1393 |
k:=mn; exit:=FALSE;
|
1401 |
kp1:=mn; exit:=FALSE;
|
| 1394 |
WHILE (k > 1) AND NOT exit DO (* FOR k:=mn TO 1 BY -1 DO *)
|
1402 |
WHILE (kp1 > 1) AND NOT exit DO (* FOR k:=mn TO 1 BY -1 DO *)
|
| 1395 |
x := -s*E^[k];
|
1403 |
x := -s*E^[kp1];
|
| 1396 |
E^[k] := c*E^[k];
|
1404 |
E^[kp1] := c*E^[kp1];
|
| 1397 |
IF (x = 0.0) THEN
|
1405 |
IF (x = 0.0) THEN
|
| 1398 |
exit:=TRUE;
|
1406 |
exit:=TRUE;
|
| 1399 |
ELSE
|
1407 |
ELSE
|
| 1400 |
IF (ABS(D[k-1]) < ABS(x)) THEN
|
1408 |
IF (ABS(D[k]) < ABS(x)) THEN
|
| 1401 |
c := D[k-1] / x;
|
1409 |
c := D[k] / x;
|
| 1402 |
y := sqrt(1.0 + c*c);
|
1410 |
y := sqrt(1.0 + c*c);
|
| 1403 |
s := 1.0 / y;
|
1411 |
s := 1.0 / y;
|
| 1404 |
c := c*s;
|
1412 |
c := c*s;
|
| 1405 |
D[k-1] := y*x;
|
1413 |
D[k] := y*x;
|
| 1406 |
ELSE
|
1414 |
ELSE
|
| 1407 |
s := x / D[k-1];
|
1415 |
s := x / D[k];
|
| 1408 |
y := sqrt(1.0 + s*s);
|
1416 |
y := sqrt(1.0 + s*s);
|
| 1409 |
c := 1.0 / y;
|
1417 |
c := 1.0 / y;
|
| 1410 |
s := s*c;
|
1418 |
s := s*c;
|
| 1411 |
D[k-1] := y*D[k-1];
|
1419 |
D[k] := y*D[k];
|
| 1412 |
END; (* IF *)
|
1420 |
END; (* IF *)
|
| 1413 |
FOR i:=1 TO N DO
|
1421 |
FOR i:=0 TO N-1 DO
|
| 1414 |
x := A[k-1,i-1];
|
1422 |
x := A[k,i];
|
| 1415 |
y := A[l-1,i-1];
|
1423 |
y := A[l-1,i];
|
| 1416 |
A[k-1,i-1] := c*x + s*y;
|
1424 |
A[k,i] := c*x + s*y;
|
| 1417 |
A[l-1,i-1] := c*y - s*x;
|
1425 |
A[l-1,i] := c*y - s*x;
|
| 1418 |
END; (* FOR *)
|
1426 |
END; (* FOR *)
|
| 1419 |
END; (* IF *)
|
1427 |
END; (* IF *)
|
| 1420 |
DEC(k);
|
1428 |
DEC(kp1);
|
| 1421 |
END; (* WHILE k *)
|
1429 |
END; (* WHILE k *)
|
| 1422 |
END; (* IF M < N *)
|
1430 |
END; (* IF M < N *)
|
| 1423 |
|
1431 |
|
| 1424 |
FOR i:=mn+1 TO N DO
|
1432 |
FOR i:=mn TO N-1 DO
|
| 1425 |
FOR j:=1 TO Nb DO
|
1433 |
FOR j:=0 TO Nb-1 DO
|
| 1426 |
B[i-1,j-1]:=0.0;
|
1434 |
B[i,j]:=0.0;
|
| 1427 |
END;
|
1435 |
END;
|
| 1428 |
D[i-1]:=0.0;
|
1436 |
D[i]:=0.0;
|
| 1429 |
END;
|
1437 |
END;
|
| 1430 |
|
1438 |
|
| 1431 |
(* singular-value decomposition of the upper- *)
|
1439 |
(* singular-value decomposition of the upper- *)
|
| 1432 |
(* bidiagonal matrix of order min(m,n) *)
|
1440 |
(* bidiagonal matrix of order min(m,n) *)
|
| 1433 |
(* L1 norm of the bidiagonal matrix *)
|
1441 |
(* L1 norm of the bidiagonal matrix *)
|
| 1434 |
bnorm := 0.0;
|
1442 |
bnorm := 0.0;
|
| 1435 |
FOR i:=1 TO mn DO
|
1443 |
FOR i:=0 TO mn-1 DO
|
| 1436 |
bnorm := DMAX(ABS(D[i-1])+ABS(E^[i-1]),bnorm);
|
1444 |
bnorm := DMAX(ABS(D[i])+ABS(E^[i]),bnorm);
|
| 1437 |
END;
|
1445 |
END;
|
| 1438 |
|
1446 |
|
| 1439 |
FOR k:=mn TO 1 BY -1 DO (* diagonalization *)
|
1447 |
FOR k:=mn TO 1 BY -1 DO (* diagonalization *)
|
| 1440 |
|
1448 |
kp1:=k+1;
|
| 1441 |
it:=0;
|
1449 |
it:=0;
|
| 1442 |
REPEAT (* L100 *)
|
1450 |
REPEAT (* L100 *)
|
| 1443 |
|
|
|
| 1444 |
(* tests for splitting *)
|
1451 |
(* tests for splitting *)
|
| 1445 |
cancellation:=TRUE;
|
1452 |
cancellation:=TRUE;
|
| 1446 |
ll:=k; (* k=mn,1,-1 *)
|
1453 |
ll:=kp1; (* k=mn,1,-1 *)
|
| 1447 |
LOOP (* FOR l:=k TO 1 BY -1 DO *)
|
1454 |
LOOP (* FOR l:=k TO 1 BY -1 DO *)
|
| 1448 |
x := bnorm + ABS(E^[ll-1]);
|
1455 |
x := bnorm + ABS(E^[ll-1]);
|
| 1449 |
IF (x = bnorm) THEN
|
1456 |
IF (x = bnorm) THEN
|
| 1450 |
cancellation:=FALSE;
|
1457 |
cancellation:=FALSE;
|
| 1451 |
EXIT;
|
1458 |
EXIT;
|
|
... |
|
... |
| 1461 |
|
1468 |
|
| 1462 |
IF cancellation THEN
|
1469 |
IF cancellation THEN
|
| 1463 |
(* cancellation of E[l] *)
|
1470 |
(* cancellation of E[l] *)
|
| 1464 |
c := 0.0;
|
1471 |
c := 0.0;
|
| 1465 |
s := 1.0;
|
1472 |
s := 1.0;
|
| 1466 |
i := l; exit := FALSE;
|
1473 |
ip1 := l; exit := FALSE;
|
| 1467 |
WHILE (l <= k) AND NOT exit DO (* FOR i:=l TO k DO *)
|
1474 |
WHILE (l <= kp1) AND NOT exit DO (* FOR i:=l TO k DO *)
|
|
|
1475 |
i:=ip1-1;
|
| 1468 |
f := s*E^[i-1];
|
1476 |
f := s*E^[i];
|
| 1469 |
E^[i-1] := c*E^[i-1];
|
1477 |
E^[i] := c*E^[i];
|
| 1470 |
x := bnorm + ABS(f);
|
1478 |
x := bnorm + ABS(f);
|
| 1471 |
IF (x = bnorm) THEN
|
1479 |
IF (x = bnorm) THEN
|
| 1472 |
exit := TRUE;
|
1480 |
exit := TRUE;
|
| 1473 |
ELSE
|
1481 |
ELSE
|
| 1474 |
g := D[i-1];
|
1482 |
g := D[i];
|
| 1475 |
IF (ABS(g) < ABS(f)) THEN
|
1483 |
IF (ABS(g) < ABS(f)) THEN
|
| 1476 |
c := g / f;
|
1484 |
c := g / f;
|
| 1477 |
h := sqrt(1.0 + c*c);
|
1485 |
h := sqrt(1.0 + c*c);
|
| 1478 |
s := -1.0 / h;
|
1486 |
s := -1.0 / h;
|
| 1479 |
c := -c*s;
|
1487 |
c := -c*s;
|
| 1480 |
D[i-1] := f*h;
|
1488 |
D[i] := f*h;
|
| 1481 |
ELSE
|
1489 |
ELSE
|
| 1482 |
s := f / g;
|
1490 |
s := f / g;
|
| 1483 |
h := sqrt(1.0 + s*s);
|
1491 |
h := sqrt(1.0 + s*s);
|
| 1484 |
c := 1.0 / h;
|
1492 |
c := 1.0 / h;
|
| 1485 |
s := -s*c;
|
1493 |
s := -s*c;
|
| 1486 |
D[i-1] := g*h;
|
1494 |
D[i] := g*h;
|
| 1487 |
END; (* IF *)
|
1495 |
END; (* IF *)
|
| 1488 |
(* update of the right-hand side *)
|
1496 |
(* update of the right-hand side *)
|
| 1489 |
FOR j:=1 TO Nb DO
|
1497 |
FOR j:=0 TO Nb-1 DO
|
| 1490 |
y := B[l-2,j-1];
|
1498 |
y := B[l-2,j];
|
| 1491 |
z := B[i-1,j-1];
|
1499 |
z := B[i,j];
|
| 1492 |
B[l-2,j-1] := y*c + z*s;
|
1500 |
B[l-2,j] := y*c + z*s;
|
| 1493 |
B[i-1,j-1] := -y*s + z*c;
|
1501 |
B[i,j] := -y*s + z*c;
|
| 1494 |
END;
|
1502 |
END;
|
| 1495 |
INC(i);
|
1503 |
INC(ip1);
|
| 1496 |
END; (* IF *)
|
1504 |
END; (* IF *)
|
| 1497 |
END; (* FOR i *)
|
1505 |
END; (* WHILE i *)
|
| 1498 |
END; (* IF cancellation *)
|
1506 |
END; (* IF cancellation *)
|
| 1499 |
|
1507 |
|
| 1500 |
(* test for convergence *)
|
1508 |
(* test for convergence *)
|
| 1501 |
z := D[k-1];
|
1509 |
z := D[k];
|
| 1502 |
IF (l # k) THEN
|
1510 |
IF (l # kp1) THEN
|
| 1503 |
(* shift from bottom principal matrix of order 2 *)
|
1511 |
(* shift from bottom principal matrix of order 2 *)
|
| 1504 |
INC(it);
|
1512 |
INC(it);
|
| 1505 |
IF (it >= maxit) THEN
|
1513 |
IF (it >= maxit) THEN
|
| 1506 |
(* no convergence to a singular value after maxit iterations *)
|
1514 |
(* no convergence to a singular value after maxit iterations *)
|
| 1507 |
iErr := k;
|
1515 |
iErr := kp1;
|
| 1508 |
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
|
1516 |
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
|
| 1509 |
Errors.Pause(999);
|
1517 |
Errors.Pause(999);
|
| 1510 |
RETURN;
|
1518 |
RETURN;
|
| 1511 |
END; (* IF *)
|
1519 |
END; (* IF *)
|
| 1512 |
x := D[l-1];
|
1520 |
x := D[l-1];
|
| 1513 |
y := D[k-2];
|
1521 |
y := D[kp1-2];
|
| 1514 |
g := E^[k-2];
|
1522 |
g := E^[kp1-2];
|
| 1515 |
h := E^[k-1];
|
1523 |
h := E^[k];
|
| 1516 |
f := 0.5*(((g + z) / h)*((g - z) / y) + y/h - h/y);
|
1524 |
f := 0.5*(((g + z) / h)*((g - z) / y) + y/h - h/y);
|
| 1517 |
IF (ABS(f) >= 1.0) THEN
|
1525 |
IF (ABS(f) >= 1.0) THEN
|
| 1518 |
f := x - (z/x)*z +
|
1526 |
f := x - (z/x)*z +
|
| 1519 |
(h/x)*(y / (f*(1.0 + sqrt(1.0 + sqr(1.0/f)))) - h);
|
1527 |
(h/x)*(y / (f*(1.0 + sqrt(1.0 + sqr(1.0/f)))) - h);
|
| 1520 |
ELSE
|
1528 |
ELSE
|
|
... |
|
... |
| 1522 |
(h/x)*(y / (f + DSIGN(sqrt(1.0 + f*f),f)) - h);
|
1530 |
(h/x)*(y / (f + DSIGN(sqrt(1.0 + f*f),f)) - h);
|
| 1523 |
END; (* IF *)
|
1531 |
END; (* IF *)
|
| 1524 |
(* QR sweep *)
|
1532 |
(* QR sweep *)
|
| 1525 |
c := 1.0;
|
1533 |
c := 1.0;
|
| 1526 |
s := 1.0;
|
1534 |
s := 1.0;
|
| 1527 |
FOR i:=l TO k-1 DO
|
1535 |
FOR i:=l-1 TO k-1 DO
|
|
|
1536 |
ip1:=i+1;
|
| 1528 |
g := E^[i];
|
1537 |
g := E^[ip1];
|
| 1529 |
y := D[i];
|
1538 |
y := D[ip1];
|
| 1530 |
h := s*g;
|
1539 |
h := s*g;
|
| 1531 |
g := c*g;
|
1540 |
g := c*g;
|
| 1532 |
IF (ABS(f) < ABS(h)) THEN
|
1541 |
IF (ABS(f) < ABS(h)) THEN
|
| 1533 |
c := f / h;
|
1542 |
c := f / h;
|
| 1534 |
z := sqrt(1.0 + c*c);
|
1543 |
z := sqrt(1.0 + c*c);
|
| 1535 |
s := 1.0 / z;
|
1544 |
s := 1.0 / z;
|
| 1536 |
c := c*s;
|
1545 |
c := c*s;
|
| 1537 |
E^[i-1] := h*z;
|
1546 |
E^[i] := h*z;
|
| 1538 |
ELSE
|
1547 |
ELSE
|
| 1539 |
s := h / f;
|
1548 |
s := h / f;
|
| 1540 |
z := sqrt(1.0 + s*s);
|
1549 |
z := sqrt(1.0 + s*s);
|
| 1541 |
c := 1.0 / z;
|
1550 |
c := 1.0 / z;
|
| 1542 |
s := s*c;
|
1551 |
s := s*c;
|
| 1543 |
E^[i-1] := f*z;
|
1552 |
E^[i] := f*z;
|
| 1544 |
END; (* IF |f| < |h| *)
|
1553 |
END; (* IF |f| < |h| *)
|
| 1545 |
f := x*c + g*s;
|
1554 |
f := x*c + g*s;
|
| 1546 |
g := -x*s + g*c;
|
1555 |
g := -x*s + g*c;
|
| 1547 |
h := y*s;
|
1556 |
h := y*s;
|
| 1548 |
y := y*c;
|
1557 |
y := y*c;
|
| 1549 |
(* update of matrix V *)
|
1558 |
(* update of matrix V *)
|
| 1550 |
FOR j:=1 TO N DO
|
1559 |
FOR j:=0 TO N-1 DO
|
| 1551 |
x := A[i-1,j-1];
|
1560 |
x := A[i,j];
|
| 1552 |
z := A[i+1-1,j-1];
|
1561 |
z := A[ip1,j];
|
| 1553 |
A[i-1,j-1] := x*c + z*s;
|
1562 |
A[i,j] := x*c + z*s;
|
| 1554 |
A[i+1-1,j-1] := -x*s + z*c;
|
1563 |
A[ip1,j] := -x*s + z*c;
|
| 1555 |
END; (* FOR j *)
|
1564 |
END; (* FOR j *)
|
| 1556 |
z := DMAX(ABS(f),ABS(h));
|
1565 |
z := DMAX(ABS(f),ABS(h));
|
| 1557 |
z := z*sqrt(sqr(f/z) + sqr(h/z));
|
1566 |
z := z*sqrt(sqr(f/z) + sqr(h/z));
|
| 1558 |
D[i-1] := z;
|
1567 |
D[i] := z;
|
| 1559 |
IF (z#0.0) THEN
|
1568 |
IF (z#0.0) THEN
|
| 1560 |
c := f / z;
|
1569 |
c := f / z;
|
| 1561 |
s := h / z;
|
1570 |
s := h / z;
|
| 1562 |
END; (* IF *)
|
1571 |
END; (* IF *)
|
| 1563 |
f := c*g + s*y;
|
1572 |
f := c*g + s*y;
|
| 1564 |
x := -s*g + c*y;
|
1573 |
x := -s*g + c*y;
|
| 1565 |
(* update of the right-hand side *)
|
1574 |
(* update of the right-hand side *)
|
| 1566 |
FOR j:=1 TO Nb DO
|
1575 |
FOR j:=0 TO Nb-1 DO
|
| 1567 |
y := B[i-1,j-1];
|
1576 |
y := B[i,j];
|
| 1568 |
z := B[i+1-1,j-1];
|
1577 |
z := B[ip1,j];
|
| 1569 |
B[i-1,j-1] := y*c + z*s;
|
1578 |
B[i,j] := y*c + z*s;
|
| 1570 |
B[i+1-1,j-1] := -y*s + z*c;
|
1579 |
B[ip1,j] := -y*s + z*c;
|
| 1571 |
END; (* FOR *)
|
1580 |
END; (* FOR *)
|
| 1572 |
END; (* FOR *)
|
1581 |
END; (* FOR i *)
|
| 1573 |
E^[l-1] := 0.0;
|
1582 |
E^[l-1] := 0.0;
|
| 1574 |
E^[k-1] := f;
|
1583 |
E^[k] := f;
|
| 1575 |
D[k-1] := x;
|
1584 |
D[k] := x;
|
| 1576 |
END; (* IF l # k *)
|
1585 |
END; (* IF l # k *)
|
| 1577 |
UNTIL (l = k); (* LOOP L100 *)
|
1586 |
UNTIL (l = kp1); (* LOOP L100 *)
|
| 1578 |
|
1587 |
|
| 1579 |
IF (z < 0.0) THEN
|
1588 |
IF (z < 0.0) THEN
|
| 1580 |
(* k-th singular value *)
|
1589 |
(* k-th singular value *)
|
| 1581 |
D[k-1] := -z;
|
1590 |
D[k] := -z;
|
| 1582 |
FOR j:=1 TO N DO
|
1591 |
FOR j:=0 TO N-1 DO
|
| 1583 |
A[k-1,j-1] := -A[k-1,j-1];
|
1592 |
A[k,j] := -A[k,j];
|
| 1584 |
END;
|
1593 |
END;
|
| 1585 |
END;
|
1594 |
END;
|
| 1586 |
|
1595 |
|
| 1587 |
END; (* FOR k *)
|
1596 |
END; (* FOR k *)
|
| 1588 |
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
|
1597 |
DEALLOCATE(E,MAX0(M,N)*TSIZE(LONGREAL));
|