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

#                    L I B Q Z . P L 

#                    doc: Thu Mar 12 15:23:15 2015

#                    dlm: Thu Mar 12 20:52:43 2015

#                    (c) 2015 A.M. Thurnherr

#                    uE-Info: 36 0 NIL 0 0 72 2 2 4 NIL ofnI

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

# adaptation of EISPACK routines

# www.netlib.org/eispack

sub eig($$)

{

	my($aR,$bR) = @_;							# args passed as refs

	my($N) = scalar(@{aR});

	croak("eig(A,B): A & B must be matching square matrices\n")

		unless (@{bR} == $N) && (@{$aR->[0]} == $N) && (@{$bR->[0]} == $N);

		ar = new double[n];

		ai = new double[n];

		beta = new double[n];

		var A = (double[,])a.Clone();

		var B = (double[,])b.Clone();

	my($matZ) = 1;

	my($iErr) = 0;

	my(@Z);

	QZhes($aR,$bR,\@Z);								# reduce A/B to upper Hessenberg/triangular forms

	QZit($aR,$bR,\@Z,\$iErr);						# reduce Hess A to quasi-triangular form

	QZval($aR,$bR,\@Z,\@alphaR,\@alphaI,\@beta);	# reduce A further

	QZvec($aR,$bR,\@Z,\@alphaR,\@alphaI,\@beta);	# compute eigenvectors & eigenvalues

}

	/// <summary>Returns the real parts of the alpha values.</summary>

	public double[] RealAlphas

	{

		get { return ar; }

	}

	/// <summary>Returns the imaginary parts of the alpha values.</summary>

	public double[] ImaginaryAlphas

	{

		get { return ai; }

	}

	/// <summary>Returns the beta values.</summary>

	public double[] Betas

	{

		get { return beta; }

	}

	/// <summary>

	///   Returns true if matrix B is singular.

	/// </summary>

	/// <remarks>

	///   This method checks if any of the generated betas is zero. It

	///   does not says that the problem is singular, but only that one

	///   of the matrices of the pencil (A,B) is singular.

	/// </remarks>

	public bool IsSingular

	{

		get

		{

			for (int i = 0; i < n; i++)

				if (beta[i] == 0)

					return true;

			return false;

		}

	}

	/// <summary>

	///   Returns true if the eigenvalue problem is degenerate (ill-posed).

	/// </summary>

	public bool IsDegenerate

	{

		get

		{

			for (int i = 0; i < n; i++)

				if (beta[i] == 0 && ar[i] == 0)

					return true;

			return false;

		}

	}

	/// <summary>Returns the real parts of the eigenvalues.</summary>

	/// <remarks>

	///   The eigenvalues are computed using the ratio alpha[i]/beta[i],

	///   which can lead to valid, but infinite eigenvalues.

	/// </remarks>

	public double[] RealEigenvalues

	{

		get

		{

			// ((alfr+i*alfi)/beta)

			double[] eval = new double[n];

			for (int i = 0; i < n; i++)

				eval[i] = ar[i] / beta[i];

			return eval;

		}

	}

	/// <summary>Returns the imaginary parts of the eigenvalues.</summary>  

	/// <remarks>

	///   The eigenvalues are computed using the ratio alpha[i]/beta[i],

	///   which can lead to valid, but infinite eigenvalues.

	/// </remarks>

	public double[] ImaginaryEigenvalues

	{

		get

		{

			// ((alfr+i*alfi)/beta)

			double[] eval = new double[n];

			for (int i = 0; i < n; i++)

				eval[i] = ai[i] / beta[i];

			return eval;

		}

	}

	/// <summary>Returns the eigenvector matrix.</summary>

	public double[,] Eigenvectors

	{

		get

		{

			return Z;

		}

	}

	/// <summary>Returns the block diagonal eigenvalue matrix.</summary>

	public double[,] DiagonalMatrix

	{

		get

		{

			double[,] x = new double[n, n];

			for (int i = 0; i < n; i++)

			{

				for (int j = 0; j < n; j++)

					x[i, j] = 0.0;

				x[i, i] = ar[i] / beta[i];

				if (ai[i] > 0)

					x[i, i + 1] = ai[i] / beta[i];

				else if (ai[i] < 0)

					x[i, i - 1] = ai[i] / beta[i];

			}

			return x;

		}

	}

	#region EISPACK Routines

	/// <summary>

	///   Adaptation of the original Fortran QZHES routine from EISPACK.

	/// </summary>

	/// <remarks>

	///   This subroutine is the first step of the qz algorithm

	///   for solving generalized matrix eigenvalue problems,

	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart.

	///

	///   This subroutine accepts a pair of real general matrices and

	///   reduces one of them to upper hessenberg form and the other

	///   to upper triangular form using orthogonal transformations.

	///   it is usually followed by  qzit,	qzval  and, possibly,  qzvec.

	///  

	///   For the full documentation, please check the original function.

	/// </remarks>

	private static int qzhes(int n, double[,] a, double[,] b, bool matz, double[,] z)

	{

		int i, j, k, l;

		double r, s, t;

		int l1;

		double u1, u2, v1, v2;

		int lb, nk1;

		double rho;

		if (matz)

		{

			// If we are interested in computing the

			//	eigenvectors, set Z to identity(n,n)

			for (j = 0; j < n; ++j)

			{

				for (i = 0; i < n; ++i)

					z[i, j] = 0.0;

				z[j, j] = 1.0;

			}

		}

		// Reduce b to upper triangular form

		if (n <= 1) return 0;

		for (l = 0; l < n - 1; ++l)

		{

			l1 = l + 1;

			s = 0.0;

			for (i = l1; i < n; ++i)

				s += (System.Math.Abs(b[i, l]));

			if (s == 0.0) continue;

			s += (System.Math.Abs(b[l, l]));

			r = 0.0;

			for (i = l; i < n; ++i)

			{

				// Computing 2nd power

				b[i, l] /= s;

				r += b[i, l] * b[i, l];

			}

			r = Special.Sign(System.Math.Sqrt(r), b[l, l]);

			b[l, l] += r;

			rho = r * b[l, l];

			for (j = l1; j < n; ++j)

			{

				t = 0.0;

				for (i = l; i < n; ++i)

					t += b[i, l] * b[i, j];

				t = -t / rho;

				for (i = l; i < n; ++i)

					b[i, j] += t * b[i, l];

			}

			for (j = 0; j < n; ++j)

			{

				t = 0.0;

				for (i = l; i < n; ++i)

					t += b[i, l] * a[i, j];

				t = -t / rho;

				for (i = l; i < n; ++i)

					a[i, j] += t * b[i, l];

			}

			b[l, l] = -s * r;

			for (i = l1; i < n; ++i)

				b[i, l] = 0.0;

		}

		// Reduce a to upper hessenberg form, while keeping b triangular

		if (n == 2) return 0;

		for (k = 0; k < n - 2; ++k)

		{

			nk1 = n - 2 - k;

			// for l=n-1 step -1 until k+1 do

			for (lb = 0; lb < nk1; ++lb)

			{

				l = n - lb - 2;

				l1 = l + 1;

				// Zero a(l+1,k)

				s = (System.Math.Abs(a[l, k])) + (System.Math.Abs(a[l1, k]));

				if (s == 0.0) continue;

				u1 = a[l, k] / s;

				u2 = a[l1, k] / s;

				r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

				v1 = -(u1 + r) / r;

				v2 = -u2 / r;

				u2 = v2 / v1;

				for (j = k; j < n; ++j)

				{

					t = a[l, j] + u2 * a[l1, j];

					a[l, j] += t * v1;

					a[l1, j] += t * v2;

				}

				a[l1, k] = 0.0;

				for (j = l; j < n; ++j)

				{

					t = b[l, j] + u2 * b[l1, j];

					b[l, j] += t * v1;

					b[l1, j] += t * v2;

				}

				// Zero b(l+1,l)

				s = (System.Math.Abs(b[l1, l1])) + (System.Math.Abs(b[l1, l]));

				if (s == 0.0) continue;

				u1 = b[l1, l1] / s;

				u2 = b[l1, l] / s;

				r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

				v1 = -(u1 + r) / r;

				v2 = -u2 / r;

				u2 = v2 / v1;

				for (i = 0; i <= l1; ++i)

				{

					t = b[i, l1] + u2 * b[i, l];

					b[i, l1] += t * v1;

					b[i, l] += t * v2;

				}

				b[l1, l] = 0.0;

				for (i = 0; i < n; ++i)

				{

					t = a[i, l1] + u2 * a[i, l];

					a[i, l1] += t * v1;

					a[i, l] += t * v2;

				}

				if (matz)

				{

					for (i = 0; i < n; ++i)

					{

						t = z[i, l1] + u2 * z[i, l];

						z[i, l1] += t * v1;

						z[i, l] += t * v2;

					}

				}

			}

		}

		return 0;

	}

	/// <summary>

	///   Adaptation of the original Fortran QZIT routine from EISPACK.

	/// </summary>

	/// <remarks>

	///   This subroutine is the second step of the qz algorithm

	///   for solving generalized matrix eigenvalue problems,

	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart,

	///   as modified in technical note nasa tn d-7305(1973) by ward.

	///  

	///   This subroutine accepts a pair of real matrices, one of them

	///   in upper hessenberg form and the other in upper triangular form.

	///   it reduces the hessenberg matrix to quasi-triangular form using

	///   orthogonal transformations while maintaining the triangular form

	///   of the other matrix.	it is usually preceded by  qzhes  and

	///   followed by  qzval  and, possibly,  qzvec.

	///  

	///   For the full documentation, please check the original function.

	/// </remarks>

	private static int qzit(int n, double[,] a, double[,] b, double eps1, bool matz, double[,] z, ref int ierr)

	{

		int i, j, k, l = 0;

		double r, s, t, a1, a2, a3 = 0;

		int k1, k2, l1, ll;

		double u1, u2, u3;

		double v1, v2, v3;

		double a11, a12, a21, a22, a33, a34, a43, a44;

		double b11, b12, b22, b33, b34, b44;

		int na, en, ld;

		double ep;

		double sh = 0;

		int km1, lm1 = 0;

		double ani, bni;

		int ish, itn, its, enm2, lor1;

		double epsa, epsb, anorm = 0, bnorm = 0;

		int enorn;

		bool notlas;

		ierr = 0;

		#region Compute epsa and epsb

		for (i = 0; i < n; ++i)

		{

			ani = 0.0;

			bni = 0.0;

			if (i != 0)

				ani = (Math.Abs(a[i, (i - 1)]));

			for (j = i; j < n; ++j)

			{

				ani += Math.Abs(a[i, j]);

				bni += Math.Abs(b[i, j]);

			}

			if (ani > anorm) anorm = ani;

			if (bni > bnorm) bnorm = bni;

		}

		if (anorm == 0.0) anorm = 1.0;

		if (bnorm == 0.0) bnorm = 1.0;

		ep = eps1;

		if (ep == 0.0)

		{

			// Use roundoff level if eps1 is zero

			ep = Special.Epslon(1.0);

		}

		epsa = ep * anorm;

		epsb = ep * bnorm;

		#endregion

		// Reduce a to quasi-triangular form, while keeping b triangular

		lor1 = 0;

		enorn = n;

		en = n - 1;

		itn = n * 30;

	// Begin QZ step

	L60:

		if (en <= 1) goto L1001;

		if (!matz) enorn = en + 1;

		its = 0;

		na = en - 1;

		enm2 = na;

	L70:

		ish = 2;

		// Check for convergence or reducibility.

		for (ll = 0; ll <= en; ++ll)

		{

			lm1 = en - ll - 1;

			l = lm1 + 1;

			if (l + 1 == 1)

				goto L95;

			if ((Math.Abs(a[l, lm1])) <= epsa)

				break;

		}

	L90:

		a[l, lm1] = 0.0;

		if (l < na) goto L95;

		// 1-by-1 or 2-by-2 block isolated

		en = lm1;

		goto L60;

	// Check for small top of b

	L95:

		ld = l;

	L100:

		l1 = l + 1;

		b11 = b[l, l];

		if (Math.Abs(b11) > epsb) goto L120;

		b[l, l] = 0.0;

		s = (Math.Abs(a[l, l]) + Math.Abs(a[l1, l]));

		u1 = a[l, l] / s;

		u2 = a[l1, l] / s;

		r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

		v1 = -(u1 + r) / r;

		v2 = -u2 / r;

		u2 = v2 / v1;

		for (j = l; j < enorn; ++j)

		{

			t = a[l, j] + u2 * a[l1, j];

			a[l, j] += t * v1;

			a[l1, j] += t * v2;

			t = b[l, j] + u2 * b[l1, j];

			b[l, j] += t * v1;

			b[l1, j] += t * v2;

		}

		if (l != 0)

			a[l, lm1] = -a[l, lm1];

		lm1 = l;

		l = l1;

		goto L90;

	L120:

		a11 = a[l, l] / b11;

		a21 = a[l1, l] / b11;

		if (ish == 1) goto L140;

		// Iteration strategy

		if (itn == 0) goto L1000;

		if (its == 10) goto L155;

		// Determine type of shift

		b22 = b[l1, l1];

		if (Math.Abs(b22) < epsb) b22 = epsb;

		b33 = b[na, na];

		if (Math.Abs(b33) < epsb) b33 = epsb;

		b44 = b[en, en];

		if (Math.Abs(b44) < epsb) b44 = epsb;

		a33 = a[na, na] / b33;

		a34 = a[na, en] / b44;

		a43 = a[en, na] / b33;

		a44 = a[en, en] / b44;

		b34 = b[na, en] / b44;

		t = (a43 * b34 - a33 - a44) * .5;

		r = t * t + a34 * a43 - a33 * a44;

		if (r < 0.0) goto L150;

		// Determine single shift zeroth column of a

		ish = 1;

		r = Math.Sqrt(r);

		sh = -t + r;

		s = -t - r;

		if (Math.Abs(s - a44) < Math.Abs(sh - a44))

			sh = s;

		// Look for two consecutive small sub-diagonal elements of a.

		for (ll = ld; ll + 1 <= enm2; ++ll)

		{

			l = enm2 + ld - ll - 1;

			if (l == ld)

				goto L140;

			lm1 = l - 1;

			l1 = l + 1;

			t = a[l + 1, l + 1];

			if (Math.Abs(b[l, l]) > epsb)

				t -= sh * b[l, l];

			if (Math.Abs(a[l, lm1]) <= (Math.Abs(t / a[l1, l])) * epsa)

				goto L100;

		}

	L140:

		a1 = a11 - sh;