Automatic Image Stitching with Accord.NET

��// Accord Math Library

// Accord.NET framework

// http://www.crsouza.com

//

// Copyright � C�sar Souza, 2009-2010

// cesarsouza at gmail.com

//



namespace Accord.Math.Decompositions

{

    using System;

    using Accord.Math;



    /// <summary>

    ///   Determines the Generalized eigenvalues and eigenvectors of two real square matrices.

    /// </summary>

    /// <remarks>

    ///   <para>

    ///     A generalized eigenvalue problem is the problem of finding a vector <c>v</c> that

    ///     obeys <c>A * v = � * B * v</c> where <c>A</c> and <c>B</c> are matrices. If <c>v</c>

    ///     obeys this equation, with some <c>�</c>, then we call <c>v</c> the generalized eigenvector

    ///     of <c>A</c> and <c>B</c>, and <c>�</c> is called the generalized eigenvalue of <c>A</c>

    ///     and <c>B</c> which corresponds to the generalized eigenvector <c>v</c>. The possible

    ///     values of <c>�</c>, must obey the identity <c>det(A - �*B) = 0</c>.</para>

    ///   <para>

    ///     Part of this code has been adapted from the original EISPACK routines in Fortran.</para>

    ///  

    ///   <para>

    ///     References:

    ///     <list type="bullet">

    ///       <item><description>

    ///         <a href="http://en.wikipedia.org/wiki/Generalized_eigenvalue_problem#Generalized_eigenvalue_problem">

    ///         http://en.wikipedia.org/wiki/Generalized_eigenvalue_problem#Generalized_eigenvalue_problem</a>

    ///       </description></item>

    ///       <item><description>

    ///         <a href="http://www.netlib.org/eispack/">

    ///         http://www.netlib.org/eispack/</a>

    ///       </description></item>

    ///     </list>

    ///   </para>

    /// </remarks>

    public sealed class GeneralizedEigenvalueDecomposition

    {

        private int n;

        private double[] ar;

        private double[] ai;

        private double[] beta;

        private double[,] Z;





        /// <summary>Construct a generalized eigenvalue decomposition.</summary>

        /// <param name="a">The first matrix of the (A,B) matrix pencil.</param>

        /// <param name="b">The second matrix of the (A,B) matrix pencil.</param>

        public GeneralizedEigenvalueDecomposition(double[,] a, double[,] b)

        {

            if (a == null) 

                throw new ArgumentNullException("a", "Matrix A cannot be null.");



            if (b == null)

                throw new ArgumentNullException("b", "Matrix B cannot be null.");



            if (a.GetLength(0) != a.GetLength(1))

                throw new ArgumentException("Matrix is not a square matrix.", "a");



            if (b.GetLength(0) != b.GetLength(1))

                throw new ArgumentException("Matrix is not a square matrix.", "b");



            if (a.GetLength(0) != b.GetLength(0) || a.GetLength(1) != b.GetLength(1))

                throw new ArgumentException("Matrix dimensions do not match", "b");





            n = a.GetLength(0);

            ar = new double[n];

            ai = new double[n];

            beta = new double[n];

            Z = new double[n, n];

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

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



            bool matz = true;

            int ierr = 0;





            // reduces A to upper hessenberg form and B to upper

            // triangular form using orthogonal transformations

            qzhes(n, A, B, matz, Z);



            // reduces the hessenberg matrix A to quasi-triangular form

            // using orthogonal transformations while maintaining the

            // triangular form of the B matrix.

            qzit(n, A, B, Special.DoubleEpsilon, matz, Z, ref ierr);



            // 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.

            qzval(n, A, B, ar, ai, beta, matz, Z);



            // computes the eigenvectors of the triangular problem and

            // transforms the results back to the original coordinate system.

            qzvec(n, A, B, ar, ai, beta, Z);



        }





        /// <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



    }

}