diff --git a/libQZ.pl b/libQZ.pl new file mode 100644 --- /dev/null +++ b/libQZ.pl @@ -0,0 +1,949 @@ +#====================================================================== +# L I B Q Z . P L +# doc: Thu Mar 12 15:23:15 2015 +# dlm: Sun Mar 15 19:26:50 2015 +# (c) 2015 A.M. Thurnherr +# uE-Info: 14 27 NIL 0 0 72 10 2 4 NIL ofnI +#====================================================================== + +# adaptation of EISPACK routines +# www.netlib.org/eispack + +# HISTORY: +# Mar 12, 2015: - created +# Mar 15, 2015: - debugging + +use strict vars; + +my($N); # size of input matrices A & B + +#---------------------------------------------------------------------- +# eig(\@A,\@B,\@evRe,\@evIm,\@V) +# - @ev{Re,Im} contain generalized eigenvalues +# - @evec contains corresponding right eigenvectors +# - A*V = B*V*D, where D is diagonal matrix of eigenvalues +#---------------------------------------------------------------------- + +sub eig($$$$$) +{ + my($aR,$bR,$erR,$eiR,$zR) = @_; # args passed as refs + my(@alphaR,@alphaI,@beta); # intermediate data + + $N = scalar(@{$aR}); + croak(sprintf("eig(A,B): A(%dx%d) & B(%dx%d) must be matching square matrices\n", + $N,scalar(@{$aR->[0]}),scalar(@{$bR}),scalar(@{$bR->[0]}))) + unless (@{$bR} == $N) && (@{$aR->[0]} == $N) && (@{$bR->[0]} == $N); + + QZhes($aR,$bR,$zR); # reduce A/B to upper Hessenberg/triangular forms + croak("QZit(): convergence failed\n") + unless (QZit($aR,$bR,$zR) == 0); # reduce Hess A to quasi-triangular form + QZval($aR,$bR,$zR,\@alphaR,\@alphaI,\@beta); # reduce A further + QZvec($aR,$bR,$zR,\@alphaR,\@alphaI,\@beta); # compute eigenvectors & eigenvalues + + for (my($i)=0; $i<$N; $i++) { + if ($beta[$i]==0 && $alphaR[$i]==0) { + $erR->[$i] = nan; + $eiR->[$i] = nan; + } else { + $erR->[$i] = $alphaR[$i] / $beta[$i]; + $eiR->[$i] = $alphaI[$i] / $beta[$i]; +# print(STDERR "gev[$i] = $erR->[$i]\n"); + } + } +} + +#---------------------------------------------------------------------- +# EISPACK routines +#---------------------------------------------------------------------- + +#---------------------------------------------------------------------- +# QZhes(\@A,\@B,\@Z) +# - first step in QZ algorithm (Moler & Stewart, SIAM JNA, 1973) +#---------------------------------------------------------------------- + +sub QZhes($$$) +{ + my($aR,$bR,$zR) = @_; + my($i,$j,$k,$l,$l1,$lb,$nk1); + my($r,$s,$t,$u1,$u2,$v1,$v2,$rho); + + croak("QZhes: need at least 3x3 matrices\n") + unless ($N >= 2); + + for ($j=0; $j<$N; $j++) { # init Z to identity matrix + for ($i=0; $i<$N; $i++) { + $zR->[$i][$j] = 0; + } + $zR->[$j][$j] = 1; + } + + for ($l=0; $l<$N-1; $l++) { + $l1 = $l + 1; + for ($s=0,$i=$l1; $i<$N; $i++) { + $s += abs($bR->[$i][$l]); + } + next if ($s == 0); + $s += abs($bR->[$l][$l]); + + for ($r=0,$i=$l; $i<$N; $i++) + { + $bR->[$i][$l] /= $s; + $r += $bR->[$i][$l]**2; + } + + $r = SIGN(sqrt($r),$bR->[$l][$l]); + $bR->[$l][$l] += $r; + $rho = $r * $bR->[$l][$l]; + + for ($j=$l1; $j<$N; $j++) { + for ($t=0,$i=$l; $i<$N; $i++) { + $t += $bR->[$i][$l] * $bR->[$i][$j]; + } + $t = -$t / $rho; + for ($i=$l; $i<$N; $i++) { + $bR->[$i][$j] += $t * $bR->[$i][$l]; + } + } + + for ($j=0; $j<$N; $j++) { + for ($t=0,$i=$l; $i<$N; $i++) { + $t += $bR->[$i][$l] * $aR->[$i][$j]; + } + $t = -$t / $rho; + for ($i=$l; $i<$N; $i++) { + $aR->[$i][$j] += $t * $bR->[$i][$l]; + } + } + + $bR->[$l][$l] = -$s * $r; + for ($i=$l1; $i<$N; $i++) { + $bR->[$i][$l] = 0; + } + } + + for ($k=0; $k<$N-2; $k++) { + $nk1 = $N - 2 - $k; + + for ($lb=0; $lb<$nk1; $lb++) { + $l = $N - $lb - 2; + $l1 = $l + 1; + + $s = (abs($aR->[$l][$k])) + (abs($aR->[$l1][$k])); + next if ($s == 0); + + $u1 = $aR->[$l][$k] / $s; + $u2 = $aR->[$l1][$k] / $s; + $r = SIGN(sqrt($u1**2+$u2**2),$u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($j=$k; $j<$N; $j++) { + $t = $aR->[$l][$j] + $u2 * $aR->[$l1][$j]; + $aR->[$l][$j] += $t * $v1; + $aR->[$l1][$j] += $t * $v2; + } + + $aR->[$l1][$k] = 0; + + for ($j=$l; $j<$N; $j++) { + $t = $bR->[$l][$j] + $u2 * $bR->[$l1][$j]; + $bR->[$l][$j] += $t * $v1; + $bR->[$l1][$j] += $t * $v2; + } + + $s = (abs($bR->[$l1][$l1])) + (abs($bR->[$l1][$l])); + next if ($s == 0); + + $u1 = $bR->[$l1][$l1] / $s; + $u2 = $bR->[$l1][$l] / $s; + $r = SIGN(sqrt($u1**2 + $u2**2),$u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($i=0; $i<=$l1; $i++) { # overwrite B with upper triangular form + $t = $bR->[$i][$l1] + $u2 * $bR->[$i][$l]; + $bR->[$i][$l1] += $t * $v1; + $bR->[$i][$l] += $t * $v2; + } + $bR->[$l1][$l] = 0; + + for ($i=0; $i<$N; $i++) { # overwrite A with upper Hessenberg form + $t = $aR->[$i][$l1] + $u2 * $aR->[$i][$l]; + $aR->[$i][$l1] += $t * $v1; + $aR->[$i][$l] += $t * $v2; + } + + for ($i=0; $i<$N; $i++) { # define matrix Z, used for eigenvectors + $t = $zR->[$i][$l1] + $u2 * $zR->[$i][$l]; + $zR->[$i][$l1] += $t * $v1; + $zR->[$i][$l] += $t * $v2; + } + } # for ($lb=0; $lb<$nk1; $lb++) + } # for ($k=0; $k<$N-2; $k++) +} + +#---------------------------------------------------------------------- +# QZit(\@A,\@B,\@Z) +# - second step in QZ algorithm (Moler & Stewart, SIAM JNA, 1973) +# - !0 return value indicates that convergence failed +#---------------------------------------------------------------------- + +sub QZit($$$$) +{ + my($aR,$bR,$zR) = @_; + + my($i,$j,$k); + my($r,$s,$t,$a1,$a2); + my($k1,$k2,$l1,$ll); + my($u1,$u2,$u3); + my($v1,$v2,$v3); + my($a11,$a12,$a21,$a22,$a33,$a34,$a43,$a44); + my($b11,$b12,$b22,$b33,$b34,$b44); + my($na,$en,$ld); + my($ep,$km1,$ani,$bni); + my($ish,$itn,$its,$enm2,$lor1,$epsa,$epsb,$enorn,$notlas); + my($l,$a3,$sh,$lm1,$anorm,$bnorm) = (0,0,0,0,0,0); + + for ($i=0; $i<$N; $i++) { + $ani = $bni = 0; + $ani = abs($aR->[$i][$i-1]) + unless ($i == 0); + for ($j=$i; $j<$N; $j++) { + $ani += abs($aR->[$i][$j]); + $bni += abs($bR->[$i][$j]); + } + $anorm = $ani if ($ani > $anorm); + $bnorm = $bni if ($bni > $bnorm); + } + $anorm = 1 if ($anorm == 0); + $bnorm = 1 if ($bnorm == 0); + + $ep = 1.11022302462516e-16; # $EPS=1; $EPS /= 2 while 0.5 + $EPS/2 > 0.5; + $epsa = $ep * $anorm; + $epsb = $ep * $bnorm; + + $lor1 = 0; + $enorn = $N; + $en = $N - 1; + $itn = $N * 30; + +L60: + goto L1001 if ($en <= 1); + $its = 0; + $na = $en - 1; + $enm2 = $na; + +L70: + $ish = 2; + for ($ll=0; $ll<=$en; $ll++) { + $lm1 = $en - $ll - 1; + $l = $lm1 + 1; + goto L95 if ($l == 0); + last if ((abs($aR->[$l][$lm1])) <= $epsa) + } + +L90: + $aR->[$l][$lm1] = 0; + goto L95 if ($l < $na); + $en = $lm1; + goto L60; + +L95: + $ld = $l; + +L100: + $l1 = $l + 1; + $b11 = $bR->[$l][$l]; + goto L120 if (abs($b11) > $epsb); + $bR->[$l][$l] = 0; + $s = (abs($aR->[$l][$l]) + abs($aR->[$l1][$l])); + $u1 = $aR->[$l][$l] / $s; + $u2 = $aR->[$l1][$l] / $s; + $r = SIGN(sqrt($u1**2 + $u2**2),$u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($j=$l; $j<$enorn; $j++) { + $t = $aR->[$l][$j] + $u2 * $aR->[$l1][$j]; + $aR->[$l][$j] += $t * $v1; + $aR->[$l1][$j] += $t * $v2; + + $t = $bR->[$l][$j] + $u2 * $bR->[$l1][$j]; + $bR->[$l][$j] += $t * $v1; + $bR->[$l1][$j] += $t * $v2; + } + $aR->[$l][$lm1] = -$aR->[$l][$lm1] + if ($l != 0); + + $lm1 = $l; + $l = $l1; + goto L90; + +L120: + $a11 = $aR->[$l][$l] / $b11; + $a21 = $aR->[$l1][$l] / $b11; + + goto L140 if ($ish == 1); + goto L1000 if ($itn == 0); + goto L155 if ($its == 10); + + $b22 = $bR->[$l1][$l1]; + $b22 = $epsb if (abs($b22) < $epsb); + $b33 = $bR->[$na][$na]; + $b33 = $epsb if (abs($b33) < $epsb); + $b44 = $bR->[$en][$en]; + $b44 = $epsb if (abs($b44) < $epsb); + + $a33 = $aR->[$na][$na] / $b33; + $a34 = $aR->[$na][$en] / $b44; + $a43 = $aR->[$en][$na] / $b33; + $a44 = $aR->[$en][$en] / $b44; + $b34 = $bR->[$na][$en] / $b44; + + $t = ($a43 * $b34 - $a33 - $a44) / 2; + $r = $t * $t + $a34 * $a43 - $a33 * $a44; + goto L150 if ($r < 0); + + $ish = 1; + $r = sqrt($r); + $sh = -$t + $r; + $s = -$t - $r; + $sh = $s if (abs($s-$a44) < abs($sh-$a44)); + + for ($ll=$ld; $ll<$enm2; $ll++) { + $l = $enm2 + $ld - $ll - 1; + goto L140 if ($l == $ld); + $lm1 = $l - 1; + $l1 = $l + 1; + $t = $aR->[$l+1][$l+1]; + $t -= $sh * $bR->[$l][$l] if (abs($bR->[$l][$l]) > $epsb); + goto L100 if (abs($aR->[$l][$lm1]) <= (abs($t / $aR->[$l1][$l])) * $epsa); + } + +L140: + $a1 = $a11 - $sh; + $a2 = $a21; + $aR->[$l][$lm1] = -$aR->[$l][$lm1] if ($l != $ld); + goto L160; + +L150: + $a12 = $aR->[$l][$l1] / $b22; + $a22 = $aR->[$l1][$l1] / $b22; + $b12 = $bR->[$l][$l1] / $b22; + $a1 = (($a33 - $a11) * ($a44 - $a11) - $a34 * $a43 + $a43 * $b34 * $a11) / $a21 + $a12 - $a11 * $b12; + $a2 = $a22 - $a11 - $a21 * $b12 - ($a33 - $a11) - ($a44 - $a11) + $a43 * $b34; + $a3 = $aR->[$l1+1][$l] / $b22; + goto L160; + +L155: + $a1 = 0; + $a2 = 1; + $a3 = 1.1605; + +L160: + $its++; + $itn--; + $lor1 = $ld; + for ($k=$l; $k<=$na; $k++) { + + $notlas = ($k!=$na) && ($ish==2); + $k1 = $k + 1; + $k2 = $k + 2; + $km1 = MAX($k,$l+1)-1; + $ll = MIN($en,$k1+$ish); + goto L190 if ($notlas); + + if ($k != $l) { + $a1 = $aR->[$k,$km1]; + $a2 = $aR->[$k1,$km1]; + } + + $s = abs($a1) + abs($a2); + goto L70 if ($s == 0); + $u1 = $a1 / $s; + $u2 = $a2 / $s; + $r = SIGN(sqrt($u1**2 + $u2**2), $u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($j=$km1; $j<$enorn; $j++) { + $t = $aR->[$k][$j] + $u2 * $aR->[$k1][$j]; + $aR->[$k][$j] += $t * $v1; + $aR->[$k1][$j] += $t * $v2; + $t = $bR->[$k][$j] + $u2 * $bR->[$k1][$j]; + $bR->[$k][$j] += $t * $v1; + $bR->[$k1][$j] += $t * $v2; + } + + $aR->[$k1,$km1] = 0 if ($k != $l); + goto L240; + + L190: + goto L200 if ($k == $l); + $a1 = $aR->[$k,$km1]; + $a2 = $aR->[$k1,$km1]; + $a3 = $aR->[$k2][$km1]; + + L200: + $s = abs($a1) + abs($a2) + abs($a3); + next if ($s == 0); + $u1 = $a1 / $s; + $u2 = $a2 / $s; + $u3 = $a3 / $s; + $r = SIGN(sqrt($u1**2 + $u2**2 + $u3**2), $u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $v3 = -$u3 / $r; + $u2 = $v2 / $v1; + $u3 = $v3 / $v1; + + for ($j=$km1; $j<$enorn; $j++) { + $t = $aR->[$k][$j] + $u2 * $aR->[$k1][$j] + $u3 * $aR->[$k2][$j]; + $aR->[$k][$j] += $t * $v1; + $aR->[$k1][$j] += $t * $v2; + $aR->[$k2][$j] += $t * $v3; + + $t = $bR->[$k][$j] + $u2 * $bR->[$k1][$j] + $u3 * $bR->[$k2][$j]; + $bR->[$k][$j] += $t * $v1; + $bR->[$k1][$j] += $t * $v2; + $bR->[$k2][$j] += $t * $v3; + } + + goto L220 if ($k == $l); + $aR->[$k1,$km1] = $aR->[$k2][$km1] = 0; + + L220: + $s = (abs($bR->[$k2][$k2])) + (abs($bR->[$k2][$k1])) + (abs($bR->[$k2][$k])); + goto L240 if ($s == 0); + $u1 = $bR->[$k2][$k2] / $s; + $u2 = $bR->[$k2][$k1] / $s; + $u3 = $bR->[$k2][$k] / $s; + $r = SIGN(sqrt($u1**2 + $u2**2 + $u3**2), $u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $v3 = -$u3 / $r; + $u2 = $v2 / $v1; + $u3 = $v3 / $v1; + + for ($i=$lor1; $i<$ll+1; $i++) { + $t = $aR->[$i][$k2] + $u2 * $aR->[$i][$k1] + $u3 * $aR->[$i][$k]; + $aR->[$i][$k2] += $t * $v1; + $aR->[$i][$k1] += $t * $v2; + $aR->[$i][$k] += $t * $v3; + $t = $bR->[$i][$k2] + $u2 * $bR->[$i][$k1] + $u3 * $bR->[$i][$k]; + $bR->[$i][$k2] += $t * $v1; + $bR->[$i][$k1] += $t * $v2; + $bR->[$i][$k] += $t * $v3; + } + + $bR->[$k2][$k] = $bR->[$k2][$k1] = 0; + + for ($i=0; $i<$N; $i++) { + $t = $zR->[$i][$k2] + $u2 * $zR->[$i][$k1] + $u3 * $zR->[$i][$k]; + $zR->[$i][$k2] += $t * $v1; + $zR->[$i][$k1] += $t * $v2; + $zR->[$i][$k] += $t * $v3; + } + + L240: + $s = (abs($bR->[$k1][$k1])) + (abs($bR->[$k1][$k])); + next if ($s == 0); + $u1 = $bR->[$k1][$k1] / $s; + $u2 = $bR->[$k1][$k] / $s; + $r = SIGN(sqrt($u1**2 + $u2**2), $u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($i=$lor1; $i<$ll+1; $i++) { + $t = $aR->[$i][$k1] + $u2 * $aR->[$i][$k]; + $aR->[$i][$k1] += $t * $v1; + $aR->[$i][$k] += $t * $v2; + $t = $bR->[$i][$k1] + $u2 * $bR->[$i][$k]; + $bR->[$i][$k1] += $t * $v1; + $bR->[$i][$k] += $t * $v2; + } + + $bR->[$k1][$k] = 0; + + for ($i=0; $i<$N; $i++) { + $t = $zR->[$i][$k1] + $u2 * $zR->[$i][$k]; + $zR->[$i][$k1] += $t * $v1; + $zR->[$i][$k] += $t * $v2; + } + + } # for loop beginning at L160 + + goto L70; # End QZ step + +L1000: # convergence failure + $bR->[$N-1][0] = $epsb; + return($en + 1); + +L1001: # convergance okay + $bR->[$N-1][0] = $epsb; + return 0; +} + +#---------------------------------------------------------------------- +# QZval(\@A,\@B,\@Z,\@alphaR,\@alphaI,\@beta); +#---------------------------------------------------------------------- + +sub QZval($$$$$$) +{ + my($aR,$bR,$zR,$alfrR,$alfiR,$betaR) = @_; + my($i,$j,$na,$en,$nn,$c,$d,$r,$s,$t,$di,$ei); + my($a1,$a2,$u1,$u2,$v1,$v2); + my($a11,$a12,$a21,$a22,$b11,$b12,$b22); + my($bn,$cq,$dr,$cz,$ti,$tr); + my($a1i,$a2i,$a11i,$a12i,$a22i,$a11r,$a12r,$a22r); + my($sqi,$ssi,$sqr,$szi,$ssr,$szr); + my($an,$e,$isw) = (0,0,1); + my($epsb) = $bR->[$N-1][0]; + + for ($nn=0; $nn<$N; $nn++) { + $en = $N - $nn - 1; + $na = $en - 1; + + goto L505 if ($isw == 2); + goto L410 if ($en == 0); + goto L420 if ($aR->[$en][$na] != 0); + + L410: + $alfrR->[$en] = ($bR->[$en][$en] < 0) ? -$alfrR->[$en] : $aR->[$en][$en]; + $betaR->[$en] = (abs($bR->[$en][$en])); + $alfiR->[$en] = 0; + next; + + L420: + goto L455 if (abs($bR->[$na][$na]) <= $epsb); + goto L430 if (abs($bR->[$en][$en]) > $epsb); + $a1 = $aR->[$en][$en]; + $a2 = $aR->[$en][$na]; + $bn = 0; + goto L435; + + L430: + $an = abs($aR->[$na][$na]) + abs($aR->[$na][$en]) + abs($aR->[$en][$na]) + abs($aR->[$en][$en]); + $bn = abs($bR->[$na][$na]) + abs($bR->[$na][$en]) + abs($bR->[$en][$en]); + $a11 = $aR->[$na][$na] / $an; + $a12 = $aR->[$na][$en] / $an; + $a21 = $aR->[$en][$na] / $an; + $a22 = $aR->[$en][$en] / $an; + $b11 = $bR->[$na][$na] / $bn; + $b12 = $bR->[$na][$en] / $bn; + $b22 = $bR->[$en][$en] / $bn; + $e = $a11 / $b11; + $ei = $a22 / $b22; + $s = $a21 / ($b11 * $b22); + $t = ($a22 - $e * $b22) / $b22; + + goto L431 if (abs($e) <= abs($ei)); + $e = $ei; + $t = ($a11 - $e * $b11) / $b11; + + L431: + $c = ($t - $s * $b12) / 2; + $d = $c**2 + $s * ($a12 - $e * $b12); + goto L480 if ($d < 0); + + $e += $c + SIGN(sqrt($d),$c); + $a11 -= $e * $b11; + $a12 -= $e * $b12; + $a22 -= $e * $b22; + + goto L432 if (abs($a11) + abs($a12) < abs($a21) + abs($a22)); + + $a1 = $a12; + $a2 = $a11; + goto L435; + + L432: + $a1 = $a22; + $a2 = $a21; + + L435: + $s = abs($a1) + abs($a2); + $u1 = $a1 / $s; + $u2 = $a2 / $s; + $r = SIGN(sqrt($u1**2 + $u2**2),$u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($i=0; $i<=$en; $i++) { + $t = $aR->[$i][$en] + $u2 * $aR->[$i][$na]; + $aR->[$i][$en] += $t * $v1; + $aR->[$i][$na] += $t * $v2; + + $t = $bR->[$i][$e] + $u2 * $bR->[$i][$na]; + $bR->[$i][$e] += $t * $v1; + $bR->[$i][$na] += $t * $v2; + } + + for ($i=0; $i<$N; $i++) { + $t = $zR->[$i][$en] + $u2 * $zR->[$i][$na]; + $zR->[$i][$en] += $t * $v1; + $zR->[$i][$na] += $t * $v2; + } + + goto L475 if ($bn == 0); + goto L455 if ($an < abs($e) * $bn); + $a1 = $bR->[$na][$na]; + $a2 = $bR->[$en][$na]; + goto L460; + + L455: + $a1 = $aR->[$na][$na]; + $a2 = $aR->[$en][$na]; + + L460: + $s = abs($a1) + abs($a2); + goto L475 if ($s == 0); + $u1 = $a1 / $s; + $u2 = $a2 / $s; + $r = SIGN(sqrt($u1**2 + $u2**2),$u1); + $v1 = -($u1 + $r) / $r; + $v2 = -$u2 / $r; + $u2 = $v2 / $v1; + + for ($j=$na; $j<$N; $j++) { + $t = $aR->[$na][$j] + $u2 * $aR->[$en][$j]; + $aR->[$na][$j] += $t * $v1; + $aR->[$en][$j] += $t * $v2; + $t = $bR->[$na][$j] + $u2 * $bR->[$en][$j]; + $bR->[$na][$j] += $t * $v1; + $bR->[$en][$j] += $t * $v2; + } + + L475: + $aR->[$en][$na] = $bR->[$en][$na] = 0; + $alfrR->[$na] = $aR->[$na][$na]; + $alfrR->[$en] = $aR->[$en][$en]; + $alfrR->[$na] = -$alfrR->[$na] + if ($bR->[$na][$na] < 0); + $alfrR->[$en] = -$alfrR->[$en] + if ($bR->[$en][$en] < 0); + $betaR->[$na] = (abs($bR->[$na][$na])); + $betaR->[$en] = (abs($bR->[$en][$en])); + $alfiR->[$en] = $alfiR->[$na] = 0; + goto L505; + + L480: + $e += $c; + $ei = sqrt(-$d); + $a11r = $a11 - $e * $b11; + $a11i = $ei * $b11; + $a12r = $a12 - $e * $b12; + $a12i = $ei * $b12; + $a22r = $a22 - $e * $b22; + $a22i = $ei * $b22; + + goto L482 + if (abs($a11r) + abs($a11i) + abs($a12r) + abs($a12i) < + abs($a21) + abs($a22r) + abs($a22i)); + + $a1 = $a12r; $a1i = $a12i; + $a2 = -$a11r; $a2i = -$a11i; + goto L485; + + L482: + $a1 = $a22r; + $a1i = $a22i; + $a2 = -$a21; + $a2i = 0; + + L485: + $cz = sqrt($a1**2 + $a1i**2); + goto L487 if ($cz == 0); + $szr = ($a1 * $a2 + $a1i * $a2i) / $cz; + $szi = ($a1 * $a2i - $a1i * $a2) / $cz; + $r = sqrt($cz**2 + $szr**2 + $szi**2); + $cz /= $r; $szr /= $r; $szi /= $r; + goto L490; + + L487: + $szr = 1; + $szi = 0; + + L490: + goto L492 if ($an < (abs($e) + $ei) * $bn); + $a1 = $cz * $b11 + $szr * $b12; + $a1i = $szi * $b12; + $a2 = $szr * $b22; + $a2i = $szi * $b22; + goto L495; + + L492: + $a1 = $cz * $a11 + $szr * $a12; + $a1i = $szi * $a12; + $a2 = $cz * $a21 + $szr * $a22; + $a2i = $szi * $a22; + + L495: + $cq = sqrt($a1**2 + $a1i**2); + goto L497 if ($cq == 0); + $sqr = ($a1 * $a2 + $a1i * $a2i) / $cq; + $sqi = ($a1 * $a2i - $a1i * $a2) / $cq; + $r = sqrt($cq**2 + $sqr**2 + $sqi**2); + $cq /= $r; + $sqr /= $r; + $sqi /= $r; + goto L500; + + L497: + $sqr = 1; + $sqi = 0; + + L500: + $ssr = $sqr * $szr + $sqi * $szi; + $ssi = $sqr * $szi - $sqi * $szr; + $i = 0; + $tr = $cq * $cz * $a11 + $cq * $szr * $a12 + $sqr * $cz * $a21 + $ssr * $a22; + $ti = $cq * $szi * $a12 - $sqi * $cz * $a21 + $ssi * $a22; + $dr = $cq * $cz * $b11 + $cq * $szr * $b12 + $ssr * $b22; + $di = $cq * $szi * $b12 + $ssi * $b22; + goto L503; + + L502: + $i = 1; + $tr = $ssr * $a11 - $sqr * $cz * $a12 - $cq * $szr * $a21 + $cq * $cz * $a22; + $ti = -$ssi * $a11 - $sqi * $cz * $a12 + $cq * $szi * $a21; + $dr = $ssr * $b11 - $sqr * $cz * $b12 + $cq * $cz * $b22; + $di = -$ssi * $b11 - $sqi * $cz * $b12; + + L503: + $t = $ti * $dr - $tr * $di; + $j = $na; + $j = $en if ($t < 0); + + $r = sqrt($dr**2 + $di**2); + $betaR->[$j] = $bn * $r; + $alfrR->[$j] = $an * ($tr * $dr + $ti * $di) / $r; + $alfiR->[$j] = $an * $t / $r; + goto L502 if ($i == 0); + + L505: + $isw = 3 - $isw; + + } # main for $nn loop + + $bR->[$N-1][0] = $epsb; + return 0; +} + +#---------------------------------------------------------------------- +# QZvec(\@A,\@B,\@Z,\@alphaR,\@alphaI,\@beta) +#---------------------------------------------------------------------- + +sub QZvec($$$$$$) +{ + my($aR,$bR,$zR,$alfrR,$alfiR,$betaR) = @_; + my($i,$j,$k,$m,$na,$ii,$en,$jj,$nn,$enm2,$d,$q,$t,$w,$y,$t1,$t2,$w1,$di); + my($ra,$dr,$sa,$ti,$rr,$tr,$alfm,$almi,$betm,$almr); + my($r,$s,$x,$x1,$z1,$zz,$isw) = (0,0,0,0,0,0,1); + my($epsb) = $bR->[$N-1][0]; + + for ($nn=0; $nn<$N; $nn++) { + $en = $N - $nn - 1; + $na = $en - 1; + goto L795 if ($isw == 2); + goto L710 if ($alfiR->[$en] != 0); + + $m = $en; + $bR->[$en][$en] = 1; + next if ($na == -1); + $alfm = $alfrR->[$m]; + $betm = $betaR->[$m]; + + for ($ii=0; $ii<=$na; $ii++) { + $i = $en - $ii - 1; + $w = $betm * $aR->[$i][$i] - $alfm * $bR->[$i][$i]; + $r = 0; + + for ($j=$m; $j<=$en; $j++) { + $r += ($betm * $aR->[$i][$j] - $alfm * $bR->[$i][$j]) * $bR->[$j][$en]; + } + + goto L630 if ($i == 0 || $isw == 2); + goto L630 if ($betm * $aR->[$i,$i-1] == 0); + + $zz = $w; + $s = $r; + goto L690; + + L630: + $m = $i; + goto L640 if ($isw == 2); + + $t = $w; + $t = $epsb if ($w == 0); + + $bR->[$i][$en] = -$r / $t; + next; + + L640: + $x = $betm * $aR->[$i][$i+1] - $alfm * $bR->[$i][$i+1]; + $y = $betm * $aR->[$i+1][$i]; + $q = $w * $zz - $x * $y; + $t = ($x * $s - $zz * $r) / $q; + $bR->[$i][$en] = $t; + goto L650 if (abs($x) <= abs($zz)); + $bR->[$i+1][$en] = (-$r - $w * $t) / $x; + goto L690; + + L650: + $bR->[$i+1][$en] = (-$s - $y * $t) / $zz; + + L690: + $isw = 3 - $isw; + + } # for ($ii inner loop + next; + + L710: + $m = $na; + $almr = $alfrR->[$m]; + $almi = $alfiR->[$m]; + $betm = $betaR->[$m]; + + $y = $betm * $aR->[$en][$na]; + $bR->[$na][$na] = -$almi * $bR->[$en][$en] / $y; + $bR->[$na][$en] = ($almr * $bR->[$en][$en] - $betm * $aR->[$en][$en]) / $y; + $bR->[$en][$na] = 0; + $bR->[$en][$en] = 1; + $enm2 = $na; + goto L795 if ($enm2 == 0); + + for ($ii=0; $ii<$enm2; $ii++) { + $i = $na - $ii - 1; + $w = $betm * $aR->[$i][$i] - $almr * $bR->[$i][$i]; + $w1 = -$almi * $bR->[$i][$i]; + $ra = $sa = 0; + + for ($j=$m; j<=$en; $j++) { + $x = $betm * $aR->[$i][$j] - $almr * $bR->[$i][$j]; + $x1 = -$almi * $bR->[$i][$j]; + $ra = $ra + $x * $bR->[$j][$na] - $x1 * $bR->[$j][$en]; + $sa = $sa + $x * $bR->[$j][$en] + $x1 * $bR->[$j][$na]; + } + + goto L770 if ($i == 0 || $isw == 2); + goto L770 if ($betm * $aR->[$i,$i-1] == 0); + + $zz = $w; $z1 = $w1; + $r = $ra; $s = $sa; + $isw = 2; + next; + + L770: + $m = $i; + goto L780 if ($isw == 2); + $tr = -$ra; $ti = -$sa; + + L773: + $dr = $w; $di = $w1; + + L775: + goto L777 if (abs($di) > abs($dr)); + $rr = $di / $dr; + $d = $dr + $di * $rr; + $t1 = ($tr + $ti * $rr) / $d; + $t2 = ($ti - $tr * $rr) / $d; + if ($isw == 1) { goto L787; } + elsif ($isw == 2) { goto L782; } + + L777: + $rr = $dr / $di; + $d = $dr * $rr + $di; + $t1 = ($tr * $rr + $ti) / $d; + $t2 = ($ti * $rr - $tr) / $d; + if ($isw == 1) { goto L787; } + elsif ($isw == 2) { goto L782; } + + L780: + $x = $betm * $aR->[$i][$i+1] - $almr * $bR->[$i][$i+1]; + $x1 = -$almi * $bR->[$i][$i+1]; + $y = $betm * $aR->[$i+1][$i]; + $tr = $y * $ra - $w * $r + $w1 * $s; + $ti = $y * $sa - $w * $s - $w1 * $r; + $dr = $w * $zz - $w1 * $z1 - $x * $y; + $di = $w * $z1 + $w1 * $zz - $x1 * $y; + $dr = $epsb if ($dr == 0 && $di == 0); + goto L775; + + L782: + $bR->[$i+1][$na] = $t1; + $bR->[$i+1][$en] = $t2; + $isw = 1; + goto L785 if (abs($y) > abs($w) + abs($w1)); + $tr = -$ra - $x * $bR->[$i+1][$na] + $x1 * $bR->[$i+1][$en]; + $ti = -$sa - $x * $bR->[$i+1][$en] - $x1 * $bR->[$i+1][$na]; + goto L773; + + L785: + $t1 = (-$r - $zz * $bR->[$i+1][$na] + $z1 * $bR->[$i+1][$en]) / $y; + $t2 = (-$s - $zz * $bR->[$i+1][$en] - $z1 * $bR->[$i+1][$na]) / $y; + + L787: + $bR->[$i][$na] = $t1; + $bR->[$i][$en] = $t2; + + } # for ($ii inner loop + + L795: + $isw = 3 - $isw; + + } # for ($nn outer loop + + for ($jj=0; $jj<$N; $jj++) { + $j = $N - $jj - 1; + for ($i=0; $i<$N; $i++) { + $zz = 0; + for ($k=0; $k<=$j; $k++) { + $zz += $zR->[$i][$k] * $bR->[$k][$j]; + } + $zR->[$i][$j] = $zz; + } + } + + for ($j=0; $j<$N; $j++) { + $d = 0; + goto L920 if ($isw == 2); + goto L945 if ($alfiR->[$j] != 0); + for ($i=0; $i<$N; $i++) { + $d = (abs($zR->[$i][$j])) if ((abs($zR->[$i][$j])) > $d); + } + for ($i=0; $i<$N; $i++) { + $zR->[$i][$j] /= $d; + } + next; + + L920: + for ($i=0; $i<$N; $i++) { + $r = abs($zR->[$i][$j-1]) + abs($zR->[$i][$j]); + if ($r != 0) { + my($u1) = $zR->[$i][$j-1] / $r; + my($u2) = $zR->[$i][$j] / $r; + $r *= sqrt($u1**2 + $u2**2); + } + $d = $r if ($r > $d); + } + + for ($i=0; $i<$N; $i++) { + $zR->[$i][$j-1] /= $d; + $zR->[$i][$j] /= $d; + } + + L945: + $isw = 3 - $isw; + + } # for ($j outer loop +} + +1; +