|
a |
|
b/SVDLib3.mod |
|
|
1 |
IMPLEMENTATION MODULE SVDLib3;
|
|
|
2 |
|
|
|
3 |
(*------------------------------------------------------------------------*)
|
|
|
4 |
(* Module provides routines for Takagi factorisation and a complex value *)
|
|
|
5 |
(* SVD routine. *)
|
|
|
6 |
(* *)
|
|
|
7 |
(* Takagi is a M2 translation of Diag library subroutines Takagi. *)
|
|
|
8 |
(* The Diag libray is developed and maintained by Thomas Hahn, *)
|
|
|
9 |
(* Max Planck Institut fuer Physik http:/www.feyarts.de/diag *)
|
|
|
10 |
(*------------------------------------------------------------------------*)
|
|
|
11 |
(* Letzte Veraenderung *)
|
|
|
12 |
(* *)
|
|
|
13 |
(* 09.08.11, ThH: last modified of Fortran orgiginal *)
|
|
|
14 |
(* 14.09.16, MRi: Erstellen der ersten uebersetzbaren Takagi Version *)
|
|
|
15 |
(* 15.09.16, MRi: Fehlerkorrektur - Takagi scheint zu funktionieren *)
|
|
|
16 |
(* 28.09.16, MRi: Umstellen der Sortierung in SVD auf externe Routine *)
|
|
|
17 |
(* 20.06.18, MRi: Erstellen der erstern uebersetzbaren Version von *)
|
|
|
18 |
(* zSVDC und der benoetigten BLAS Routinen *)
|
|
|
19 |
(* 21.06.18, MRi: Korrekturen in dznrm2 und zdotc *)
|
|
|
20 |
(* Testmatrix wird mit zSVDc korrekt berechnet *)
|
|
|
21 |
(*------------------------------------------------------------------------*)
|
|
|
22 |
(* Testroutinen fuer zSVDc in TstSVDLib3a, fuer Takagi in TstSVDLib3b *)
|
|
|
23 |
(*------------------------------------------------------------------------*)
|
|
|
24 |
(* Offene Punkte *)
|
|
|
25 |
(* *)
|
|
|
26 |
(* - Weitere Verbesserung der Indizierung in zSVDc ([i-1] Problem) *)
|
|
|
27 |
(*------------------------------------------------------------------------*)
|
|
|
28 |
(* Implementation : Michael Riedl *)
|
|
|
29 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
|
|
30 |
(*------------------------------------------------------------------------*)
|
|
|
31 |
|
|
|
32 |
(* $Id: SVDLib3.mod,v 1.3 2018/07/28 07:06:02 mriedl Exp mriedl $ *)
|
|
|
33 |
|
|
|
34 |
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
|
|
|
35 |
FROM Deklera IMPORT SIZELONGCMPLX;
|
|
|
36 |
IMPORT Errors;
|
|
|
37 |
FROM LongMath IMPORT sqrt;
|
|
|
38 |
IMPORT LongComplexMath;
|
|
|
39 |
FROM LMathLib IMPORT MachEps,CardPot,sign2;
|
|
|
40 |
FROM LibDBlas IMPORT drotg,zswap,zdotc,dznrm2,zscal,zaxpy,zdrot;
|
|
|
41 |
|
|
|
42 |
CONST CABS = LongComplexMath.abs;
|
|
|
43 |
conj = LongComplexMath.conj;
|
|
|
44 |
scalarMult = LongComplexMath.scalarMult;
|
|
|
45 |
|
|
|
46 |
zero = LongComplexMath.zero;
|
|
|
47 |
one = LongComplexMath.one;
|
|
|
48 |
|
|
|
49 |
PROCEDURE zSVDc(VAR X : ARRAY OF ARRAY OF LONGCOMPLEX;
|
|
|
50 |
N,P : INTEGER;
|
|
|
51 |
VAR S : ARRAY OF LONGCOMPLEX;
|
|
|
52 |
VAR E : ARRAY OF LONGCOMPLEX;
|
|
|
53 |
VAR U : ARRAY OF ARRAY OF LONGCOMPLEX;
|
|
|
54 |
VAR V : ARRAY OF ARRAY OF LONGCOMPLEX;
|
|
|
55 |
VAR Work : ARRAY OF LONGCOMPLEX;
|
|
|
56 |
Job : INTEGER;
|
|
|
57 |
VAR Info : INTEGER);
|
|
|
58 |
|
|
|
59 |
(* X = U * s * conj(tran(V)) *)
|
|
|
60 |
(* X[P,N], V[P,P], U[N,N] *)
|
|
|
61 |
(* Benutzt BLAS zaxpy,zdotc,zscal,zswap,dznrm2,drotg *)
|
|
|
62 |
|
|
|
63 |
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
|
|
|
64 |
|
|
|
65 |
PROCEDURE MAX0(a,b : CARDINAL) : CARDINAL;
|
|
|
66 |
|
|
|
67 |
BEGIN
|
|
|
68 |
IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
|
|
|
69 |
END MAX0;
|
|
|
70 |
|
|
|
71 |
PROCEDURE MIN0(a,b : CARDINAL) : CARDINAL;
|
|
|
72 |
|
|
|
73 |
BEGIN
|
|
|
74 |
IF (a >= b) THEN RETURN b; ELSE RETURN a; END;
|
|
|
75 |
END MIN0;
|
|
|
76 |
|
|
|
77 |
PROCEDURE DMAX1(a,b : LONGREAL) : LONGREAL;
|
|
|
78 |
|
|
|
79 |
BEGIN
|
|
|
80 |
IF (a >= b) THEN RETURN a; ELSE RETURN b; END;
|
|
|
81 |
END DMAX1;
|
|
|
82 |
|
|
|
83 |
PROCEDURE CABS1(x : LONGCOMPLEX) : LONGREAL;
|
|
|
84 |
|
|
|
85 |
BEGIN
|
|
|
86 |
RETURN ABS(RE(x)) + ABS(IM(x));
|
|
|
87 |
END CABS1;
|
|
|
88 |
|
|
|
89 |
PROCEDURE CSIGN(a,b : LONGCOMPLEX) : LONGCOMPLEX;
|
|
|
90 |
|
|
|
91 |
BEGIN
|
|
|
92 |
RETURN scalarMult(1.0/CABS(b), scalarMult(CABS(a),b));
|
|
|
93 |
END CSIGN;
|
|
|
94 |
CONST maxit = 30;
|
|
|
95 |
|
|
|
96 |
VAR i,j,k,l,m,kk,ll,mm : INTEGER;
|
|
|
97 |
iter,jobu,kase,nctp1,nrtp1 : INTEGER;
|
|
|
98 |
lls,lm1,lp1,ls,lu,mm1,mp1,nct,ncu,nrt : INTEGER;
|
|
|
99 |
t,r : LONGCOMPLEX;
|
|
|
100 |
b,c,cs,sn,el,emm1,f,g : LONGREAL;
|
|
|
101 |
scale,shift,test,ztest : LONGREAL;
|
|
|
102 |
sl,sm,smm1,t1 : LONGREAL;
|
|
|
103 |
wantu,wantv : BOOLEAN;
|
|
|
104 |
BEGIN
|
|
|
105 |
(* determine what is to be computed. *)
|
|
|
106 |
|
|
|
107 |
wantu := FALSE;
|
|
|
108 |
wantv := FALSE;
|
|
|
109 |
jobu := (Job MOD 100) DIV 10;
|
|
|
110 |
ncu := N;
|
|
|
111 |
IF (jobu > 1) THEN
|
|
|
112 |
IF (N <= P) THEN ncu:=N; ELSE ncu:=P; END;
|
|
|
113 |
END;
|
|
|
114 |
IF (jobu # 0) THEN
|
|
|
115 |
wantu := TRUE;
|
|
|
116 |
END;
|
|
|
117 |
IF ((Job MOD 10) # 0) THEN
|
|
|
118 |
wantv := TRUE;
|
|
|
119 |
END;
|
|
|
120 |
|
|
|
121 |
(* reduce x to bidiagonal form, storing the diagonal *)
|
|
|
122 |
(* elements in s and the super-diagonal elements in e. *)
|
|
|
123 |
|
|
|
124 |
Info := 0;
|
|
|
125 |
nct := MIN0(N-1,P);
|
|
|
126 |
nrt := MAX0(0,MIN0(P-2,N));
|
|
|
127 |
lu := MAX0(nct,nrt);
|
|
|
128 |
IF (lu >= 1) THEN
|
|
|
129 |
FOR l:=1 TO lu DO
|
|
|
130 |
lp1 := l + 1;
|
|
|
131 |
IF (l <= nct) THEN
|
|
|
132 |
(* compute the transformation for the l-th column *)
|
|
|
133 |
(* and place the l-th diagonal in s(l). *)
|
|
|
134 |
S[l-1] := CMPLX(dznrm2(N-l+1,X[l-1,l-1],1),0.0);
|
|
|
135 |
IF (CABS1(S[l-1]) # 0.0) THEN
|
|
|
136 |
IF (CABS1(X[l-1,l-1]) # 0.0) THEN
|
|
|
137 |
S[l-1] := CSIGN(S[l-1],X[l-1,l-1]);
|
|
|
138 |
END;
|
|
|
139 |
zscal(N-l+1,1.0/S[l-1],X[l-1,l-1],1);
|
|
|
140 |
X[l-1,l-1] := one + X[l-1,l-1];
|
|
|
141 |
END; (* IF *)
|
|
|
142 |
S[l-1] := -S[l-1];
|
|
|
143 |
END; (* IF *)
|
|
|
144 |
IF (P >= lp1) THEN
|
|
|
145 |
FOR j:=lp1 TO P DO
|
|
|
146 |
IF (l <= nct) THEN
|
|
|
147 |
IF (CABS1(S[l-1]) # 0.0) THEN
|
|
|
148 |
(* apply the transformation. *)
|
|
|
149 |
t := -zdotc(N-l+1,X[l-1,l-1],1,X[j-1,l-1],1)/X[l-1,l-1];
|
|
|
150 |
zaxpy(N-l+1,t,X[l-1,l-1],1,X[j-1,l-1],1);
|
|
|
151 |
END; (* IF *)
|
|
|
152 |
END; (* IF *)
|
|
|
153 |
(* place the l-th row of x into e for the *)
|
|
|
154 |
(* subsequent calculation of the row transformation. *)
|
|
|
155 |
E[j-1] := conj(X[j-1,l-1]);
|
|
|
156 |
END; (* FOR *)
|
|
|
157 |
END; (* IF *)
|
|
|
158 |
IF NOT ((NOT wantu) OR (l > nct)) THEN
|
|
|
159 |
(* place the transformation in u for *)
|
|
|
160 |
(* subsequent back multiplication. *)
|
|
|
161 |
FOR i:=l TO N DO
|
|
|
162 |
U[l-1,i-1] := X[l-1,i-1];
|
|
|
163 |
END; (* FOR *)
|
|
|
164 |
END; (* IF *)
|
|
|
165 |
IF (l <= nrt) THEN
|
|
|
166 |
(* compute the l-th row transformation and *)
|
|
|
167 |
(* place the l-th super-diagonal in e(l). *)
|
|
|
168 |
E[l-1] := CMPLX(dznrm2(P-l,E[lp1-1],1),0.0);
|
|
|
169 |
IF (CABS1(E[l-1]) # 0.0) THEN
|
|
|
170 |
IF (CABS1(E[lp1-1]) # 0.0) THEN
|
|
|
171 |
E[l-1] := CSIGN(E[l-1],E[lp1-1]);
|
|
|
172 |
END; (* IF *)
|
|
|
173 |
zscal(P-l,1.0/E[l-1],E[lp1-1],1);
|
|
|
174 |
E[lp1-1] := one + E[lp1-1];
|
|
|
175 |
END; (* IF *)
|
|
|
176 |
E[l-1] := -conj(E[l-1]);
|
|
|
177 |
IF (lp1 <= N) AND (CABS1(E[l-1]) # 0.0) THEN
|
|
|
178 |
(* apply the transformation. *)
|
|
|
179 |
FOR i:=lp1 TO N DO Work[i-1] := zero; END;
|
|
|
180 |
FOR j:=lp1 TO P DO
|
|
|
181 |
zaxpy(N-l,E[j-1],X[j-1,lp1-1],1,Work[lp1-1],1);
|
|
|
182 |
END;
|
|
|
183 |
FOR j:=lp1 TO P DO
|
|
|
184 |
zaxpy(N-l,conj(-E[j-1]/E[lp1-1]),Work[lp1-1],1,X[j-1,lp1-1],1);
|
|
|
185 |
END; (* FOR *)
|
|
|
186 |
END; (* IF *)
|
|
|
187 |
IF wantv THEN
|
|
|
188 |
(* place the transformation in V for *)
|
|
|
189 |
(* subsequent back multiplication. *)
|
|
|
190 |
FOR i:=lp1 TO P DO
|
|
|
191 |
V[l-1,i-1] := E[i-1];
|
|
|
192 |
END;
|
|
|
193 |
END; (* IF *)
|
|
|
194 |
END; (* IF *)
|
|
|
195 |
END; (* FOR *)
|
|
|
196 |
END; (* IF *)
|
|
|
197 |
|
|
|
198 |
(* set up the final bidiagonal matrix or order m *)
|
|
|
199 |
|
|
|
200 |
m := MIN0(P,N+1);
|
|
|
201 |
nctp1 := nct + 1;
|
|
|
202 |
nrtp1 := nrt + 1;
|
|
|
203 |
IF (nct < P) THEN
|
|
|
204 |
S[nctp1-1] := X[nctp1-1,nctp1-1];
|
|
|
205 |
END; (* IF *)
|
|
|
206 |
IF (N < m) THEN
|
|
|
207 |
S[m-1] := zero;
|
|
|
208 |
END; (* IF *)
|
|
|
209 |
IF (nrtp1 < m) THEN
|
|
|
210 |
E[nrtp1-1] := X[m-1,nrtp1-1];
|
|
|
211 |
END; (* IF *)
|
|
|
212 |
E[m-1] := zero;
|
|
|
213 |
|
|
|
214 |
IF (wantu) THEN (* if required, generate U *)
|
|
|
215 |
IF (ncu >= nctp1) THEN
|
|
|
216 |
FOR j:=nctp1 TO ncu DO
|
|
|
217 |
FOR i:=1 TO N DO U[j-1,i-1] := zero; END;
|
|
|
218 |
U[j-1,j-1] := one;
|
|
|
219 |
END; (* FOR *)
|
|
|
220 |
END; (* IF *)
|
|
|
221 |
IF (nct >= 1) THEN
|
|
|
222 |
FOR ll:=1 TO nct DO
|
|
|
223 |
l := nct - ll + 1;
|
|
|
224 |
IF (CABS1(S[l-1]) = 0.0) THEN
|
|
|
225 |
FOR i:=1 TO N DO
|
|
|
226 |
U[l-1,i-1] := zero;
|
|
|
227 |
END; (* FOR *)
|
|
|
228 |
U[l-1,l-1] := one;
|
|
|
229 |
ELSE
|
|
|
230 |
lp1 := l + 1;
|
|
|
231 |
IF (ncu >= lp1) THEN
|
|
|
232 |
FOR j:=lp1 TO ncu DO
|
|
|
233 |
t := -zdotc(N-l+1,U[l-1,l-1],1,U[j-1,l-1],1) /
|
|
|
234 |
U[l-1,l-1];
|
|
|
235 |
zaxpy(N-l+1,t,U[l-1,l-1],1,U[j-1,l-1],1);
|
|
|
236 |
END; (* FOR *)
|
|
|
237 |
END; (* IF *)
|
|
|
238 |
zscal(N-l+1,-one,U[l-1,l-1],1);
|
|
|
239 |
U[l-1,l-1] := one + U[l-1,l-1];
|
|
|
240 |
lm1 := l - 1;
|
|
|
241 |
IF (lm1 >= 1) THEN
|
|
|
242 |
FOR i := 1 TO lm1 DO U[l-1,i-1] := zero; END;
|
|
|
243 |
END; (* IF *)
|
|
|
244 |
END; (* IF *)
|
|
|
245 |
END; (* FOR *)
|
|
|
246 |
END; (* IF *)
|
|
|
247 |
END; (* IF *)
|
|
|
248 |
|
|
|
249 |
IF (wantv) THEN (* if it is required, generate V *)
|
|
|
250 |
FOR ll:=1 TO P DO
|
|
|
251 |
l := P - ll + 1;
|
|
|
252 |
lp1 := l + 1;
|
|
|
253 |
IF (l <= nrt) THEN
|
|
|
254 |
IF (CABS1(E[l-1]) # 0.0) THEN
|
|
|
255 |
FOR j:=lp1 TO P DO
|
|
|
256 |
t := -zdotc(P-l,V[l-1,lp1-1],1,V[j-1,lp1-1],1)/V[l-1,lp1-1];
|
|
|
257 |
zaxpy(P-l,t,V[l-1,lp1-1],1,V[j-1,lp1-1],1);
|
|
|
258 |
END;
|
|
|
259 |
END; (* IF *)
|
|
|
260 |
END; (* IF *)
|
|
|
261 |
FOR i:=0 TO P-1 DO V[l-1,i] := zero; END;
|
|
|
262 |
V[l-1,l-1] := one;
|
|
|
263 |
END; (* FOR *)
|
|
|
264 |
END; (* IF *)
|
|
|
265 |
|
|
|
266 |
|
|
|
267 |
(* transform s and e so that they are double precision. *)
|
|
|
268 |
|
|
|
269 |
FOR i:=1 TO m DO
|
|
|
270 |
IF (CABS1(S[i-1]) # 0.0) THEN
|
|
|
271 |
t := CMPLX(CABS(S[i-1]),0.0);
|
|
|
272 |
r := S[i-1] / t;
|
|
|
273 |
S[i-1] := t;
|
|
|
274 |
IF (i < m) THEN
|
|
|
275 |
E[i-1] := E[i-1] / r;
|
|
|
276 |
END;
|
|
|
277 |
IF (wantu) THEN
|
|
|
278 |
zscal(N,r,U[i-1,1 -1],1);
|
|
|
279 |
END; (* IF *)
|
|
|
280 |
END; (* IF *)
|
|
|
281 |
|
|
|
282 |
IF (i # m) THEN
|
|
|
283 |
IF (CABS1(E[i-1]) # 0.0) THEN
|
|
|
284 |
t := CMPLX(CABS(E[i-1]),0.0);
|
|
|
285 |
r := t / E[i-1];
|
|
|
286 |
E[i-1] := t;
|
|
|
287 |
S[i+1-1] := S[i+1-1]*r;
|
|
|
288 |
IF (wantv) THEN
|
|
|
289 |
zscal(P,r,V[i+1-1,1-1],1);
|
|
|
290 |
END;
|
|
|
291 |
END; (* IF *)
|
|
|
292 |
END; (* IF *)
|
|
|
293 |
END; (* FOR *)
|
|
|
294 |
|
|
|
295 |
(* main iteration loop for the singular values. *)
|
|
|
296 |
|
|
|
297 |
mm := m;
|
|
|
298 |
iter := 0;
|
|
|
299 |
|
|
|
300 |
(* quit if all the singular values have been found. *)
|
|
|
301 |
|
|
|
302 |
(* exit ??? *)
|
|
|
303 |
|
|
|
304 |
WHILE (m # 0) DO
|
|
|
305 |
|
|
|
306 |
(* if too many iterations have been performed, set *)
|
|
|
307 |
(* flag and return. *)
|
|
|
308 |
|
|
|
309 |
IF (iter < maxit) THEN
|
|
|
310 |
|
|
|
311 |
(* This section of the program inspects for *)
|
|
|
312 |
(* negligible elements in the s and e arrays. On *)
|
|
|
313 |
(* completion the variables kase and l are set as follows. *)
|
|
|
314 |
(* *)
|
|
|
315 |
(* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m *)
|
|
|
316 |
(* kase = 2 if s(l) is negligible and l.lt.m *)
|
|
|
317 |
(* kase = 3 if e(l-1) is negligible, l.lt.m, and *)
|
|
|
318 |
(* s(l), ..., s(m) are not negligible (qr step). *)
|
|
|
319 |
(* kase = 4 if e(m-1) is negligible (convergence). *)
|
|
|
320 |
|
|
|
321 |
ll:=1;
|
|
|
322 |
LOOP
|
|
|
323 |
l := m - ll;
|
|
|
324 |
IF (l = 0) THEN EXIT; END;
|
|
|
325 |
test := CABS(S[l-1]) + CABS(S[l+1-1]);
|
|
|
326 |
ztest := test + CABS(E[l-1]);
|
|
|
327 |
IF (ztest = test) THEN
|
|
|
328 |
E[l-1] := zero;
|
|
|
329 |
EXIT;
|
|
|
330 |
END;
|
|
|
331 |
INC(ll);
|
|
|
332 |
END;
|
|
|
333 |
|
|
|
334 |
IF (l # m-1) THEN
|
|
|
335 |
lp1 := l + 1;
|
|
|
336 |
mp1 := m + 1;
|
|
|
337 |
|
|
|
338 |
ls := l; (* Wg. Compilerwarnung *)
|
|
|
339 |
lls:=lp1;
|
|
|
340 |
LOOP
|
|
|
341 |
IF (lls > mp1) THEN EXIT; END;
|
|
|
342 |
ls := m - lls + lp1;
|
|
|
343 |
IF (ls = l) THEN
|
|
|
344 |
EXIT;
|
|
|
345 |
END;
|
|
|
346 |
test := 0.0;
|
|
|
347 |
IF (ls # m) THEN
|
|
|
348 |
test:= CABS(E[ls-1]); (* + test *)
|
|
|
349 |
END;
|
|
|
350 |
IF (ls # l+1) THEN
|
|
|
351 |
test := test + CABS(E[ls-1-1]);
|
|
|
352 |
END;
|
|
|
353 |
ztest := test + CABS(S[ls-1]);
|
|
|
354 |
IF (ztest = test) THEN
|
|
|
355 |
S[ls-1-1] := zero;
|
|
|
356 |
EXIT;
|
|
|
357 |
END;
|
|
|
358 |
INC(lls);
|
|
|
359 |
END;
|
|
|
360 |
|
|
|
361 |
IF (ls = l) THEN
|
|
|
362 |
kase := 3;
|
|
|
363 |
ELSIF (ls # m) THEN
|
|
|
364 |
kase := 2;
|
|
|
365 |
l := ls;
|
|
|
366 |
ELSE
|
|
|
367 |
kase := 1;
|
|
|
368 |
END;
|
|
|
369 |
ELSE
|
|
|
370 |
kase := 4;
|
|
|
371 |
END; (* IF *)
|
|
|
372 |
INC(l);
|
|
|
373 |
|
|
|
374 |
(* perform the task indicated by kase. *)
|
|
|
375 |
|
|
|
376 |
IF (kase = 1) THEN (* deflate negligible s(m). *)
|
|
|
377 |
|
|
|
378 |
mm1 := m - 1;
|
|
|
379 |
f := RE(E[m-1-1]);
|
|
|
380 |
E[m-1-1] := zero;
|
|
|
381 |
FOR kk:=l TO mm1 DO
|
|
|
382 |
k := mm1 - kk + l;
|
|
|
383 |
t1 := RE(S[k-1]);
|
|
|
384 |
drotg(t1,f,cs,sn);
|
|
|
385 |
S[k-1] := CMPLX(t1,0.0);
|
|
|
386 |
IF (k # l) THEN
|
|
|
387 |
f := -sn*RE(E[k-1-1]);
|
|
|
388 |
E[k-1-1] := scalarMult(cs,E[k-1-1]);
|
|
|
389 |
END;
|
|
|
390 |
IF (wantv) THEN
|
|
|
391 |
zdrot(P,V[k-1,1 -1],1,V[m-1,1 -1],1,cs,sn);
|
|
|
392 |
END; (* IF *)
|
|
|
393 |
END; (* FOR *)
|
|
|
394 |
|
|
|
395 |
ELSIF (kase = 2) THEN (* split at negligible s(l). *)
|
|
|
396 |
|
|
|
397 |
f := RE(E[l-1-1]);
|
|
|
398 |
E[l-1-1] := zero;
|
|
|
399 |
FOR k:=l TO m DO
|
|
|
400 |
t1 := RE(S[k-1]);
|
|
|
401 |
drotg(t1,f,cs,sn);
|
|
|
402 |
S[k-1] := CMPLX(t1,0.0);
|
|
|
403 |
f := -sn*RE(E[k-1]);
|
|
|
404 |
E[k-1] := scalarMult(cs,E[k-1]);
|
|
|
405 |
IF (wantu) THEN
|
|
|
406 |
zdrot(N,U[k-1,1 -1],1,U[l-1 -1,1 -1],1,cs,sn);
|
|
|
407 |
END;
|
|
|
408 |
END; (* FOR *)
|
|
|
409 |
|
|
|
410 |
ELSIF (kase = 3) THEN (* perform one qr step *)
|
|
|
411 |
|
|
|
412 |
(* calculate the shift. *)
|
|
|
413 |
scale := DMAX1(DMAX1(CABS(S[m -1]),CABS(S[m-1-1])),
|
|
|
414 |
DMAX1(CABS(E[m-1-1]),CABS(S[l -1])));
|
|
|
415 |
scale := DMAX1(scale,CABS(E[l-1]));
|
|
|
416 |
sm := RE(S[m -1]) / scale;
|
|
|
417 |
smm1 := RE(S[m-1-1]) / scale;
|
|
|
418 |
emm1 := RE(E[m-1-1]) / scale;
|
|
|
419 |
sl := RE(S[l -1]) / scale;
|
|
|
420 |
el := RE(E[l -1]) / scale;
|
|
|
421 |
b := ((smm1 + sm)*(smm1 - sm) + emm1*emm1) / 2.0;
|
|
|
422 |
c := sqr(sm*emm1);
|
|
|
423 |
shift := 0.0;
|
|
|
424 |
IF (b # 0.0) OR (c # 0.0) THEN
|
|
|
425 |
shift := sqrt(b*b + c);
|
|
|
426 |
IF (b < 0.0) THEN
|
|
|
427 |
shift := -shift;
|
|
|
428 |
END;
|
|
|
429 |
shift := c/(b+shift);
|
|
|
430 |
END; (* IF *)
|
|
|
431 |
f := (sl+sm)*(sl-sm) + shift;
|
|
|
432 |
g := sl*el;
|
|
|
433 |
|
|
|
434 |
(* chase zeros. *)
|
|
|
435 |
|
|
|
436 |
mm1 := m - 1;
|
|
|
437 |
FOR k:=l TO mm1 DO
|
|
|
438 |
drotg(f,g,cs,sn);
|
|
|
439 |
IF (k # l) THEN
|
|
|
440 |
E[k-1-1] := CMPLX(f,0.0);
|
|
|
441 |
END;
|
|
|
442 |
f := cs*RE(S[k-1]) + sn*RE(E[k-1]);
|
|
|
443 |
E[k-1] := scalarMult(cs,E[k-1]) - scalarMult(sn,S[k-1]);
|
|
|
444 |
g := sn*RE(S[k+1-1]);
|
|
|
445 |
S[k+1-1] := scalarMult(cs,S[k+1-1]);
|
|
|
446 |
IF wantv THEN
|
|
|
447 |
zdrot(P,V[k-1,1 -1],1,V[k+1 -1,1 -1],1,cs,sn);
|
|
|
448 |
END;
|
|
|
449 |
drotg(f,g,cs,sn);
|
|
|
450 |
S[k-1] := CMPLX(f,0.0);
|
|
|
451 |
f := cs*RE(E[k-1]) + sn*RE(S[k+1-1]);
|
|
|
452 |
S[k+1-1] := -scalarMult(sn,E[k-1]) + scalarMult(cs,S[k+1-1]);
|
|
|
453 |
g := sn*RE(E[k+1-1]);
|
|
|
454 |
E[k+1-1] := scalarMult(cs,E[k+1-1]);
|
|
|
455 |
IF wantu AND (k < N) THEN
|
|
|
456 |
zdrot(N,U[k-1,1 -1],1,U[k+1 -1,1 -1],1,cs,sn);
|
|
|
457 |
END;
|
|
|
458 |
END; (* FOR *)
|
|
|
459 |
E[m-1-1] := CMPLX(f,0.0);
|
|
|
460 |
INC(iter);
|
|
|
461 |
|
|
|
462 |
ELSIF (kase = 4) THEN
|
|
|
463 |
|
|
|
464 |
(* convergence, make the singular value positive *)
|
|
|
465 |
|
|
|
466 |
IF (RE(S[l-1]) < 0.0) THEN
|
|
|
467 |
S[l-1] := -S[l-1];
|
|
|
468 |
IF (wantv) THEN
|
|
|
469 |
zscal(P,CMPLX(-1.0,0.0),V[l-1,1 -1],1)
|
|
|
470 |
END;
|
|
|
471 |
END;
|
|
|
472 |
(* order the singular value. *)
|
|
|
473 |
WHILE (l # mm) AND (RE(S[l-1]) < RE(S[l+1-1])) DO
|
|
|
474 |
t := S[l-1];
|
|
|
475 |
S[l-1] := S[l+1-1];
|
|
|
476 |
S[l+1-1] := t;
|
|
|
477 |
IF wantv AND (l < P) THEN
|
|
|
478 |
zswap(P,V[l-1,1 -1],1,V[l+1 -1,1 -1],1)
|
|
|
479 |
END; (* IF *)
|
|
|
480 |
IF wantu AND (l < N) THEN
|
|
|
481 |
zswap(N,U[l-1,1 -1],1,U[l+1 -1,1 -1],1)
|
|
|
482 |
END; (* IF *)
|
|
|
483 |
INC(l);
|
|
|
484 |
END; (* WHILE *)
|
|
|
485 |
iter := 0;
|
|
|
486 |
DEC(m);
|
|
|
487 |
|
|
|
488 |
END; (* IF kase *)
|
|
|
489 |
|
|
|
490 |
ELSE
|
|
|
491 |
Info := m;
|
|
|
492 |
RETURN;
|
|
|
493 |
END; (* IF *)
|
|
|
494 |
END; (* FOR *)
|
|
|
495 |
END zSVDc;
|
|
|
496 |
|
|
|
497 |
PROCEDURE Takagi( N : INTEGER;
|
|
|
498 |
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
|
|
|
499 |
VAR D : ARRAY OF LONGREAL;
|
|
|
500 |
VAR U : ARRAY OF ARRAY OF LONGCOMPLEX;
|
|
|
501 |
sort : INTEGER);
|
|
|
502 |
|
|
|
503 |
PROCEDURE SQ(c : LONGCOMPLEX) : LONGREAL;
|
|
|
504 |
|
|
|
505 |
BEGIN
|
|
|
506 |
RETURN RE(c*conj(c));
|
|
|
507 |
END SQ;
|
|
|
508 |
|
|
|
509 |
CONST eps = 1.0E+01*MachEps;
|
|
|
510 |
MaxI = MAX(INTEGER);
|
|
|
511 |
|
|
|
512 |
VAR p,q,j : INTEGER;
|
|
|
513 |
red,off,thresh : LONGREAL;
|
|
|
514 |
sqp,sqq,t,invc : LONGREAL;
|
|
|
515 |
f,x,y : LONGCOMPLEX;
|
|
|
516 |
ev1,ev2 : POINTER TO ARRAY [0..MaxI-1] OF LONGCOMPLEX;
|
|
|
517 |
sweep : INTEGER;
|
|
|
518 |
BEGIN
|
|
|
519 |
ALLOCATE(ev1,N*SIZELONGCMPLX);
|
|
|
520 |
ALLOCATE(ev2,N*SIZELONGCMPLX);
|
|
|
521 |
|
|
|
522 |
IF (ev1 = NIL) OR (ev2 = NIL) THEN
|
|
|
523 |
Errors.ErrOut("Kein Freispeicher vorhanden (Tagati)");
|
|
|
524 |
IF (ev1 # NIL) THEN DEALLOCATE(ev1,N*SIZELONGCMPLX); END;
|
|
|
525 |
RETURN;
|
|
|
526 |
END;
|
|
|
527 |
|
|
|
528 |
FOR p:=0 TO N-1 DO
|
|
|
529 |
ev1^[p] := zero;
|
|
|
530 |
ev2^[p] := A[p,p];
|
|
|
531 |
END;
|
|
|
532 |
|
|
|
533 |
FOR p:=0 TO N-1 DO
|
|
|
534 |
FOR q:=0 TO N-1 DO
|
|
|
535 |
U[p,q] := zero;
|
|
|
536 |
END;
|
|
|
537 |
U[p,p] := one;
|
|
|
538 |
END;
|
|
|
539 |
|
|
|
540 |
red := 0.04 / CardPot(VAL(LONGREAL,N),4);
|
|
|
541 |
|
|
|
542 |
sweep:=0;
|
|
|
543 |
LOOP
|
|
|
544 |
INC(sweep);
|
|
|
545 |
|
|
|
546 |
IF (sweep > 50) THEN
|
|
|
547 |
Errors.ErrOut("Bad convergence in TakagiFactor");
|
|
|
548 |
EXIT;
|
|
|
549 |
END;
|
|
|
550 |
|
|
|
551 |
off := 0.0;
|
|
|
552 |
FOR q:=1 TO N-1 DO
|
|
|
553 |
FOR p:=0 TO q-1 DO
|
|
|
554 |
off:=off + SQ(A[q,p]);
|
|
|
555 |
END;
|
|
|
556 |
END;
|
|
|
557 |
IF (off <= eps*eps) THEN
|
|
|
558 |
EXIT;
|
|
|
559 |
END;
|
|
|
560 |
thresh := 0.0;
|
|
|
561 |
IF (sweep < 4) THEN
|
|
|
562 |
thresh := off*red
|
|
|
563 |
END;
|
|
|
564 |
|
|
|
565 |
FOR q:=1 TO N-1 DO
|
|
|
566 |
FOR p:=0 TO q-1 DO
|
|
|
567 |
off := SQ(A[q,p]);
|
|
|
568 |
sqp := SQ(ev2^[p]);
|
|
|
569 |
sqq := SQ(ev2^[q]);
|
|
|
570 |
IF (sweep > 4) AND (off < eps*(sqp+sqq)) THEN
|
|
|
571 |
A[q,p] := zero;
|
|
|
572 |
ELSIF (off > thresh) THEN
|
|
|
573 |
t := 0.5*ABS(sqp - sqq);
|
|
|
574 |
IF (t > 0.0) THEN
|
|
|
575 |
f := scalarMult(sign2(1.0,sqp-sqq),(ev2^[q]*conj(A[q,p]) +
|
|
|
576 |
conj(ev2^[p])*A[q,p]));
|
|
|
577 |
ELSE
|
|
|
578 |
f := one;
|
|
|
579 |
IF (sqp # 0) THEN
|
|
|
580 |
f := LongComplexMath.sqrt(ev2^[q] / ev2^[p])
|
|
|
581 |
END;
|
|
|
582 |
END;
|
|
|
583 |
t:=t + sqrt(t*t + SQ(f));
|
|
|
584 |
f:=f / CMPLX(t,0.0);
|
|
|
585 |
|
|
|
586 |
ev1^[p] := ev1^[p] + A[q,p]*conj(f);
|
|
|
587 |
ev2^[p] := A[p,p] + ev1^[p];
|
|
|
588 |
ev1^[q] := ev1^[q] - A[q,p]*f;
|
|
|
589 |
ev2^[q] := A[q,q] + ev1^[q];
|
|
|
590 |
|
|
|
591 |
t := SQ(f);
|
|
|
592 |
invc := sqrt(t +1.0);
|
|
|
593 |
f:=f / CMPLX(invc,0.0);
|
|
|
594 |
t:=t / (invc*(invc + 1.0));
|
|
|
595 |
|
|
|
596 |
FOR j:=0 TO p-1 DO
|
|
|
597 |
x := A[p,j];
|
|
|
598 |
y := A[q,j];
|
|
|
599 |
A[p,j] := x + (conj(f)*y - CMPLX(t,0.0)*x);
|
|
|
600 |
A[q,j] := y - (f *x + CMPLX(t,0.0)*y);
|
|
|
601 |
END;
|
|
|
602 |
|
|
|
603 |
FOR j:=p+1 TO q-1 DO
|
|
|
604 |
x := A[j,p];
|
|
|
605 |
y := A[q,j];
|
|
|
606 |
A[j,p] := x + (conj(f)*y - CMPLX(t,0.0)*x);
|
|
|
607 |
A[q,j] := y - (f *x + CMPLX(t,0.0)*y);
|
|
|
608 |
END;
|
|
|
609 |
|
|
|
610 |
FOR j:=q+1 TO N-1 DO
|
|
|
611 |
x := A[j,p];
|
|
|
612 |
y := A[j,q];
|
|
|
613 |
A[j,p] := x + (conj(f)*y - CMPLX(t,0.0)*x);
|
|
|
614 |
A[j,q] := y - (f *x + CMPLX(t,0.0)*y);
|
|
|
615 |
END;
|
|
|
616 |
|
|
|
617 |
A[q,p] := zero;
|
|
|
618 |
|
|
|
619 |
FOR j:=0 TO N-1 DO
|
|
|
620 |
x := U[j,p];
|
|
|
621 |
y := U[j,q];
|
|
|
622 |
U[j,p] := x + (f *y - CMPLX(t,0.0)*x);
|
|
|
623 |
U[j,q] := y - (conj(f)*x + CMPLX(t,0.0)*y);
|
|
|
624 |
END;
|
|
|
625 |
|
|
|
626 |
END; (* IF *)
|
|
|
627 |
END;
|
|
|
628 |
END;
|
|
|
629 |
|
|
|
630 |
FOR p:=0 TO N-1 DO
|
|
|
631 |
ev1^[p] := zero;
|
|
|
632 |
A[p,p] := ev2^[p];
|
|
|
633 |
END;
|
|
|
634 |
|
|
|
635 |
END; (* LOOP *)
|
|
|
636 |
|
|
|
637 |
(* make the diagonal elements nonnegative *)
|
|
|
638 |
|
|
|
639 |
FOR p:=0 TO N-1 DO
|
|
|
640 |
D[p] := CABS(A[p,p]);
|
|
|
641 |
IF (D[p] > eps) AND (D[p] # RE(A[p,p])) THEN
|
|
|
642 |
f := LongComplexMath.sqrt(A[p,p] / CMPLX(D[p],0.0));
|
|
|
643 |
FOR q:=0 TO N-1 DO
|
|
|
644 |
U[q,p] := U[q,p]*f;
|
|
|
645 |
END;
|
|
|
646 |
END; (* IF *)
|
|
|
647 |
END;
|
|
|
648 |
|
|
|
649 |
IF (sort # 0) THEN (* sort the eigenvalues *)
|
|
|
650 |
FOR p:=0 TO N-2 DO
|
|
|
651 |
j := p;
|
|
|
652 |
t := D[p];
|
|
|
653 |
FOR q := p+1 TO N-1 DO
|
|
|
654 |
IF (VAL(LONGREAL,sort)*(t - D[q]) > 0.0) THEN
|
|
|
655 |
j := q;
|
|
|
656 |
t := D[q];
|
|
|
657 |
END; (* IF *)
|
|
|
658 |
END;
|
|
|
659 |
IF (j # p) THEN
|
|
|
660 |
D[j] := D[p];
|
|
|
661 |
D[p] := t;
|
|
|
662 |
FOR q:=0 TO N-1 DO
|
|
|
663 |
x := U[q,p];
|
|
|
664 |
U[q,p] := U[q,j];
|
|
|
665 |
U[q,j] := x;
|
|
|
666 |
END;
|
|
|
667 |
END; (* IF *)
|
|
|
668 |
END;
|
|
|
669 |
END; (* IF sort *)
|
|
|
670 |
|
|
|
671 |
DEALLOCATE(ev2,N*SIZELONGCMPLX);
|
|
|
672 |
DEALLOCATE(ev1,N*SIZELONGCMPLX);
|
|
|
673 |
END Takagi;
|
|
|
674 |
|
|
|
675 |
END SVDLib3.
|