17
|
1 |
#======================================================================
|
|
2 |
# L I B Q Z . P L
|
|
3 |
# doc: Thu Mar 12 15:23:15 2015
|
|
4 |
# dlm: Sun Mar 15 19:26:50 2015
|
|
5 |
# (c) 2015 A.M. Thurnherr
|
|
6 |
# uE-Info: 14 27 NIL 0 0 72 10 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# adaptation of EISPACK routines
|
|
10 |
# www.netlib.org/eispack
|
|
11 |
|
|
12 |
# HISTORY:
|
|
13 |
# Mar 12, 2015: - created
|
|
14 |
# Mar 15, 2015: - debugging
|
|
15 |
|
|
16 |
use strict vars;
|
|
17 |
|
|
18 |
my($N); # size of input matrices A & B
|
|
19 |
|
|
20 |
#----------------------------------------------------------------------
|
|
21 |
# eig(\@A,\@B,\@evRe,\@evIm,\@V)
|
|
22 |
# - @ev{Re,Im} contain generalized eigenvalues
|
|
23 |
# - @evec contains corresponding right eigenvectors
|
|
24 |
# - A*V = B*V*D, where D is diagonal matrix of eigenvalues
|
|
25 |
#----------------------------------------------------------------------
|
|
26 |
|
|
27 |
sub eig($$$$$)
|
|
28 |
{
|
|
29 |
my($aR,$bR,$erR,$eiR,$zR) = @_; # args passed as refs
|
|
30 |
my(@alphaR,@alphaI,@beta); # intermediate data
|
|
31 |
|
|
32 |
$N = scalar(@{$aR});
|
|
33 |
croak(sprintf("eig(A,B): A(%dx%d) & B(%dx%d) must be matching square matrices\n",
|
|
34 |
$N,scalar(@{$aR->[0]}),scalar(@{$bR}),scalar(@{$bR->[0]})))
|
|
35 |
unless (@{$bR} == $N) && (@{$aR->[0]} == $N) && (@{$bR->[0]} == $N);
|
|
36 |
|
|
37 |
QZhes($aR,$bR,$zR); # reduce A/B to upper Hessenberg/triangular forms
|
|
38 |
croak("QZit(): convergence failed\n")
|
|
39 |
unless (QZit($aR,$bR,$zR) == 0); # reduce Hess A to quasi-triangular form
|
|
40 |
QZval($aR,$bR,$zR,\@alphaR,\@alphaI,\@beta); # reduce A further
|
|
41 |
QZvec($aR,$bR,$zR,\@alphaR,\@alphaI,\@beta); # compute eigenvectors & eigenvalues
|
|
42 |
|
|
43 |
for (my($i)=0; $i<$N; $i++) {
|
|
44 |
if ($beta[$i]==0 && $alphaR[$i]==0) {
|
|
45 |
$erR->[$i] = nan;
|
|
46 |
$eiR->[$i] = nan;
|
|
47 |
} else {
|
|
48 |
$erR->[$i] = $alphaR[$i] / $beta[$i];
|
|
49 |
$eiR->[$i] = $alphaI[$i] / $beta[$i];
|
|
50 |
# print(STDERR "gev[$i] = $erR->[$i]\n");
|
|
51 |
}
|
|
52 |
}
|
|
53 |
}
|
|
54 |
|
|
55 |
#----------------------------------------------------------------------
|
|
56 |
# EISPACK routines
|
|
57 |
#----------------------------------------------------------------------
|
|
58 |
|
|
59 |
#----------------------------------------------------------------------
|
|
60 |
# QZhes(\@A,\@B,\@Z)
|
|
61 |
# - first step in QZ algorithm (Moler & Stewart, SIAM JNA, 1973)
|
|
62 |
#----------------------------------------------------------------------
|
|
63 |
|
|
64 |
sub QZhes($$$)
|
|
65 |
{
|
|
66 |
my($aR,$bR,$zR) = @_;
|
|
67 |
my($i,$j,$k,$l,$l1,$lb,$nk1);
|
|
68 |
my($r,$s,$t,$u1,$u2,$v1,$v2,$rho);
|
|
69 |
|
|
70 |
croak("QZhes: need at least 3x3 matrices\n")
|
|
71 |
unless ($N >= 2);
|
|
72 |
|
|
73 |
for ($j=0; $j<$N; $j++) { # init Z to identity matrix
|
|
74 |
for ($i=0; $i<$N; $i++) {
|
|
75 |
$zR->[$i][$j] = 0;
|
|
76 |
}
|
|
77 |
$zR->[$j][$j] = 1;
|
|
78 |
}
|
|
79 |
|
|
80 |
for ($l=0; $l<$N-1; $l++) {
|
|
81 |
$l1 = $l + 1;
|
|
82 |
for ($s=0,$i=$l1; $i<$N; $i++) {
|
|
83 |
$s += abs($bR->[$i][$l]);
|
|
84 |
}
|
|
85 |
next if ($s == 0);
|
|
86 |
$s += abs($bR->[$l][$l]);
|
|
87 |
|
|
88 |
for ($r=0,$i=$l; $i<$N; $i++)
|
|
89 |
{
|
|
90 |
$bR->[$i][$l] /= $s;
|
|
91 |
$r += $bR->[$i][$l]**2;
|
|
92 |
}
|
|
93 |
|
|
94 |
$r = SIGN(sqrt($r),$bR->[$l][$l]);
|
|
95 |
$bR->[$l][$l] += $r;
|
|
96 |
$rho = $r * $bR->[$l][$l];
|
|
97 |
|
|
98 |
for ($j=$l1; $j<$N; $j++) {
|
|
99 |
for ($t=0,$i=$l; $i<$N; $i++) {
|
|
100 |
$t += $bR->[$i][$l] * $bR->[$i][$j];
|
|
101 |
}
|
|
102 |
$t = -$t / $rho;
|
|
103 |
for ($i=$l; $i<$N; $i++) {
|
|
104 |
$bR->[$i][$j] += $t * $bR->[$i][$l];
|
|
105 |
}
|
|
106 |
}
|
|
107 |
|
|
108 |
for ($j=0; $j<$N; $j++) {
|
|
109 |
for ($t=0,$i=$l; $i<$N; $i++) {
|
|
110 |
$t += $bR->[$i][$l] * $aR->[$i][$j];
|
|
111 |
}
|
|
112 |
$t = -$t / $rho;
|
|
113 |
for ($i=$l; $i<$N; $i++) {
|
|
114 |
$aR->[$i][$j] += $t * $bR->[$i][$l];
|
|
115 |
}
|
|
116 |
}
|
|
117 |
|
|
118 |
$bR->[$l][$l] = -$s * $r;
|
|
119 |
for ($i=$l1; $i<$N; $i++) {
|
|
120 |
$bR->[$i][$l] = 0;
|
|
121 |
}
|
|
122 |
}
|
|
123 |
|
|
124 |
for ($k=0; $k<$N-2; $k++) {
|
|
125 |
$nk1 = $N - 2 - $k;
|
|
126 |
|
|
127 |
for ($lb=0; $lb<$nk1; $lb++) {
|
|
128 |
$l = $N - $lb - 2;
|
|
129 |
$l1 = $l + 1;
|
|
130 |
|
|
131 |
$s = (abs($aR->[$l][$k])) + (abs($aR->[$l1][$k]));
|
|
132 |
next if ($s == 0);
|
|
133 |
|
|
134 |
$u1 = $aR->[$l][$k] / $s;
|
|
135 |
$u2 = $aR->[$l1][$k] / $s;
|
|
136 |
$r = SIGN(sqrt($u1**2+$u2**2),$u1);
|
|
137 |
$v1 = -($u1 + $r) / $r;
|
|
138 |
$v2 = -$u2 / $r;
|
|
139 |
$u2 = $v2 / $v1;
|
|
140 |
|
|
141 |
for ($j=$k; $j<$N; $j++) {
|
|
142 |
$t = $aR->[$l][$j] + $u2 * $aR->[$l1][$j];
|
|
143 |
$aR->[$l][$j] += $t * $v1;
|
|
144 |
$aR->[$l1][$j] += $t * $v2;
|
|
145 |
}
|
|
146 |
|
|
147 |
$aR->[$l1][$k] = 0;
|
|
148 |
|
|
149 |
for ($j=$l; $j<$N; $j++) {
|
|
150 |
$t = $bR->[$l][$j] + $u2 * $bR->[$l1][$j];
|
|
151 |
$bR->[$l][$j] += $t * $v1;
|
|
152 |
$bR->[$l1][$j] += $t * $v2;
|
|
153 |
}
|
|
154 |
|
|
155 |
$s = (abs($bR->[$l1][$l1])) + (abs($bR->[$l1][$l]));
|
|
156 |
next if ($s == 0);
|
|
157 |
|
|
158 |
$u1 = $bR->[$l1][$l1] / $s;
|
|
159 |
$u2 = $bR->[$l1][$l] / $s;
|
|
160 |
$r = SIGN(sqrt($u1**2 + $u2**2),$u1);
|
|
161 |
$v1 = -($u1 + $r) / $r;
|
|
162 |
$v2 = -$u2 / $r;
|
|
163 |
$u2 = $v2 / $v1;
|
|
164 |
|
|
165 |
for ($i=0; $i<=$l1; $i++) { # overwrite B with upper triangular form
|
|
166 |
$t = $bR->[$i][$l1] + $u2 * $bR->[$i][$l];
|
|
167 |
$bR->[$i][$l1] += $t * $v1;
|
|
168 |
$bR->[$i][$l] += $t * $v2;
|
|
169 |
}
|
|
170 |
$bR->[$l1][$l] = 0;
|
|
171 |
|
|
172 |
for ($i=0; $i<$N; $i++) { # overwrite A with upper Hessenberg form
|
|
173 |
$t = $aR->[$i][$l1] + $u2 * $aR->[$i][$l];
|
|
174 |
$aR->[$i][$l1] += $t * $v1;
|
|
175 |
$aR->[$i][$l] += $t * $v2;
|
|
176 |
}
|
|
177 |
|
|
178 |
for ($i=0; $i<$N; $i++) { # define matrix Z, used for eigenvectors
|
|
179 |
$t = $zR->[$i][$l1] + $u2 * $zR->[$i][$l];
|
|
180 |
$zR->[$i][$l1] += $t * $v1;
|
|
181 |
$zR->[$i][$l] += $t * $v2;
|
|
182 |
}
|
|
183 |
} # for ($lb=0; $lb<$nk1; $lb++)
|
|
184 |
} # for ($k=0; $k<$N-2; $k++)
|
|
185 |
}
|
|
186 |
|
|
187 |
#----------------------------------------------------------------------
|
|
188 |
# QZit(\@A,\@B,\@Z)
|
|
189 |
# - second step in QZ algorithm (Moler & Stewart, SIAM JNA, 1973)
|
|
190 |
# - !0 return value indicates that convergence failed
|
|
191 |
#----------------------------------------------------------------------
|
|
192 |
|
|
193 |
sub QZit($$$$)
|
|
194 |
{
|
|
195 |
my($aR,$bR,$zR) = @_;
|
|
196 |
|
|
197 |
my($i,$j,$k);
|
|
198 |
my($r,$s,$t,$a1,$a2);
|
|
199 |
my($k1,$k2,$l1,$ll);
|
|
200 |
my($u1,$u2,$u3);
|
|
201 |
my($v1,$v2,$v3);
|
|
202 |
my($a11,$a12,$a21,$a22,$a33,$a34,$a43,$a44);
|
|
203 |
my($b11,$b12,$b22,$b33,$b34,$b44);
|
|
204 |
my($na,$en,$ld);
|
|
205 |
my($ep,$km1,$ani,$bni);
|
|
206 |
my($ish,$itn,$its,$enm2,$lor1,$epsa,$epsb,$enorn,$notlas);
|
|
207 |
my($l,$a3,$sh,$lm1,$anorm,$bnorm) = (0,0,0,0,0,0);
|
|
208 |
|
|
209 |
for ($i=0; $i<$N; $i++) {
|
|
210 |
$ani = $bni = 0;
|
|
211 |
$ani = abs($aR->[$i][$i-1])
|
|
212 |
unless ($i == 0);
|
|
213 |
for ($j=$i; $j<$N; $j++) {
|
|
214 |
$ani += abs($aR->[$i][$j]);
|
|
215 |
$bni += abs($bR->[$i][$j]);
|
|
216 |
}
|
|
217 |
$anorm = $ani if ($ani > $anorm);
|
|
218 |
$bnorm = $bni if ($bni > $bnorm);
|
|
219 |
}
|
|
220 |
$anorm = 1 if ($anorm == 0);
|
|
221 |
$bnorm = 1 if ($bnorm == 0);
|
|
222 |
|
|
223 |
$ep = 1.11022302462516e-16; # $EPS=1; $EPS /= 2 while 0.5 + $EPS/2 > 0.5;
|
|
224 |
$epsa = $ep * $anorm;
|
|
225 |
$epsb = $ep * $bnorm;
|
|
226 |
|
|
227 |
$lor1 = 0;
|
|
228 |
$enorn = $N;
|
|
229 |
$en = $N - 1;
|
|
230 |
$itn = $N * 30;
|
|
231 |
|
|
232 |
L60:
|
|
233 |
goto L1001 if ($en <= 1);
|
|
234 |
$its = 0;
|
|
235 |
$na = $en - 1;
|
|
236 |
$enm2 = $na;
|
|
237 |
|
|
238 |
L70:
|
|
239 |
$ish = 2;
|
|
240 |
for ($ll=0; $ll<=$en; $ll++) {
|
|
241 |
$lm1 = $en - $ll - 1;
|
|
242 |
$l = $lm1 + 1;
|
|
243 |
goto L95 if ($l == 0);
|
|
244 |
last if ((abs($aR->[$l][$lm1])) <= $epsa)
|
|
245 |
}
|
|
246 |
|
|
247 |
L90:
|
|
248 |
$aR->[$l][$lm1] = 0;
|
|
249 |
goto L95 if ($l < $na);
|
|
250 |
$en = $lm1;
|
|
251 |
goto L60;
|
|
252 |
|
|
253 |
L95:
|
|
254 |
$ld = $l;
|
|
255 |
|
|
256 |
L100:
|
|
257 |
$l1 = $l + 1;
|
|
258 |
$b11 = $bR->[$l][$l];
|
|
259 |
goto L120 if (abs($b11) > $epsb);
|
|
260 |
$bR->[$l][$l] = 0;
|
|
261 |
$s = (abs($aR->[$l][$l]) + abs($aR->[$l1][$l]));
|
|
262 |
$u1 = $aR->[$l][$l] / $s;
|
|
263 |
$u2 = $aR->[$l1][$l] / $s;
|
|
264 |
$r = SIGN(sqrt($u1**2 + $u2**2),$u1);
|
|
265 |
$v1 = -($u1 + $r) / $r;
|
|
266 |
$v2 = -$u2 / $r;
|
|
267 |
$u2 = $v2 / $v1;
|
|
268 |
|
|
269 |
for ($j=$l; $j<$enorn; $j++) {
|
|
270 |
$t = $aR->[$l][$j] + $u2 * $aR->[$l1][$j];
|
|
271 |
$aR->[$l][$j] += $t * $v1;
|
|
272 |
$aR->[$l1][$j] += $t * $v2;
|
|
273 |
|
|
274 |
$t = $bR->[$l][$j] + $u2 * $bR->[$l1][$j];
|
|
275 |
$bR->[$l][$j] += $t * $v1;
|
|
276 |
$bR->[$l1][$j] += $t * $v2;
|
|
277 |
}
|
|
278 |
$aR->[$l][$lm1] = -$aR->[$l][$lm1]
|
|
279 |
if ($l != 0);
|
|
280 |
|
|
281 |
$lm1 = $l;
|
|
282 |
$l = $l1;
|
|
283 |
goto L90;
|
|
284 |
|
|
285 |
L120:
|
|
286 |
$a11 = $aR->[$l][$l] / $b11;
|
|
287 |
$a21 = $aR->[$l1][$l] / $b11;
|
|
288 |
|
|
289 |
goto L140 if ($ish == 1);
|
|
290 |
goto L1000 if ($itn == 0);
|
|
291 |
goto L155 if ($its == 10);
|
|
292 |
|
|
293 |
$b22 = $bR->[$l1][$l1];
|
|
294 |
$b22 = $epsb if (abs($b22) < $epsb);
|
|
295 |
$b33 = $bR->[$na][$na];
|
|
296 |
$b33 = $epsb if (abs($b33) < $epsb);
|
|
297 |
$b44 = $bR->[$en][$en];
|
|
298 |
$b44 = $epsb if (abs($b44) < $epsb);
|
|
299 |
|
|
300 |
$a33 = $aR->[$na][$na] / $b33;
|
|
301 |
$a34 = $aR->[$na][$en] / $b44;
|
|
302 |
$a43 = $aR->[$en][$na] / $b33;
|
|
303 |
$a44 = $aR->[$en][$en] / $b44;
|
|
304 |
$b34 = $bR->[$na][$en] / $b44;
|
|
305 |
|
|
306 |
$t = ($a43 * $b34 - $a33 - $a44) / 2;
|
|
307 |
$r = $t * $t + $a34 * $a43 - $a33 * $a44;
|
|
308 |
goto L150 if ($r < 0);
|
|
309 |
|
|
310 |
$ish = 1;
|
|
311 |
$r = sqrt($r);
|
|
312 |
$sh = -$t + $r;
|
|
313 |
$s = -$t - $r;
|
|
314 |
$sh = $s if (abs($s-$a44) < abs($sh-$a44));
|
|
315 |
|
|
316 |
for ($ll=$ld; $ll<$enm2; $ll++) {
|
|
317 |
$l = $enm2 + $ld - $ll - 1;
|
|
318 |
goto L140 if ($l == $ld);
|
|
319 |
$lm1 = $l - 1;
|
|
320 |
$l1 = $l + 1;
|
|
321 |
$t = $aR->[$l+1][$l+1];
|
|
322 |
$t -= $sh * $bR->[$l][$l] if (abs($bR->[$l][$l]) > $epsb);
|
|
323 |
goto L100 if (abs($aR->[$l][$lm1]) <= (abs($t / $aR->[$l1][$l])) * $epsa);
|
|
324 |
}
|
|
325 |
|
|
326 |
L140:
|
|
327 |
$a1 = $a11 - $sh;
|
|
328 |
$a2 = $a21;
|
|
329 |
$aR->[$l][$lm1] = -$aR->[$l][$lm1] if ($l != $ld);
|
|
330 |
goto L160;
|
|
331 |
|
|
332 |
L150:
|
|
333 |
$a12 = $aR->[$l][$l1] / $b22;
|
|
334 |
$a22 = $aR->[$l1][$l1] / $b22;
|
|
335 |
$b12 = $bR->[$l][$l1] / $b22;
|
|
336 |
$a1 = (($a33 - $a11) * ($a44 - $a11) - $a34 * $a43 + $a43 * $b34 * $a11) / $a21 + $a12 - $a11 * $b12;
|
|
337 |
$a2 = $a22 - $a11 - $a21 * $b12 - ($a33 - $a11) - ($a44 - $a11) + $a43 * $b34;
|
|
338 |
$a3 = $aR->[$l1+1][$l] / $b22;
|
|
339 |
goto L160;
|
|
340 |
|
|
341 |
L155:
|
|
342 |
$a1 = 0;
|
|
343 |
$a2 = 1;
|
|
344 |
$a3 = 1.1605;
|
|
345 |
|
|
346 |
L160:
|
|
347 |
$its++;
|
|
348 |
$itn--;
|
|
349 |
$lor1 = $ld;
|
|
350 |
for ($k=$l; $k<=$na; $k++) {
|
|
351 |
|
|
352 |
$notlas = ($k!=$na) && ($ish==2);
|
|
353 |
$k1 = $k + 1;
|
|
354 |
$k2 = $k + 2;
|
|
355 |
$km1 = MAX($k,$l+1)-1;
|
|
356 |
$ll = MIN($en,$k1+$ish);
|
|
357 |
goto L190 if ($notlas);
|
|
358 |
|
|
359 |
if ($k != $l) {
|
|
360 |
$a1 = $aR->[$k,$km1];
|
|
361 |
$a2 = $aR->[$k1,$km1];
|
|
362 |
}
|
|
363 |
|
|
364 |
$s = abs($a1) + abs($a2);
|
|
365 |
goto L70 if ($s == 0);
|
|
366 |
$u1 = $a1 / $s;
|
|
367 |
$u2 = $a2 / $s;
|
|
368 |
$r = SIGN(sqrt($u1**2 + $u2**2), $u1);
|
|
369 |
$v1 = -($u1 + $r) / $r;
|
|
370 |
$v2 = -$u2 / $r;
|
|
371 |
$u2 = $v2 / $v1;
|
|
372 |
|
|
373 |
for ($j=$km1; $j<$enorn; $j++) {
|
|
374 |
$t = $aR->[$k][$j] + $u2 * $aR->[$k1][$j];
|
|
375 |
$aR->[$k][$j] += $t * $v1;
|
|
376 |
$aR->[$k1][$j] += $t * $v2;
|
|
377 |
$t = $bR->[$k][$j] + $u2 * $bR->[$k1][$j];
|
|
378 |
$bR->[$k][$j] += $t * $v1;
|
|
379 |
$bR->[$k1][$j] += $t * $v2;
|
|
380 |
}
|
|
381 |
|
|
382 |
$aR->[$k1,$km1] = 0 if ($k != $l);
|
|
383 |
goto L240;
|
|
384 |
|
|
385 |
L190:
|
|
386 |
goto L200 if ($k == $l);
|
|
387 |
$a1 = $aR->[$k,$km1];
|
|
388 |
$a2 = $aR->[$k1,$km1];
|
|
389 |
$a3 = $aR->[$k2][$km1];
|
|
390 |
|
|
391 |
L200:
|
|
392 |
$s = abs($a1) + abs($a2) + abs($a3);
|
|
393 |
next if ($s == 0);
|
|
394 |
$u1 = $a1 / $s;
|
|
395 |
$u2 = $a2 / $s;
|
|
396 |
$u3 = $a3 / $s;
|
|
397 |
$r = SIGN(sqrt($u1**2 + $u2**2 + $u3**2), $u1);
|
|
398 |
$v1 = -($u1 + $r) / $r;
|
|
399 |
$v2 = -$u2 / $r;
|
|
400 |
$v3 = -$u3 / $r;
|
|
401 |
$u2 = $v2 / $v1;
|
|
402 |
$u3 = $v3 / $v1;
|
|
403 |
|
|
404 |
for ($j=$km1; $j<$enorn; $j++) {
|
|
405 |
$t = $aR->[$k][$j] + $u2 * $aR->[$k1][$j] + $u3 * $aR->[$k2][$j];
|
|
406 |
$aR->[$k][$j] += $t * $v1;
|
|
407 |
$aR->[$k1][$j] += $t * $v2;
|
|
408 |
$aR->[$k2][$j] += $t * $v3;
|
|
409 |
|
|
410 |
$t = $bR->[$k][$j] + $u2 * $bR->[$k1][$j] + $u3 * $bR->[$k2][$j];
|
|
411 |
$bR->[$k][$j] += $t * $v1;
|
|
412 |
$bR->[$k1][$j] += $t * $v2;
|
|
413 |
$bR->[$k2][$j] += $t * $v3;
|
|
414 |
}
|
|
415 |
|
|
416 |
goto L220 if ($k == $l);
|
|
417 |
$aR->[$k1,$km1] = $aR->[$k2][$km1] = 0;
|
|
418 |
|
|
419 |
L220:
|
|
420 |
$s = (abs($bR->[$k2][$k2])) + (abs($bR->[$k2][$k1])) + (abs($bR->[$k2][$k]));
|
|
421 |
goto L240 if ($s == 0);
|
|
422 |
$u1 = $bR->[$k2][$k2] / $s;
|
|
423 |
$u2 = $bR->[$k2][$k1] / $s;
|
|
424 |
$u3 = $bR->[$k2][$k] / $s;
|
|
425 |
$r = SIGN(sqrt($u1**2 + $u2**2 + $u3**2), $u1);
|
|
426 |
$v1 = -($u1 + $r) / $r;
|
|
427 |
$v2 = -$u2 / $r;
|
|
428 |
$v3 = -$u3 / $r;
|
|
429 |
$u2 = $v2 / $v1;
|
|
430 |
$u3 = $v3 / $v1;
|
|
431 |
|
|
432 |
for ($i=$lor1; $i<$ll+1; $i++) {
|
|
433 |
$t = $aR->[$i][$k2] + $u2 * $aR->[$i][$k1] + $u3 * $aR->[$i][$k];
|
|
434 |
$aR->[$i][$k2] += $t * $v1;
|
|
435 |
$aR->[$i][$k1] += $t * $v2;
|
|
436 |
$aR->[$i][$k] += $t * $v3;
|
|
437 |
$t = $bR->[$i][$k2] + $u2 * $bR->[$i][$k1] + $u3 * $bR->[$i][$k];
|
|
438 |
$bR->[$i][$k2] += $t * $v1;
|
|
439 |
$bR->[$i][$k1] += $t * $v2;
|
|
440 |
$bR->[$i][$k] += $t * $v3;
|
|
441 |
}
|
|
442 |
|
|
443 |
$bR->[$k2][$k] = $bR->[$k2][$k1] = 0;
|
|
444 |
|
|
445 |
for ($i=0; $i<$N; $i++) {
|
|
446 |
$t = $zR->[$i][$k2] + $u2 * $zR->[$i][$k1] + $u3 * $zR->[$i][$k];
|
|
447 |
$zR->[$i][$k2] += $t * $v1;
|
|
448 |
$zR->[$i][$k1] += $t * $v2;
|
|
449 |
$zR->[$i][$k] += $t * $v3;
|
|
450 |
}
|
|
451 |
|
|
452 |
L240:
|
|
453 |
$s = (abs($bR->[$k1][$k1])) + (abs($bR->[$k1][$k]));
|
|
454 |
next if ($s == 0);
|
|
455 |
$u1 = $bR->[$k1][$k1] / $s;
|
|
456 |
$u2 = $bR->[$k1][$k] / $s;
|
|
457 |
$r = SIGN(sqrt($u1**2 + $u2**2), $u1);
|
|
458 |
$v1 = -($u1 + $r) / $r;
|
|
459 |
$v2 = -$u2 / $r;
|
|
460 |
$u2 = $v2 / $v1;
|
|
461 |
|
|
462 |
for ($i=$lor1; $i<$ll+1; $i++) {
|
|
463 |
$t = $aR->[$i][$k1] + $u2 * $aR->[$i][$k];
|
|
464 |
$aR->[$i][$k1] += $t * $v1;
|
|
465 |
$aR->[$i][$k] += $t * $v2;
|
|
466 |
$t = $bR->[$i][$k1] + $u2 * $bR->[$i][$k];
|
|
467 |
$bR->[$i][$k1] += $t * $v1;
|
|
468 |
$bR->[$i][$k] += $t * $v2;
|
|
469 |
}
|
|
470 |
|
|
471 |
$bR->[$k1][$k] = 0;
|
|
472 |
|
|
473 |
for ($i=0; $i<$N; $i++) {
|
|
474 |
$t = $zR->[$i][$k1] + $u2 * $zR->[$i][$k];
|
|
475 |
$zR->[$i][$k1] += $t * $v1;
|
|
476 |
$zR->[$i][$k] += $t * $v2;
|
|
477 |
}
|
|
478 |
|
|
479 |
} # for loop beginning at L160
|
|
480 |
|
|
481 |
goto L70; # End QZ step
|
|
482 |
|
|
483 |
L1000: # convergence failure
|
|
484 |
$bR->[$N-1][0] = $epsb;
|
|
485 |
return($en + 1);
|
|
486 |
|
|
487 |
L1001: # convergance okay
|
|
488 |
$bR->[$N-1][0] = $epsb;
|
|
489 |
return 0;
|
|
490 |
}
|
|
491 |
|
|
492 |
#----------------------------------------------------------------------
|
|
493 |
# QZval(\@A,\@B,\@Z,\@alphaR,\@alphaI,\@beta);
|
|
494 |
#----------------------------------------------------------------------
|
|
495 |
|
|
496 |
sub QZval($$$$$$)
|
|
497 |
{
|
|
498 |
my($aR,$bR,$zR,$alfrR,$alfiR,$betaR) = @_;
|
|
499 |
my($i,$j,$na,$en,$nn,$c,$d,$r,$s,$t,$di,$ei);
|
|
500 |
my($a1,$a2,$u1,$u2,$v1,$v2);
|
|
501 |
my($a11,$a12,$a21,$a22,$b11,$b12,$b22);
|
|
502 |
my($bn,$cq,$dr,$cz,$ti,$tr);
|
|
503 |
my($a1i,$a2i,$a11i,$a12i,$a22i,$a11r,$a12r,$a22r);
|
|
504 |
my($sqi,$ssi,$sqr,$szi,$ssr,$szr);
|
|
505 |
my($an,$e,$isw) = (0,0,1);
|
|
506 |
my($epsb) = $bR->[$N-1][0];
|
|
507 |
|
|
508 |
for ($nn=0; $nn<$N; $nn++) {
|
|
509 |
$en = $N - $nn - 1;
|
|
510 |
$na = $en - 1;
|
|
511 |
|
|
512 |
goto L505 if ($isw == 2);
|
|
513 |
goto L410 if ($en == 0);
|
|
514 |
goto L420 if ($aR->[$en][$na] != 0);
|
|
515 |
|
|
516 |
L410:
|
|
517 |
$alfrR->[$en] = ($bR->[$en][$en] < 0) ? -$alfrR->[$en] : $aR->[$en][$en];
|
|
518 |
$betaR->[$en] = (abs($bR->[$en][$en]));
|
|
519 |
$alfiR->[$en] = 0;
|
|
520 |
next;
|
|
521 |
|
|
522 |
L420:
|
|
523 |
goto L455 if (abs($bR->[$na][$na]) <= $epsb);
|
|
524 |
goto L430 if (abs($bR->[$en][$en]) > $epsb);
|
|
525 |
$a1 = $aR->[$en][$en];
|
|
526 |
$a2 = $aR->[$en][$na];
|
|
527 |
$bn = 0;
|
|
528 |
goto L435;
|
|
529 |
|
|
530 |
L430:
|
|
531 |
$an = abs($aR->[$na][$na]) + abs($aR->[$na][$en]) + abs($aR->[$en][$na]) + abs($aR->[$en][$en]);
|
|
532 |
$bn = abs($bR->[$na][$na]) + abs($bR->[$na][$en]) + abs($bR->[$en][$en]);
|
|
533 |
$a11 = $aR->[$na][$na] / $an;
|
|
534 |
$a12 = $aR->[$na][$en] / $an;
|
|
535 |
$a21 = $aR->[$en][$na] / $an;
|
|
536 |
$a22 = $aR->[$en][$en] / $an;
|
|
537 |
$b11 = $bR->[$na][$na] / $bn;
|
|
538 |
$b12 = $bR->[$na][$en] / $bn;
|
|
539 |
$b22 = $bR->[$en][$en] / $bn;
|
|
540 |
$e = $a11 / $b11;
|
|
541 |
$ei = $a22 / $b22;
|
|
542 |
$s = $a21 / ($b11 * $b22);
|
|
543 |
$t = ($a22 - $e * $b22) / $b22;
|
|
544 |
|
|
545 |
goto L431 if (abs($e) <= abs($ei));
|
|
546 |
$e = $ei;
|
|
547 |
$t = ($a11 - $e * $b11) / $b11;
|
|
548 |
|
|
549 |
L431:
|
|
550 |
$c = ($t - $s * $b12) / 2;
|
|
551 |
$d = $c**2 + $s * ($a12 - $e * $b12);
|
|
552 |
goto L480 if ($d < 0);
|
|
553 |
|
|
554 |
$e += $c + SIGN(sqrt($d),$c);
|
|
555 |
$a11 -= $e * $b11;
|
|
556 |
$a12 -= $e * $b12;
|
|
557 |
$a22 -= $e * $b22;
|
|
558 |
|
|
559 |
goto L432 if (abs($a11) + abs($a12) < abs($a21) + abs($a22));
|
|
560 |
|
|
561 |
$a1 = $a12;
|
|
562 |
$a2 = $a11;
|
|
563 |
goto L435;
|
|
564 |
|
|
565 |
L432:
|
|
566 |
$a1 = $a22;
|
|
567 |
$a2 = $a21;
|
|
568 |
|
|
569 |
L435:
|
|
570 |
$s = abs($a1) + abs($a2);
|
|
571 |
$u1 = $a1 / $s;
|
|
572 |
$u2 = $a2 / $s;
|
|
573 |
$r = SIGN(sqrt($u1**2 + $u2**2),$u1);
|
|
574 |
$v1 = -($u1 + $r) / $r;
|
|
575 |
$v2 = -$u2 / $r;
|
|
576 |
$u2 = $v2 / $v1;
|
|
577 |
|
|
578 |
for ($i=0; $i<=$en; $i++) {
|
|
579 |
$t = $aR->[$i][$en] + $u2 * $aR->[$i][$na];
|
|
580 |
$aR->[$i][$en] += $t * $v1;
|
|
581 |
$aR->[$i][$na] += $t * $v2;
|
|
582 |
|
|
583 |
$t = $bR->[$i][$e] + $u2 * $bR->[$i][$na];
|
|
584 |
$bR->[$i][$e] += $t * $v1;
|
|
585 |
$bR->[$i][$na] += $t * $v2;
|
|
586 |
}
|
|
587 |
|
|
588 |
for ($i=0; $i<$N; $i++) {
|
|
589 |
$t = $zR->[$i][$en] + $u2 * $zR->[$i][$na];
|
|
590 |
$zR->[$i][$en] += $t * $v1;
|
|
591 |
$zR->[$i][$na] += $t * $v2;
|
|
592 |
}
|
|
593 |
|
|
594 |
goto L475 if ($bn == 0);
|
|
595 |
goto L455 if ($an < abs($e) * $bn);
|
|
596 |
$a1 = $bR->[$na][$na];
|
|
597 |
$a2 = $bR->[$en][$na];
|
|
598 |
goto L460;
|
|
599 |
|
|
600 |
L455:
|
|
601 |
$a1 = $aR->[$na][$na];
|
|
602 |
$a2 = $aR->[$en][$na];
|
|
603 |
|
|
604 |
L460:
|
|
605 |
$s = abs($a1) + abs($a2);
|
|
606 |
goto L475 if ($s == 0);
|
|
607 |
$u1 = $a1 / $s;
|
|
608 |
$u2 = $a2 / $s;
|
|
609 |
$r = SIGN(sqrt($u1**2 + $u2**2),$u1);
|
|
610 |
$v1 = -($u1 + $r) / $r;
|
|
611 |
$v2 = -$u2 / $r;
|
|
612 |
$u2 = $v2 / $v1;
|
|
613 |
|
|
614 |
for ($j=$na; $j<$N; $j++) {
|
|
615 |
$t = $aR->[$na][$j] + $u2 * $aR->[$en][$j];
|
|
616 |
$aR->[$na][$j] += $t * $v1;
|
|
617 |
$aR->[$en][$j] += $t * $v2;
|
|
618 |
$t = $bR->[$na][$j] + $u2 * $bR->[$en][$j];
|
|
619 |
$bR->[$na][$j] += $t * $v1;
|
|
620 |
$bR->[$en][$j] += $t * $v2;
|
|
621 |
}
|
|
622 |
|
|
623 |
L475:
|
|
624 |
$aR->[$en][$na] = $bR->[$en][$na] = 0;
|
|
625 |
$alfrR->[$na] = $aR->[$na][$na];
|
|
626 |
$alfrR->[$en] = $aR->[$en][$en];
|
|
627 |
$alfrR->[$na] = -$alfrR->[$na]
|
|
628 |
if ($bR->[$na][$na] < 0);
|
|
629 |
$alfrR->[$en] = -$alfrR->[$en]
|
|
630 |
if ($bR->[$en][$en] < 0);
|
|
631 |
$betaR->[$na] = (abs($bR->[$na][$na]));
|
|
632 |
$betaR->[$en] = (abs($bR->[$en][$en]));
|
|
633 |
$alfiR->[$en] = $alfiR->[$na] = 0;
|
|
634 |
goto L505;
|
|
635 |
|
|
636 |
L480:
|
|
637 |
$e += $c;
|
|
638 |
$ei = sqrt(-$d);
|
|
639 |
$a11r = $a11 - $e * $b11;
|
|
640 |
$a11i = $ei * $b11;
|
|
641 |
$a12r = $a12 - $e * $b12;
|
|
642 |
$a12i = $ei * $b12;
|
|
643 |
$a22r = $a22 - $e * $b22;
|
|
644 |
$a22i = $ei * $b22;
|
|
645 |
|
|
646 |
goto L482
|
|
647 |
if (abs($a11r) + abs($a11i) + abs($a12r) + abs($a12i) <
|
|
648 |
abs($a21) + abs($a22r) + abs($a22i));
|
|
649 |
|
|
650 |
$a1 = $a12r; $a1i = $a12i;
|
|
651 |
$a2 = -$a11r; $a2i = -$a11i;
|
|
652 |
goto L485;
|
|
653 |
|
|
654 |
L482:
|
|
655 |
$a1 = $a22r;
|
|
656 |
$a1i = $a22i;
|
|
657 |
$a2 = -$a21;
|
|
658 |
$a2i = 0;
|
|
659 |
|
|
660 |
L485:
|
|
661 |
$cz = sqrt($a1**2 + $a1i**2);
|
|
662 |
goto L487 if ($cz == 0);
|
|
663 |
$szr = ($a1 * $a2 + $a1i * $a2i) / $cz;
|
|
664 |
$szi = ($a1 * $a2i - $a1i * $a2) / $cz;
|
|
665 |
$r = sqrt($cz**2 + $szr**2 + $szi**2);
|
|
666 |
$cz /= $r; $szr /= $r; $szi /= $r;
|
|
667 |
goto L490;
|
|
668 |
|
|
669 |
L487:
|
|
670 |
$szr = 1;
|
|
671 |
$szi = 0;
|
|
672 |
|
|
673 |
L490:
|
|
674 |
goto L492 if ($an < (abs($e) + $ei) * $bn);
|
|
675 |
$a1 = $cz * $b11 + $szr * $b12;
|
|
676 |
$a1i = $szi * $b12;
|
|
677 |
$a2 = $szr * $b22;
|
|
678 |
$a2i = $szi * $b22;
|
|
679 |
goto L495;
|
|
680 |
|
|
681 |
L492:
|
|
682 |
$a1 = $cz * $a11 + $szr * $a12;
|
|
683 |
$a1i = $szi * $a12;
|
|
684 |
$a2 = $cz * $a21 + $szr * $a22;
|
|
685 |
$a2i = $szi * $a22;
|
|
686 |
|
|
687 |
L495:
|
|
688 |
$cq = sqrt($a1**2 + $a1i**2);
|
|
689 |
goto L497 if ($cq == 0);
|
|
690 |
$sqr = ($a1 * $a2 + $a1i * $a2i) / $cq;
|
|
691 |
$sqi = ($a1 * $a2i - $a1i * $a2) / $cq;
|
|
692 |
$r = sqrt($cq**2 + $sqr**2 + $sqi**2);
|
|
693 |
$cq /= $r;
|
|
694 |
$sqr /= $r;
|
|
695 |
$sqi /= $r;
|
|
696 |
goto L500;
|
|
697 |
|
|
698 |
L497:
|
|
699 |
$sqr = 1;
|
|
700 |
$sqi = 0;
|
|
701 |
|
|
702 |
L500:
|
|
703 |
$ssr = $sqr * $szr + $sqi * $szi;
|
|
704 |
$ssi = $sqr * $szi - $sqi * $szr;
|
|
705 |
$i = 0;
|
|
706 |
$tr = $cq * $cz * $a11 + $cq * $szr * $a12 + $sqr * $cz * $a21 + $ssr * $a22;
|
|
707 |
$ti = $cq * $szi * $a12 - $sqi * $cz * $a21 + $ssi * $a22;
|
|
708 |
$dr = $cq * $cz * $b11 + $cq * $szr * $b12 + $ssr * $b22;
|
|
709 |
$di = $cq * $szi * $b12 + $ssi * $b22;
|
|
710 |
goto L503;
|
|
711 |
|
|
712 |
L502:
|
|
713 |
$i = 1;
|
|
714 |
$tr = $ssr * $a11 - $sqr * $cz * $a12 - $cq * $szr * $a21 + $cq * $cz * $a22;
|
|
715 |
$ti = -$ssi * $a11 - $sqi * $cz * $a12 + $cq * $szi * $a21;
|
|
716 |
$dr = $ssr * $b11 - $sqr * $cz * $b12 + $cq * $cz * $b22;
|
|
717 |
$di = -$ssi * $b11 - $sqi * $cz * $b12;
|
|
718 |
|
|
719 |
L503:
|
|
720 |
$t = $ti * $dr - $tr * $di;
|
|
721 |
$j = $na;
|
|
722 |
$j = $en if ($t < 0);
|
|
723 |
|
|
724 |
$r = sqrt($dr**2 + $di**2);
|
|
725 |
$betaR->[$j] = $bn * $r;
|
|
726 |
$alfrR->[$j] = $an * ($tr * $dr + $ti * $di) / $r;
|
|
727 |
$alfiR->[$j] = $an * $t / $r;
|
|
728 |
goto L502 if ($i == 0);
|
|
729 |
|
|
730 |
L505:
|
|
731 |
$isw = 3 - $isw;
|
|
732 |
|
|
733 |
} # main for $nn loop
|
|
734 |
|
|
735 |
$bR->[$N-1][0] = $epsb;
|
|
736 |
return 0;
|
|
737 |
}
|
|
738 |
|
|
739 |
#----------------------------------------------------------------------
|
|
740 |
# QZvec(\@A,\@B,\@Z,\@alphaR,\@alphaI,\@beta)
|
|
741 |
#----------------------------------------------------------------------
|
|
742 |
|
|
743 |
sub QZvec($$$$$$)
|
|
744 |
{
|
|
745 |
my($aR,$bR,$zR,$alfrR,$alfiR,$betaR) = @_;
|
|
746 |
my($i,$j,$k,$m,$na,$ii,$en,$jj,$nn,$enm2,$d,$q,$t,$w,$y,$t1,$t2,$w1,$di);
|
|
747 |
my($ra,$dr,$sa,$ti,$rr,$tr,$alfm,$almi,$betm,$almr);
|
|
748 |
my($r,$s,$x,$x1,$z1,$zz,$isw) = (0,0,0,0,0,0,1);
|
|
749 |
my($epsb) = $bR->[$N-1][0];
|
|
750 |
|
|
751 |
for ($nn=0; $nn<$N; $nn++) {
|
|
752 |
$en = $N - $nn - 1;
|
|
753 |
$na = $en - 1;
|
|
754 |
goto L795 if ($isw == 2);
|
|
755 |
goto L710 if ($alfiR->[$en] != 0);
|
|
756 |
|
|
757 |
$m = $en;
|
|
758 |
$bR->[$en][$en] = 1;
|
|
759 |
next if ($na == -1);
|
|
760 |
$alfm = $alfrR->[$m];
|
|
761 |
$betm = $betaR->[$m];
|
|
762 |
|
|
763 |
for ($ii=0; $ii<=$na; $ii++) {
|
|
764 |
$i = $en - $ii - 1;
|
|
765 |
$w = $betm * $aR->[$i][$i] - $alfm * $bR->[$i][$i];
|
|
766 |
$r = 0;
|
|
767 |
|
|
768 |
for ($j=$m; $j<=$en; $j++) {
|
|
769 |
$r += ($betm * $aR->[$i][$j] - $alfm * $bR->[$i][$j]) * $bR->[$j][$en];
|
|
770 |
}
|
|
771 |
|
|
772 |
goto L630 if ($i == 0 || $isw == 2);
|
|
773 |
goto L630 if ($betm * $aR->[$i,$i-1] == 0);
|
|
774 |
|
|
775 |
$zz = $w;
|
|
776 |
$s = $r;
|
|
777 |
goto L690;
|
|
778 |
|
|
779 |
L630:
|
|
780 |
$m = $i;
|
|
781 |
goto L640 if ($isw == 2);
|
|
782 |
|
|
783 |
$t = $w;
|
|
784 |
$t = $epsb if ($w == 0);
|
|
785 |
|
|
786 |
$bR->[$i][$en] = -$r / $t;
|
|
787 |
next;
|
|
788 |
|
|
789 |
L640:
|
|
790 |
$x = $betm * $aR->[$i][$i+1] - $alfm * $bR->[$i][$i+1];
|
|
791 |
$y = $betm * $aR->[$i+1][$i];
|
|
792 |
$q = $w * $zz - $x * $y;
|
|
793 |
$t = ($x * $s - $zz * $r) / $q;
|
|
794 |
$bR->[$i][$en] = $t;
|
|
795 |
goto L650 if (abs($x) <= abs($zz));
|
|
796 |
$bR->[$i+1][$en] = (-$r - $w * $t) / $x;
|
|
797 |
goto L690;
|
|
798 |
|
|
799 |
L650:
|
|
800 |
$bR->[$i+1][$en] = (-$s - $y * $t) / $zz;
|
|
801 |
|
|
802 |
L690:
|
|
803 |
$isw = 3 - $isw;
|
|
804 |
|
|
805 |
} # for ($ii inner loop
|
|
806 |
next;
|
|
807 |
|
|
808 |
L710:
|
|
809 |
$m = $na;
|
|
810 |
$almr = $alfrR->[$m];
|
|
811 |
$almi = $alfiR->[$m];
|
|
812 |
$betm = $betaR->[$m];
|
|
813 |
|
|
814 |
$y = $betm * $aR->[$en][$na];
|
|
815 |
$bR->[$na][$na] = -$almi * $bR->[$en][$en] / $y;
|
|
816 |
$bR->[$na][$en] = ($almr * $bR->[$en][$en] - $betm * $aR->[$en][$en]) / $y;
|
|
817 |
$bR->[$en][$na] = 0;
|
|
818 |
$bR->[$en][$en] = 1;
|
|
819 |
$enm2 = $na;
|
|
820 |
goto L795 if ($enm2 == 0);
|
|
821 |
|
|
822 |
for ($ii=0; $ii<$enm2; $ii++) {
|
|
823 |
$i = $na - $ii - 1;
|
|
824 |
$w = $betm * $aR->[$i][$i] - $almr * $bR->[$i][$i];
|
|
825 |
$w1 = -$almi * $bR->[$i][$i];
|
|
826 |
$ra = $sa = 0;
|
|
827 |
|
|
828 |
for ($j=$m; j<=$en; $j++) {
|
|
829 |
$x = $betm * $aR->[$i][$j] - $almr * $bR->[$i][$j];
|
|
830 |
$x1 = -$almi * $bR->[$i][$j];
|
|
831 |
$ra = $ra + $x * $bR->[$j][$na] - $x1 * $bR->[$j][$en];
|
|
832 |
$sa = $sa + $x * $bR->[$j][$en] + $x1 * $bR->[$j][$na];
|
|
833 |
}
|
|
834 |
|
|
835 |
goto L770 if ($i == 0 || $isw == 2);
|
|
836 |
goto L770 if ($betm * $aR->[$i,$i-1] == 0);
|
|
837 |
|
|
838 |
$zz = $w; $z1 = $w1;
|
|
839 |
$r = $ra; $s = $sa;
|
|
840 |
$isw = 2;
|
|
841 |
next;
|
|
842 |
|
|
843 |
L770:
|
|
844 |
$m = $i;
|
|
845 |
goto L780 if ($isw == 2);
|
|
846 |
$tr = -$ra; $ti = -$sa;
|
|
847 |
|
|
848 |
L773:
|
|
849 |
$dr = $w; $di = $w1;
|
|
850 |
|
|
851 |
L775:
|
|
852 |
goto L777 if (abs($di) > abs($dr));
|
|
853 |
$rr = $di / $dr;
|
|
854 |
$d = $dr + $di * $rr;
|
|
855 |
$t1 = ($tr + $ti * $rr) / $d;
|
|
856 |
$t2 = ($ti - $tr * $rr) / $d;
|
|
857 |
if ($isw == 1) { goto L787; }
|
|
858 |
elsif ($isw == 2) { goto L782; }
|
|
859 |
|
|
860 |
L777:
|
|
861 |
$rr = $dr / $di;
|
|
862 |
$d = $dr * $rr + $di;
|
|
863 |
$t1 = ($tr * $rr + $ti) / $d;
|
|
864 |
$t2 = ($ti * $rr - $tr) / $d;
|
|
865 |
if ($isw == 1) { goto L787; }
|
|
866 |
elsif ($isw == 2) { goto L782; }
|
|
867 |
|
|
868 |
L780:
|
|
869 |
$x = $betm * $aR->[$i][$i+1] - $almr * $bR->[$i][$i+1];
|
|
870 |
$x1 = -$almi * $bR->[$i][$i+1];
|
|
871 |
$y = $betm * $aR->[$i+1][$i];
|
|
872 |
$tr = $y * $ra - $w * $r + $w1 * $s;
|
|
873 |
$ti = $y * $sa - $w * $s - $w1 * $r;
|
|
874 |
$dr = $w * $zz - $w1 * $z1 - $x * $y;
|
|
875 |
$di = $w * $z1 + $w1 * $zz - $x1 * $y;
|
|
876 |
$dr = $epsb if ($dr == 0 && $di == 0);
|
|
877 |
goto L775;
|
|
878 |
|
|
879 |
L782:
|
|
880 |
$bR->[$i+1][$na] = $t1;
|
|
881 |
$bR->[$i+1][$en] = $t2;
|
|
882 |
$isw = 1;
|
|
883 |
goto L785 if (abs($y) > abs($w) + abs($w1));
|
|
884 |
$tr = -$ra - $x * $bR->[$i+1][$na] + $x1 * $bR->[$i+1][$en];
|
|
885 |
$ti = -$sa - $x * $bR->[$i+1][$en] - $x1 * $bR->[$i+1][$na];
|
|
886 |
goto L773;
|
|
887 |
|
|
888 |
L785:
|
|
889 |
$t1 = (-$r - $zz * $bR->[$i+1][$na] + $z1 * $bR->[$i+1][$en]) / $y;
|
|
890 |
$t2 = (-$s - $zz * $bR->[$i+1][$en] - $z1 * $bR->[$i+1][$na]) / $y;
|
|
891 |
|
|
892 |
L787:
|
|
893 |
$bR->[$i][$na] = $t1;
|
|
894 |
$bR->[$i][$en] = $t2;
|
|
895 |
|
|
896 |
} # for ($ii inner loop
|
|
897 |
|
|
898 |
L795:
|
|
899 |
$isw = 3 - $isw;
|
|
900 |
|
|
901 |
} # for ($nn outer loop
|
|
902 |
|
|
903 |
for ($jj=0; $jj<$N; $jj++) {
|
|
904 |
$j = $N - $jj - 1;
|
|
905 |
for ($i=0; $i<$N; $i++) {
|
|
906 |
$zz = 0;
|
|
907 |
for ($k=0; $k<=$j; $k++) {
|
|
908 |
$zz += $zR->[$i][$k] * $bR->[$k][$j];
|
|
909 |
}
|
|
910 |
$zR->[$i][$j] = $zz;
|
|
911 |
}
|
|
912 |
}
|
|
913 |
|
|
914 |
for ($j=0; $j<$N; $j++) {
|
|
915 |
$d = 0;
|
|
916 |
goto L920 if ($isw == 2);
|
|
917 |
goto L945 if ($alfiR->[$j] != 0);
|
|
918 |
for ($i=0; $i<$N; $i++) {
|
|
919 |
$d = (abs($zR->[$i][$j])) if ((abs($zR->[$i][$j])) > $d);
|
|
920 |
}
|
|
921 |
for ($i=0; $i<$N; $i++) {
|
|
922 |
$zR->[$i][$j] /= $d;
|
|
923 |
}
|
|
924 |
next;
|
|
925 |
|
|
926 |
L920:
|
|
927 |
for ($i=0; $i<$N; $i++) {
|
|
928 |
$r = abs($zR->[$i][$j-1]) + abs($zR->[$i][$j]);
|
|
929 |
if ($r != 0) {
|
|
930 |
my($u1) = $zR->[$i][$j-1] / $r;
|
|
931 |
my($u2) = $zR->[$i][$j] / $r;
|
|
932 |
$r *= sqrt($u1**2 + $u2**2);
|
|
933 |
}
|
|
934 |
$d = $r if ($r > $d);
|
|
935 |
}
|
|
936 |
|
|
937 |
for ($i=0; $i<$N; $i++) {
|
|
938 |
$zR->[$i][$j-1] /= $d;
|
|
939 |
$zR->[$i][$j] /= $d;
|
|
940 |
}
|
|
941 |
|
|
942 |
L945:
|
|
943 |
$isw = 3 - $isw;
|
|
944 |
|
|
945 |
} # for ($j outer loop
|
|
946 |
}
|
|
947 |
|
|
948 |
1;
|
|
949 |
|