new file mode 100644
--- /dev/null
+++ b/libQZ.neardist
@@ -0,0 +1,1355 @@
+#======================================================================
+# 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;
+ a2 = a21;
+ if (l != ld)
+ a[l, lm1] = -a[l, lm1];
+ goto L160;
+
+ // Determine double shift zeroth column of a
+ L150:
+ a12 = a[l, l1] / b22;
+ a22 = a[l1, l1] / b22;
+ b12 = b[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 = a[l1 + 1, l1] / b22;
+ goto L160;
+
+ // Ad hoc shift
+ L155:
+ a1 = 0.0;
+ a2 = 1.0;
+ a3 = 1.1605;
+
+ L160:
+ ++its;
+ --itn;
+
+ if (!matz) lor1 = ld;
+
+ // Main loop
+ for (k = l; k <= na; ++k)
+ {
+ notlas = k != na && ish == 2;
+ k1 = k + 1;
+ k2 = k + 2;
+
+ km1 = Math.Max(k, l + 1) - 1; // Computing MAX
+ ll = Math.Min(en, k1 + ish); // Computing MIN
+
+ if (notlas) goto L190;
+
+ // Zero a(k+1,k-1)
+ if (k == l) goto L170;
+ a1 = a[k, km1];
+ a2 = a[k1, km1];
+
+ L170:
+ s = Math.Abs(a1) + Math.Abs(a2);
+ if (s == 0.0) goto L70;
+ u1 = a1 / s;
+ u2 = a2 / s;
+ r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);
+ v1 = -(u1 + r) / r;
+ v2 = -u2 / r;
+ u2 = v2 / v1;
+
+ for (j = km1; j < enorn; ++j)
+ {
+ t = a[k, j] + u2 * a[k1, j];
+ a[k, j] += t * v1;
+ a[k1, j] += t * v2;
+
+ t = b[k, j] + u2 * b[k1, j];
+ b[k, j] += t * v1;
+ b[k1, j] += t * v2;
+ }
+
+ if (k != l)
+ a[k1, km1] = 0.0;
+ goto L240;
+
+ // Zero a(k+1,k-1) and a(k+2,k-1)
+ L190:
+ if (k == l) goto L200;
+ a1 = a[k, km1];
+ a2 = a[k1, km1];
+ a3 = a[k2, km1];
+
+ L200:
+ s = Math.Abs(a1) + Math.Abs(a2) + Math.Abs(a3);
+ if (s == 0.0) goto L260;
+ u1 = a1 / s;
+ u2 = a2 / s;
+ u3 = a3 / s;
+ r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2 + u3 * u3), u1);
+ v1 = -(u1 + r) / r;
+ v2 = -u2 / r;
+ v3 = -u3 / r;
+ u2 = v2 / v1;
+ u3 = v3 / v1;
+
+ for (j = km1; j < enorn; ++j)
+ {
+ t = a[k, j] + u2 * a[k1, j] + u3 * a[k2, j];
+ a[k, j] += t * v1;
+ a[k1, j] += t * v2;
+ a[k2, j] += t * v3;
+
+ t = b[k, j] + u2 * b[k1, j] + u3 * b[k2, j];
+ b[k, j] += t * v1;
+ b[k1, j] += t * v2;
+ b[k2, j] += t * v3;
+ }
+
+ if (k == l) goto L220;
+ a[k1, km1] = 0.0;
+ a[k2, km1] = 0.0;
+
+ // Zero b(k+2,k+1) and b(k+2,k)
+ L220:
+ s = (Math.Abs(b[k2, k2])) + (Math.Abs(b[k2, k1])) + (Math.Abs(b[k2, k]));
+ if (s == 0.0) goto L240;
+ u1 = b[k2, k2] / s;
+ u2 = b[k2, k1] / s;
+ u3 = b[k2, k] / s;
+ r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2 + u3 * u3), 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 = a[i, k2] + u2 * a[i, k1] + u3 * a[i, k];
+ a[i, k2] += t * v1;
+ a[i, k1] += t * v2;
+ a[i, k] += t * v3;
+
+ t = b[i, k2] + u2 * b[i, k1] + u3 * b[i, k];
+ b[i, k2] += t * v1;
+ b[i, k1] += t * v2;
+ b[i, k] += t * v3;
+ }
+
+ b[k2, k] = 0.0;
+ b[k2, k1] = 0.0;
+
+ if (matz)
+ {
+ for (i = 0; i < n; ++i)
+ {
+ t = z[i, k2] + u2 * z[i, k1] + u3 * z[i, k];
+ z[i, k2] += t * v1;
+ z[i, k1] += t * v2;
+ z[i, k] += t * v3;
+ }
+ }
+
+ // Zero b(k+1,k)
+ L240:
+ s = (Math.Abs(b[k1, k1])) + (Math.Abs(b[k1, k]));
+ if (s == 0.0) goto L260;
+ u1 = b[k1, k1] / s;
+ u2 = b[k1, k] / s;
+ r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);
+ v1 = -(u1 + r) / r;
+ v2 = -u2 / r;
+ u2 = v2 / v1;
+
+ for (i = lor1; i < ll + 1; ++i)
+ {
+ t = a[i, k1] + u2 * a[i, k];
+ a[i, k1] += t * v1;
+ a[i, k] += t * v2;
+
+ t = b[i, k1] + u2 * b[i, k];
+ b[i, k1] += t * v1;
+ b[i, k] += t * v2;
+ }
+
+ b[k1, k] = 0.0;
+
+ if (matz)
+ {
+ for (i = 0; i < n; ++i)
+ {
+ t = z[i, k1] + u2 * z[i, k];
+ z[i, k1] += t * v1;
+ z[i, k] += t * v2;
+ }
+ }
+
+ L260:
+ ;
+ }
+
+ goto L70; // End QZ step
+
+ // Set error -- all eigenvalues have not converged after 30*n iterations
+ L1000:
+ ierr = en + 1;
+
+ // Save epsb for use by qzval and qzvec
+ L1001:
+ if (n > 1)
+ b[n - 1, 0] = epsb;
+ return 0;
+ }
+
+ /// <summary>
+ /// Adaptation of the original Fortran QZVAL routine from EISPACK.
+ /// </summary>
+ /// <remarks>
+ /// This subroutine is the third 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 matrices, one of them
+ /// in quasi-triangular form and the other in upper triangular form.
+ /// it reduces the quasi-triangular matrix further, so that any
+ /// remaining 2-by-2 blocks correspond to pairs of complex
+ /// eigenvalues, and returns quantities whose ratios give the
+ /// generalized eigenvalues. it is usually preceded by qzhes
+ /// and qzit and may be followed by qzvec.
+ ///
+ /// For the full documentation, please check the original function.
+ /// </remarks>
+ private static int qzval(int n, double[,] a, double[,] b, double[] alfr, double[] alfi, double[] beta, bool matz, double[,] z)
+ {
+ int i, j;
+ int na, en, nn;
+ double c, d, e = 0;
+ double r, s, t;
+ double a1, a2, u1, u2, v1, v2;
+ double a11, a12, a21, a22;
+ double b11, b12, b22;
+ double di, ei;
+ double an = 0, bn;
+ double cq, dr;
+ double cz, ti, tr;
+ double a1i, a2i, a11i, a12i, a22i, a11r, a12r, a22r;
+ double sqi, ssi, sqr, szi, ssr, szr;
+
+ double epsb = b[n - 1, 0];
+ int isw = 1;
+
+
+ // Find eigenvalues of quasi-triangular matrices.
+ for (nn = 0; nn < n; ++nn)
+ {
+ en = n - nn - 1;
+ na = en - 1;
+
+ if (isw == 2) goto L505;
+ if (en == 0) goto L410;
+ if (a[en, na] != 0.0) goto L420;
+
+ // 1-by-1 block, one real root
+ L410:
+ alfr[en] = a[en, en];
+ if (b[en, en] < 0.0)
+ {
+ alfr[en] = -alfr[en];
+ }
+ beta[en] = (Math.Abs(b[en, en]));
+ alfi[en] = 0.0;
+ goto L510;
+
+ // 2-by-2 block
+ L420:
+ if (Math.Abs(b[na, na]) <= epsb) goto L455;
+ if (Math.Abs(b[en, en]) > epsb) goto L430;
+ a1 = a[en, en];
+ a2 = a[en, na];
+ bn = 0.0;
+ goto L435;
+
+ L430:
+ an = Math.Abs(a[na, na]) + Math.Abs(a[na, en]) + Math.Abs(a[en, na]) + Math.Abs(a[en, en]);
+ bn = Math.Abs(b[na, na]) + Math.Abs(b[na, en]) + Math.Abs(b[en, en]);
+ a11 = a[na, na] / an;
+ a12 = a[na, en] / an;
+ a21 = a[en, na] / an;
+ a22 = a[en, en] / an;
+ b11 = b[na, na] / bn;
+ b12 = b[na, en] / bn;
+ b22 = b[en, en] / bn;
+ e = a11 / b11;
+ ei = a22 / b22;
+ s = a21 / (b11 * b22);
+ t = (a22 - e * b22) / b22;
+
+ if (Math.Abs(e) <= Math.Abs(ei))
+ goto L431;
+
+ e = ei;
+ t = (a11 - e * b11) / b11;
+
+ L431:
+ c = (t - s * b12) * .5;
+ d = c * c + s * (a12 - e * b12);
+ if (d < 0.0) goto L480;
+
+ // Two real roots. Zero both a(en,na) and b(en,na)
+ e += c + Special.Sign(Math.Sqrt(d), c);
+ a11 -= e * b11;
+ a12 -= e * b12;
+ a22 -= e * b22;
+
+ if (Math.Abs(a11) + Math.Abs(a12) < Math.Abs(a21) + Math.Abs(a22))
+ goto L432;
+
+ a1 = a12;
+ a2 = a11;
+ goto L435;
+
+ L432:
+ a1 = a22;
+ a2 = a21;
+
+ // Choose and apply real z
+ L435:
+ s = Math.Abs(a1) + Math.Abs(a2);
+ u1 = a1 / s;
+ u2 = a2 / 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 <= en; ++i)
+ {
+ t = a[i, en] + u2 * a[i, na];
+ a[i, en] += t * v1;
+ a[i, na] += t * v2;
+
+ t = b[i, en] + u2 * b[i, na];
+ b[i, en] += t * v1;
+ b[i, na] += t * v2;
+ }
+
+ if (matz)
+ {
+ for (i = 0; i < n; ++i)
+ {
+ t = z[i, en] + u2 * z[i, na];
+ z[i, en] += t * v1;
+ z[i, na] += t * v2;
+ }
+ }
+
+ if (bn == 0.0) goto L475;
+ if (an < System.Math.Abs(e) * bn) goto L455;
+ a1 = b[na, na];
+ a2 = b[en, na];
+ goto L460;
+
+ L455:
+ a1 = a[na, na];
+ a2 = a[en, na];
+
+ // Choose and apply real q
+ L460:
+ s = System.Math.Abs(a1) + System.Math.Abs(a2);
+ if (s == 0.0) goto L475;
+ u1 = a1 / s;
+ u2 = a2 / s;
+ r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);
+ v1 = -(u1 + r) / r;
+ v2 = -u2 / r;
+ u2 = v2 / v1;
+
+ for (j = na; j < n; ++j)
+ {
+ t = a[na, j] + u2 * a[en, j];
+ a[na, j] += t * v1;
+ a[en, j] += t * v2;
+
+ t = b[na, j] + u2 * b[en, j];
+ b[na, j] += t * v1;
+ b[en, j] += t * v2;
+ }
+
+ L475:
+ a[en, na] = 0.0;
+ b[en, na] = 0.0;
+ alfr[na] = a[na, na];
+ alfr[en] = a[en, en];
+
+ if (b[na, na] < 0.0)
+ alfr[na] = -alfr[na];
+
+ if (b[en, en] < 0.0)
+ alfr[en] = -alfr[en];
+
+ beta[na] = (System.Math.Abs(b[na, na]));
+ beta[en] = (System.Math.Abs(b[en, en]));
+ alfi[en] = 0.0;
+ alfi[na] = 0.0;
+ goto L505;
+
+ // Two complex roots
+ L480:
+ e += c;
+ ei = System.Math.Sqrt(-d);
+ a11r = a11 - e * b11;
+ a11i = ei * b11;
+ a12r = a12 - e * b12;
+ a12i = ei * b12;
+ a22r = a22 - e * b22;
+ a22i = ei * b22;
+
+ if (System.Math.Abs(a11r) + System.Math.Abs(a11i) +
+ System.Math.Abs(a12r) + System.Math.Abs(a12i) <
+ System.Math.Abs(a21) + System.Math.Abs(a22r)
+ + System.Math.Abs(a22i))
+ goto L482;
+
+ a1 = a12r;
+ a1i = a12i;
+ a2 = -a11r;
+ a2i = -a11i;
+ goto L485;
+
+ L482:
+ a1 = a22r;
+ a1i = a22i;
+ a2 = -a21;
+ a2i = 0.0;
+
+ // Choose complex z
+ L485:
+ cz = System.Math.Sqrt(a1 * a1 + a1i * a1i);
+ if (cz == 0.0) goto L487;
+ szr = (a1 * a2 + a1i * a2i) / cz;
+ szi = (a1 * a2i - a1i * a2) / cz;
+ r = System.Math.Sqrt(cz * cz + szr * szr + szi * szi);
+ cz /= r;
+ szr /= r;
+ szi /= r;
+ goto L490;
+
+ L487:
+ szr = 1.0;
+ szi = 0.0;
+
+ L490:
+ if (an < (System.Math.Abs(e) + ei) * bn) goto L492;
+ 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;
+
+ // Choose complex q
+ L495:
+ cq = System.Math.Sqrt(a1 * a1 + a1i * a1i);
+ if (cq == 0.0) goto L497;
+ sqr = (a1 * a2 + a1i * a2i) / cq;
+ sqi = (a1 * a2i - a1i * a2) / cq;
+ r = System.Math.Sqrt(cq * cq + sqr * sqr + sqi * sqi);
+ cq /= r;
+ sqr /= r;
+ sqi /= r;
+ goto L500;
+
+ L497:
+ sqr = 1.0;
+ sqi = 0.0;
+
+ // Compute diagonal elements that would result if transformations were applied
+ 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;
+
+ if (t < 0.0)
+ j = en;
+
+ r = Math.Sqrt(dr * dr + di * di);
+ beta[j] = bn * r;
+ alfr[j] = an * (tr * dr + ti * di) / r;
+ alfi[j] = an * t / r;
+ if (i == 0) goto L502;
+
+ L505:
+ isw = 3 - isw;
+
+ L510:
+ ;
+ }
+
+ b[n - 1, 0] = epsb;
+
+ return 0;
+ }
+
+ /// <summary>
+ /// Adaptation of the original Fortran QZVEC routine from EISPACK.
+ /// </summary>
+ /// <remarks>
+ /// This subroutine is the optional fourth 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 matrices, one of them in
+ /// quasi-triangular form (in which each 2-by-2 block corresponds to
+ /// a pair of complex eigenvalues) and the other in upper triangular
+ /// form. It computes the eigenvectors of the triangular problem and
+ /// transforms the results back to the original coordinate system.
+ /// it is usually preceded by qzhes, qzit, and qzval.
+ ///
+ /// For the full documentation, please check the original function.
+ /// </remarks>
+ private static int qzvec(int n, double[,] a, double[,] b, double[] alfr, double[] alfi, double[] beta, double[,] z)
+ {
+ int i, j, k, m;
+ int na, ii, en, jj, nn, enm2;
+ double d, q;
+ double r = 0, s = 0, t, w, x = 0, y, t1, t2, w1, x1 = 0, z1 = 0, di;
+ double ra, dr, sa;
+ double ti, rr, tr, zz = 0;
+ double alfm, almi, betm, almr;
+
+ double epsb = b[n - 1, 0];
+ int isw = 1;
+
+
+ // for en=n step -1 until 1 do --
+ for (nn = 0; nn < n; ++nn)
+ {
+ en = n - nn - 1;
+ na = en - 1;
+ if (isw == 2) goto L795;
+ if (alfi[en] != 0.0) goto L710;
+
+ // Real vector
+ m = en;
+ b[en, en] = 1.0;
+ if (na == -1) goto L800;
+ alfm = alfr[m];
+ betm = beta[m];
+
+ // for i=en-1 step -1 until 1 do --
+ for (ii = 0; ii <= na; ++ii)
+ {
+ i = en - ii - 1;
+ w = betm * a[i, i] - alfm * b[i, i];
+ r = 0.0;
+
+ for (j = m; j <= en; ++j)
+ r += (betm * a[i, j] - alfm * b[i, j]) * b[j, en];
+
+ if (i == 0 || isw == 2)
+ goto L630;
+
+ if (betm * a[i, i - 1] == 0.0)
+ goto L630;
+
+ zz = w;
+ s = r;
+ goto L690;
+
+ L630:
+ m = i;
+ if (isw == 2) goto L640;
+
+ // Real 1-by-1 block
+ t = w;
+ if (w == 0.0)
+ t = epsb;
+ b[i, en] = -r / t;
+ goto L700;
+
+ // Real 2-by-2 block
+ L640:
+ x = betm * a[i, i + 1] - alfm * b[i, i + 1];
+ y = betm * a[i + 1, i];
+ q = w * zz - x * y;
+ t = (x * s - zz * r) / q;
+ b[i, en] = t;
+ if (Math.Abs(x) <= Math.Abs(zz)) goto L650;
+ b[i + 1, en] = (-r - w * t) / x;
+ goto L690;
+
+ L650:
+ b[i + 1, en] = (-s - y * t) / zz;
+
+ L690:
+ isw = 3 - isw;
+
+ L700:
+ ;
+ }
+ // End real vector
+ goto L800;
+
+ // Complex vector
+ L710:
+ m = na;
+ almr = alfr[m];
+ almi = alfi[m];
+ betm = beta[m];
+
+ // last vector component chosen imaginary so that eigenvector matrix is triangular
+ y = betm * a[en, na];
+ b[na, na] = -almi * b[en, en] / y;
+ b[na, en] = (almr * b[en, en] - betm * a[en, en]) / y;
+ b[en, na] = 0.0;
+ b[en, en] = 1.0;
+ enm2 = na;
+ if (enm2 == 0) goto L795;
+
+ // for i=en-2 step -1 until 1 do --
+ for (ii = 0; ii < enm2; ++ii)
+ {
+ i = na - ii - 1;
+ w = betm * a[i, i] - almr * b[i, i];
+ w1 = -almi * b[i, i];
+ ra = 0.0;
+ sa = 0.0;
+
+ for (j = m; j <= en; ++j)
+ {
+ x = betm * a[i, j] - almr * b[i, j];
+ x1 = -almi * b[i, j];
+ ra = ra + x * b[j, na] - x1 * b[j, en];
+ sa = sa + x * b[j, en] + x1 * b[j, na];
+ }
+
+ if (i == 0 || isw == 2) goto L770;
+ if (betm * a[i, i - 1] == 0.0) goto L770;
+
+ zz = w;
+ z1 = w1;
+ r = ra;
+ s = sa;
+ isw = 2;
+ goto L790;
+
+ L770:
+ m = i;
+ if (isw == 2) goto L780;
+
+ // Complex 1-by-1 block
+ tr = -ra;
+ ti = -sa;
+
+ L773:
+ dr = w;
+ di = w1;
+
+ // Complex divide (t1,t2) = (tr,ti) / (dr,di)
+ L775:
+ if (Math.Abs(di) > Math.Abs(dr)) goto L777;
+ rr = di / dr;
+ d = dr + di * rr;
+ t1 = (tr + ti * rr) / d;
+ t2 = (ti - tr * rr) / d;
+
+ switch (isw)
+ {
+ case 1: goto L787;
+ case 2: goto L782;
+ }
+
+ L777:
+ rr = dr / di;
+ d = dr * rr + di;
+ t1 = (tr * rr + ti) / d;
+ t2 = (ti * rr - tr) / d;
+ switch (isw)
+ {
+ case 1: goto L787;
+ case 2: goto L782;
+ }
+
+ // Complex 2-by-2 block
+ L780:
+ x = betm * a[i, i + 1] - almr * b[i, i + 1];
+ x1 = -almi * b[i, i + 1];
+ y = betm * a[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;
+ if (dr == 0.0 && di == 0.0)
+ dr = epsb;
+ goto L775;
+
+ L782:
+ b[i + 1, na] = t1;
+ b[i + 1, en] = t2;
+ isw = 1;
+ if (Math.Abs(y) > Math.Abs(w) + Math.Abs(w1))
+ goto L785;
+ tr = -ra - x * b[(i + 1), na] + x1 * b[(i + 1), en];
+ ti = -sa - x * b[(i + 1), en] - x1 * b[(i + 1), na];
+ goto L773;
+
+ L785:
+ t1 = (-r - zz * b[(i + 1), na] + z1 * b[(i + 1), en]) / y;
+ t2 = (-s - zz * b[(i + 1), en] - z1 * b[(i + 1), na]) / y;
+
+ L787:
+ b[i, na] = t1;
+ b[i, en] = t2;
+
+ L790:
+ ;
+ }
+
+ // End complex vector
+ L795:
+ isw = 3 - isw;
+
+ L800:
+ ;
+ }
+
+ // End back substitution. Transform to original coordinate system.
+ for (jj = 0; jj < n; ++jj)
+ {
+ j = n - jj - 1;
+
+ for (i = 0; i < n; ++i)
+ {
+ zz = 0.0;
+ for (k = 0; k <= j; ++k)
+ zz += z[i, k] * b[k, j];
+ z[i, j] = zz;
+ }
+ }
+
+ // Normalize so that modulus of largest component of each vector is 1.
+ // (isw is 1 initially from before)
+ for (j = 0; j < n; ++j)
+ {
+ d = 0.0;
+ if (isw == 2) goto L920;
+ if (alfi[j] != 0.0) goto L945;
+
+ for (i = 0; i < n; ++i)
+ {
+ if ((Math.Abs(z[i, j])) > d)
+ d = (Math.Abs(z[i, j]));
+ }
+
+ for (i = 0; i < n; ++i)
+ z[i, j] /= d;
+
+ goto L950;
+
+ L920:
+ for (i = 0; i < n; ++i)
+ {
+ r = System.Math.Abs(z[i, j - 1]) + System.Math.Abs(z[i, j]);
+ if (r != 0.0)
+ {
+ // Computing 2nd power
+ double u1 = z[i, j - 1] / r;
+ double u2 = z[i, j] / r;
+ r *= Math.Sqrt(u1 * u1 + u2 * u2);
+ }
+ if (r > d)
+ d = r;
+ }
+
+ for (i = 0; i < n; ++i)
+ {
+ z[i, j - 1] /= d;
+ z[i, j] /= d;
+ }
+
+ L945:
+ isw = 3 - isw;
+
+ L950:
+ ;
+ }
+
+ return 0;
+ }
+
+ #endregion
+
+
+
+ #region ICloneable Members
+
+ private GeneralizedEigenvalueDecomposition()
+ {
+ }
+
+ /// <summary>
+ /// Creates a new object that is a copy of the current instance.
+ /// </summary>
+ /// <returns>
+ /// A new object that is a copy of this instance.
+ /// </returns>
+ public object Clone()
+ {
+ var clone = new GeneralizedEigenvalueDecomposition();
+ clone.ai = (double[])ai.Clone();
+ clone.ar = (double[])ar.Clone();
+ clone.beta = (double[])beta.Clone();
+ clone.n = n;
+ clone.Z = (double[,])Z.Clone();
+ return clone;
+ }
+
+ #endregion
+
+}
+