new file mode 100644
--- /dev/null
+++ b/HISTORY
@@ -0,0 +1,11 @@
+#======================================================================
+# H I S T O R Y
+# doc: Thu May 7 13:12:05 2015
+# dlm: Thu May 7 13:12:52 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 10 10 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+May 7, 2015:
+ - V6.0 [ants.pl] [.hg/hgrc] published for release of LADCPproc V1.2 (Slocum/Explorer processing
+
--- a/INDEX
+++ b/INDEX
@@ -1,9 +1,9 @@
======================================================================
I N D E X
doc: Wed Jun 18 09:46:58 1997
- dlm: Mon Nov 3 12:42:00 2014
+ dlm: Fri Mar 13 20:55:01 2015
(c) 1997 Andreas Thurnherr
- uE-Info: 186 41 NIL 0 0 72 2 2 4 NIL ofnI
+ uE-Info: 195 48 NIL 0 0 72 2 2 4 NIL ofnI
======================================================================
NOTES:
@@ -189,7 +189,10 @@
[libvec.pl] vector stuff, including distance on globe
[libWOCE.pl] WOCE conversions
-=Numerical Recipes Routines=
+=Numerical Routines=
+
+[libQZ.pl] generalized eigenvalue problem (eig(A,B))
+[libSVD.pl] singular value decomposition
[covsrt.pl] for [lmfit]
[fft.pl] FFT
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
+
+}
+
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;
+
new file mode 100644
--- /dev/null
+++ b/libSVD.pl
@@ -0,0 +1,353 @@
+#======================================================================
+# L I B S V D . P L
+# doc: Sat Jul 31 22:47:03 1999
+# dlm: Fri Mar 13 20:53:26 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 305 1 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Jul 31, 1999: - created
+# Jul 19, 2001: - done *something* (message only?)
+# Mar 11, 2015: - fixed syntax errors (code had never been used before)
+# Mar 13, 2015: - combined from [sbbksb.pl] [svdcmp.pl] [pythag.pl] [svdfit.pl]
+
+require "$ANTS/nrutil.pl";
+use strict;
+
+#----------------------------------------------------------------------
+# SVBKSB routine from Numerical Recipes adapted to ANTS
+#
+# solves Ax = b for x, given b
+#
+# Notes:
+# - A = U W V' as done in [svdcmp.pl]
+#----------------------------------------------------------------------
+
+sub svbksb($$$$$)
+{
+ my($uR,$wR,$vR,$bR,$xR) = @_;
+ my($jj,$j,$i); # int
+ my($s); # float
+ my(@tmp); # float[]
+
+ &vector(\@tmp,1,$#{$wR});
+ for ($j=1; $j<=$#{$wR}; $j++) {
+ $s = 0;
+ if ($wR->[$j]) {
+ for ($i=1; $i<=$#{$uR}; $i++) {
+ $s += $uR->[$i][$j] * $bR->[$i];
+ }
+ $s /= $wR->[$j];
+ }
+ $tmp[$j]=$s;
+ }
+ for ($j=1; $j<=$#{$wR}; $j++) {
+ $s = 0;
+ for ($jj=1; $jj<=$#{$wR}; $jj++) {
+ $s += $vR->[$j][$jj] * $tmp[$jj];
+ }
+ $xR->[$j] = $s;
+ }
+}
+
+#----------------------------------------------------------------------
+# PYTHAG routine
+#----------------------------------------------------------------------
+
+sub pythag($$)
+{
+ my($a,$b) = @_; # params
+ my($absa,$absb); # float
+
+ $absa = abs($a);
+ $absb = abs($b);
+ return $absa*sqrt(1.0+SQR($absb/$absa))
+ if ($absa > $absb);
+ return ($absb == 0 ? 0 : $absb*sqrt(1+$absa*$absa/$absb/$absb));
+}
+
+
+#----------------------------------------------------------------------
+# SVDCMP routine from Numerical Recipes adapted to ANTS
+#----------------------------------------------------------------------
+
+sub svdcmp($$$)
+{
+ my($aR,$wR,$vR) = @_; # params
+ my($flag,$i,$its,$j,$jj,$k,$l,$nm); # int
+ my($anorm,$c,$f,$g,$h,$s,$scale,$x,$y,$z); # float
+ my(@rv1); # float[]
+
+ vector(\@rv1,1,$#{$vR});
+ $g = $scale = $anorm = 0;
+ for ($i=1; $i<=$#{$vR}; $i++) {
+ $l = $i+1;
+ $rv1[$i] = $scale*$g;
+ $g = $s = $scale = 0;
+ if ($i <= $#{$aR}) {
+ for ($k=$i; $k<=$#{$aR}; $k++) {
+ $scale += abs($aR->[$k][$i]);
+ }
+ if ($scale) {
+ for ($k=$i; $k<=$#{$aR}; $k++) {
+ $aR->[$k][$i] /= $scale;
+ $s += $aR->[$k][$i]*$aR->[$k][$i];
+ }
+ $f = $aR->[$i][$i];
+ $g = -&SIGN(sqrt($s),$f);
+ $h = $f*$g-$s;
+ $aR->[$i][$i] = $f-$g;
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ for ($s=0,$k=$i; $k<=$#{$aR}; $k++) {
+ $s += $aR->[$k][$i]*$aR->[$k][$j];
+ }
+ $f = $s/$h;
+ for ($k=$i; $k<=$#{$aR}; $k++) {
+ $aR->[$k][$j] += $f*$aR->[$k][$i];
+ }
+ }
+ for ($k=$i; $k<=$#{$aR}; $k++) {
+ $aR->[$k][$i] *= $scale;
+ }
+ }
+ }
+ $wR->[$i] = $scale * $g;
+ $g = 0; $s = 0; $scale = 0;
+ if ($i <= $#{$aR} && $i != $#{$vR}) {
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $scale += abs($aR->[$i][$k]);
+ }
+ if ($scale) {
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $aR->[$i][$k] /= $scale;
+ $s += $aR->[$i][$k]*$aR->[$i][$k];
+ }
+ $f = $aR->[$i][$l];
+ $g = -&SIGN(sqrt($s),$f);
+ $h = $f*$g-$s;
+ $aR->[$i][$l] = $f-$g;
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $rv1[$k] = $aR->[$i][$k]/$h;
+ }
+ for ($j=$l; $j<=$#{$aR}; $j++) {
+ for ($s=0,$k=$l; $k<=$#{$vR}; $k++) {
+ $s += $aR->[$j][$k]*$aR->[$i][$k];
+ }
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $aR->[$j][$k] += $s*$rv1[$k];
+ }
+ }
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $aR->[$i][$k] *= $scale;
+ }
+ }
+ }
+ $anorm = MAX($anorm,(abs($wR->[$i])+abs($rv1[$i])));
+ }
+ for ($i=$#{$vR}; $i>=1; $i--) {
+ if ($i < $#{$vR}) {
+ if ($g) {
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ $vR->[$j][$i] = ($aR->[$i][$j]/$aR->[$i][$l])/$g;
+ }
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ for ($s=0,$k=$l; $k<=$#{$vR}; $k++) {
+ $s += $aR->[$i][$k]*$vR->[$k][$j];
+ }
+ for ($k=$l; $k<=$#{$vR}; $k++) {
+ $vR->[$k][$j] += $s*$vR->[$k][$i];
+ }
+ }
+ }
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ $vR->[$i][$j] = 0; $vR->[$j][$i] = 0;
+ }
+ }
+ $vR->[$i][$i] = 1;
+ $g = $rv1[$i];
+ $l = $i;
+ }
+ for ($i=MIN($#{$aR},$#{$vR}); $i>=1; $i--) {
+ $l = $i+1;
+ $g = $wR->[$i];
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ $aR->[$i][$j] = 0;
+ }
+ if ($g) {
+ $g = 1/$g;
+ for ($j=$l; $j<=$#{$vR}; $j++) {
+ for ($s=0,$k=$l; $k<=$#{$aR}; $k++) {
+ $s += $aR->[$k][$i]*$aR->[$k][$j];
+ }
+ $f = ($s/$aR->[$i][$i])*$g;
+ for ($k=$i; $k<=$#{$aR}; $k++) {
+ $aR->[$k][$j] += $f*$aR->[$k][$i];
+ }
+ }
+ for ($j=$i; $j<=$#{$aR}; $j++) {
+ $aR->[$j][$i] *= $g;
+ }
+ } else {
+ for ($j=$i; $j<=$#{$aR}; $j++) {
+ $aR->[$j][$i] = 0;
+ }
+ }
+ ++$aR->[$i][$i];
+ }
+ for ($k=$#{$vR}; $k>=1; $k--) {
+ for ($its=1; $its<=30; $its++) {
+ $flag = 1;
+ for ($l=$k; $l>=1; $l--) {
+ $nm = $l-1;
+ if ((abs($rv1[$l])+$anorm) == $anorm) {
+ $flag = 0;
+ last;
+ }
+ last if ($nm == 0) || ((abs($wR->[$nm])+$anorm) == $anorm); ## $nm == 0 test not in original code
+ }
+ if ($flag) {
+ $c = 0;
+ $s = 1;
+ for ($i=$l; $i<=$k; $i++) {
+ $f = $s*$rv1[$i];
+ $rv1[$i] = $c*$rv1[$i];
+ last if ((abs($f)+$anorm) == $anorm);
+ $g = $wR->[$i];
+ $h = &pythag($f,$g);
+ $wR->[$i] = $h;
+ $h = 1/$h;
+ $c = $g*$h;
+ $s = -$f*$h;
+ for ($j=1; $j<=$#{$aR}; $j++) {
+ $y = $aR->[$j][$nm];
+ $z = $aR->[$j][$i];
+ $aR->[$j][$nm] = $y*$c+$z*$s;
+ $aR->[$j][$i] = $z*$c-$y*$s;
+ }
+ }
+ }
+ $z = $wR->[$k];
+ if ($l == $k) {
+ if ($z < 0) {
+ $wR->[$k] = -$z;
+ for ($j=1; $j<=$#{$vR}; $j++) {
+ $vR->[$j][$k] = -$vR->[$j][$k];
+ }
+ }
+# print(STDERR "its($k) = $its\n");
+ last;
+ }
+ croak("no convergence in 30 svdcmp iterations\n") if ($its == 30);
+ $x = $wR->[$l];
+ $nm = $k-1;
+ $y = $wR->[$nm];
+ $g = $rv1[$nm];
+ $h = $rv1[$k];
+ $f = (($y-$z)*($y+$z)+($g-$h)*($g+$h))/(2.0*$h*$y);
+ $g = &pythag($f,1);
+ $f = (($x-$z)*($x+$z)+$h*(($y/($f+&SIGN($g,$f)))-$h))/$x;
+ $c = 1; $s = 1;
+ for ($j=$l; $j<=$nm; $j++) {
+ $i = $j+1;
+ $g = $rv1[$i];
+ $y = $wR->[$i];
+ $h = $s*$g;
+ $g = $c*$g;
+ $z = &pythag($f,$h);
+ $rv1[$j] = $z;
+ $c = $f/$z;
+ $s = $h/$z;
+ $f = $x*$c+$g*$s;
+ $g = $g*$c-$x*$s;
+ $h = $y*$s;
+ $y *= $c;
+ for ($jj=1; $jj<=$#{$vR}; $jj++) {
+ $x = $vR->[$jj][$j];
+ $z = $vR->[$jj][$i];
+ $vR->[$jj][$j] = $x*$c+$z*$s;
+ $vR->[$jj][$i] = $z*$c-$x*$s;
+ }
+ $z = &pythag($f,$h);
+ $wR->[$j] = $z;
+ if ($z) {
+ $z = 1/$z;
+ $c = $f*$z;
+ $s = $h*$z;
+ }
+ $f = $c*$g+$s*$y;
+ $x = $c*$y-$s*$g;
+ for ($jj=1; $jj<=$#{$aR}; $jj++) {
+ $y = $aR->[$jj][$j];
+ $z = $aR->[$jj][$i];
+ $aR->[$jj][$j] = $y*$c+$z*$s;
+ $aR->[$jj][$i] = $z*$c-$y*$s;
+ }
+ }
+ $rv1[$l] = 0;
+ $rv1[$k] = $f;
+ $wR->[$k] = $x;
+ }
+ }
+}
+
+#----------------------------------------------------------------------
+# SVDFIT routine from Numerical Recipes adapted to ANTS
+#
+# UNTESTED CODE!!!!!
+#
+# Notes:
+# - x,y,sig are field numbers for data in $ants_
+# - if sig is a negative number, -sig is used as constant input stddev
+# - @a, @u, @v, @w, &funcs passed as refs
+# - chi square is returned
+#----------------------------------------------------------------------
+#
+#{ # BEGIN static scope
+#
+#my($TOL) = 1.0e-5;
+#
+#sub svdfit($$$$$$$$)
+#{
+# die("untested code");
+# my($xfnr,$yfnr,$sig,$aR,$uR,$vR,$wR,$funcsR) = @_;
+# my($j,$i); # int
+# my($chisq,$wmax,$tmp,$thresh,$sum); # float
+# my(@b,@afunc); # float[]
+#
+# &vector(\@b,1,$#ants_);
+# &vector(\@afunc,1,$#{$aR});
+# for ($i=0; $i<=$#ants_; $i++) {
+# next if ($antsFlagged[$i]);
+# &$funcsR($ants_[$i][$xfnr],\@afunc);
+# $tmp = 1.0 / (($sig > 0) ? $ants_[$i][$sig] : -$sig);
+# for ($j=1; $j<=$#{$aR}; $j++) {
+# $uR->[$i][$j] = $afunc[$j]*$tmp;
+# }
+# $b[$i] = $ants_[$i][$yfnr]*$tmp;
+# }
+# &svdcmp($uR,$wR,$vR);
+# for ($j=1; $j<=$#{$aR}; $j++) {
+# $wmax = $wR->[$j] if ($wR->[$j] > $wmax);
+# }
+# $thresh = $TOL*$wmax;
+# for ($j=1; $j<=$#{$aR}; $j++) {
+# $wR->[$j] = 0 if ($wR->[$j] < $thresh);
+# }
+# &svbksb($uR,$wR,$vR,\@b,$aR);
+# for ($i=0; $i<=$#ants_; $i++) {
+# next if ($antsFlagged[$i]);
+# &$funcsR($ants_[$i][$xfnr],\@afunc);
+# for ($j=1; $j<=$#{$aR}; $j++) {
+# $sum += $aR->[$j]*$afunc[$j];
+# }
+# $tmp = ($ants_[$i][$yfnr] - $sum) /
+# (($sig > 0) ? $ants_[$i][$sig] : -$sig);
+# $chisq += $tmp * $tmp;
+# }
+# return $chisq;
+#}
+#
+#} # END static scope
+
+1;
new file mode 100644
--- /dev/null
+++ b/libloopfit.pl
@@ -0,0 +1,114 @@
+#======================================================================
+# L I B L O O P F I T . P L
+# doc: Thu Mar 12 08:05:26 2015
+# dlm: Sun Mar 15 19:46:49 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 47 22 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Mar 12, 2015: - created
+# Mar 15, 2015: - added [libQZ.pl]
+
+require "$ANTS/nrutil.pl";
+require "$ANTS/libSVD.pl"; # for circular fit
+require "$ANTS/libQZ.pl"; # for ellipse fit
+
+#----------------------------------------------------------------------
+# fitCircle(\@x,\@y) => ($r,$cx,$cy)
+# \@x,\@y references to list of data values
+# $r circle radius
+# $cx,$cy circle center
+#----------------------------------------------------------------------
+
+sub fitCircle($$)
+{
+ my($xR,$yR) = @_;
+ my(@A,@V,@W);
+ my($nSamp) = scalar(@{$xR});
+ my($nParam) = 3;
+
+ matrix(\@A,1,$nSamp,1,$nParam); # rows, then columns
+ matrix(\@V,1,$nParam,1,$nParam); # 3x3 matrix for circle fitting
+ vector(\@W,1,$nParam); # diag matrix of singular values output as vector
+ vector(\@b,1,$nSamp); # rhs
+ vector(\@coeff,1,$nParam); # fit coefficients
+
+ for ($r=1; $r<=$nSamp; $r++) { # circle equation c1*x + c2*y + c3 = -(x^2 + y^2)
+ $A[$r][1] = $xR->[$r-1]; # r = sqrt((c1^2 + c2^2)/4 - c3)
+ $A[$r][2] = $yR->[$r-1]; # cx = -c1/2
+ $A[$r][3] = 1; # cy = -c2/2
+ $b[$r] = -($xR->[$r-1]**2+$yR->[$r-1]**2);
+ }
+
+ svdcmp(\@A,\@W,\@V); # solve Ax = b, with x == coeff
+ svbksb(\@A,\@W,\@V,\@b,\@coeff);
+
+ return ($coeff[1],$coeff[2],$coeff[3]);
+ my($r) = sqrt(($coeff[1]**2+$coeff[2]**2)/4-$coeff[3]);
+ my($cx) = -0.5*$coeff[1];
+ my($cy) = -0.5*$coeff[2];
+ return ($r,$cx,$cy);
+}
+
+#----------------------------------------------------------------------
+# fitEllipse(\@x,\@y) => ?
+# \@x,\@y references to list of data values
+#
+# Direct Least Square Fitting of Ellipses
+# Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher, MAY 1999
+# IEEE transactions on pattern analysis and machine intelligence 21(5): 476-480
+#----------------------------------------------------------------------
+
+sub fitEllipse($$)
+{
+ my($xR,$yR) = @_;
+ my($nSamp) = scalar(@{$xR});
+
+ # design matrix D = [ x.*x x.*y y.*y x y ones(size(x)) ];
+ my(@D);
+ for (my($r)=0; $r<$nSamp; $r++) {
+ $D[$r][0] = $xR->[$r]**2;
+ $D[$r][1] = $xR->[$r] * $yR->[$r];
+ $D[$r][2] = $yR->[$r]**2;
+ $D[$r][3] = $xR->[$r];
+ $D[$r][4] = $yR->[$r];
+ $D[$r][5] = 1;
+ }
+
+ # scatter matrix S = D' * D;
+ my(@S);
+ for (my($i)=0; $i<6; $i++) {
+ for (my($j)=0; $j<6; $j++) {
+ for (my($k)=0; $k<$nSamp; $k++) {
+ $S[$i][$j] += $D[$k][$j] * $D[$k][$i];
+ }
+ }
+ }
+
+ # 6x6 constraint matrix C
+ my(@C);
+ for (my($r)=0; $r<6; $r++) {
+ for (my($c)=0; $c<6; $c++) {
+ $C[$r][$c] = 0;
+ }
+ }
+ $C[0][2] = -2; $C[1][1] = 1; $C[2][0] = -2;
+
+ # solve generalized eigensystem
+ my(@geRe,@geImg,@geVec);
+ eig(\@S,\@C,\@geRe,\@geImg,\@geVec);
+
+ # find only negative eigenvalue
+ my($i);
+ for ($i=0; $i<@geRe; $i++) {
+ last if (numberp($geRe[$i]) && ($geRe[$i] < 0));
+ }
+ croak("internal error") unless ($geRe[$i] < 0);
+
+ # get fitted parameters
+ return @{$geVec[$i]};
+
+}
+
+1;
--- a/pythag.pl
+++ b/pythag.pl
@@ -1,15 +1,16 @@
#======================================================================
-# . / P Y T H A G . P L
+# P Y T H A G . P L
# doc: Sun Aug 1 10:41:34 1999
-# dlm: Sun Aug 1 10:46:43 1999
+# dlm: Wed Mar 11 14:53:21 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 23 65 NIL 0 0 72 0 2 4 ofnI
+# uE-Info: 27 2 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# PYTHAG routine from Numerical Recipes adapted to ANTS
# HISTORY:
# Aug 01, 1999: - manually converted from c-source
+# Mar 11, 2015: - BUG: syntax errors (code had never been used)
sub pythag($$)
{
@@ -20,6 +21,7 @@
$absb = abs($b);
return $absa*sqrt(1.0+SQR($absb/$absa))
if ($absa > $absb);
- return ($absb == 0 ? 0 : $absb*sqrt(1+$absa*$absa/$absb/$absb)));
+ return ($absb == 0 ? 0 : $absb*sqrt(1+$absa*$absa/$absb/$absb));
}
+1;
--- a/svbksb.pl
+++ b/svbksb.pl
@@ -1,21 +1,26 @@
#======================================================================
-# . / S V B K S B . P L
+# S V B K S B . P L
# doc: Sat Jul 31 22:47:03 1999
-# dlm: Sat Jul 31 23:06:40 1999
+# dlm: Wed Mar 11 19:29:39 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 30 32 NIL 0 0 72 2 2 4 ofnI
+# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# SVBKSB routine from Numerical Recipes adapted to ANTS
+#
+# solves Ax = b for x, given b
+#
+# Notes:
+# - A = U W V' as done in [svdcmp.pl]
# HISTORY:
# Jul 31, 1999: - manually converted from c-source
-
-# Notes:
-# - everything passed as refs
+# Mar 11, 2015: - fixed syntax errors (code had never been used before)
require "$ANTS/nrutil.pl";
+use strict;
+
sub svbksb($$$$$)
{
my($uR,$wR,$vR,$bR,$xR) = @_;
@@ -37,9 +42,9 @@
for ($j=1; $j<=$#{$wR}; $j++) {
$s = 0;
for ($jj=1; $jj<=$#{$wR}; $jj++) {
- $s += $vR->[$j][$jj] * tmp[$jj];
+ $s += $vR->[$j][$jj] * $tmp[$jj];
}
- $x->[$j] = $s;
+ $xR->[$j] = $s;
}
}
--- a/svdcmp.pl
+++ b/svdcmp.pl
@@ -1,15 +1,18 @@
#======================================================================
# S V D C M P . P L
# doc: Sun Aug 1 09:51:37 1999
-# dlm: Thu Jul 19 09:45:52 2001
+# dlm: Wed Mar 11 16:54:26 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 184 18 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 188 1 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# SVDCMP routine from Numerical Recipes adapted to ANTS
# HISTORY:
# Aug 01, 1999: - manually converted from c-source
+# Jul 19, 2001: - done *something* (message only?)
+# Mar 11, 2015: - picked up this never-used code because of need to fit circles
+
# Notes:
# - everything passed as refs
@@ -17,18 +20,21 @@
require "$ANTS/nrutil.pl";
require "$ANTS/pythag.pl";
+use strict;
+
sub svdcmp($$$)
-{
+{
my($aR,$wR,$vR) = @_; # params
my($flag,$i,$its,$j,$jj,$k,$l,$nm); # int
my($anorm,$c,$f,$g,$h,$s,$scale,$x,$y,$z); # float
my(@rv1); # float[]
vector(\@rv1,1,$#{$vR});
+ $g = $scale = $anorm = 0;
for ($i=1; $i<=$#{$vR}; $i++) {
$l = $i+1;
$rv1[$i] = $scale*$g;
- $g = 0; $s = 0; $scale = 0;
+ $g = $s = $scale = 0;
if ($i <= $#{$aR}) {
for ($k=$i; $k<=$#{$aR}; $k++) {
$scale += abs($aR->[$k][$i]);
@@ -87,7 +93,7 @@
}
}
}
- $anorm = &FMAX($anorm,(abs($wR->[$i])+abs($rv1[$i])));
+ $anorm = MAX($anorm,(abs($wR->[$i])+abs($rv1[$i])));
}
for ($i=$#{$vR}; $i>=1; $i--) {
if ($i < $#{$vR}) {
@@ -104,7 +110,7 @@
}
}
}
- for ($j=$l; $j<=$#{$vR; $j++) {
+ for ($j=$l; $j<=$#{$vR}; $j++) {
$vR->[$i][$j] = 0; $vR->[$j][$i] = 0;
}
}
@@ -112,7 +118,7 @@
$g = $rv1[$i];
$l = $i;
}
- for ($i=IMIN($#{$aR},$#{$vR}); $i>=1; $i--) {
+ for ($i=MIN($#{$aR},$#{$vR}); $i>=1; $i--) {
$l = $i+1;
$g = $wR->[$i];
for ($j=$l; $j<=$#{$vR}; $j++) {
@@ -146,9 +152,9 @@
$nm = $l-1;
if ((abs($rv1[$l])+$anorm) == $anorm) {
$flag = 0;
- break;
+ last;
}
- break if ((abs($wR->[$nm])+$anorm) == $anorm);
+ last if ($nm == 0) || ((abs($wR->[$nm])+$anorm) == $anorm); ## $nm == 0 test not in original code
}
if ($flag) {
$c = 0;
@@ -156,7 +162,7 @@
for ($i=$l; $i<=$k; $i++) {
$f = $s*$rv1[$i];
$rv1[$i] = $c*$rv1[$i];
- break if ((abs($f)+$anorm) == $anorm);
+ last if ((abs($f)+$anorm) == $anorm);
$g = $wR->[$i];
$h = &pythag($f,$g);
$wR->[$i] = $h;
@@ -179,7 +185,8 @@
$vR->[$j][$k] = -$vR->[$j][$k];
}
}
- break;
+# print(STDERR "its($k) = $its\n");
+ last;
}
croak("no convergence in 30 svdcmp iterations\n") if ($its == 30);
$x = $wR->[$l];