#======================================================================

#                    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;