Click here to Skip to main content
Click here to Skip to main content
Add your own
alternative version

Yet Another Math Parser (YAMP)

, 30 Sep 2012
Constructing a fast math parser using Reflection to do numerics like Matlab.
YAMP_Demo.zip
YAMPCompare
LLMathParser
MathFormula
MathParserNet
Exceptions
MathParserTK
MathParser
YAMPCompare.pidb
YAMPConsole
YAMPConsole.pidb
YAMP
Exceptions
Expressions
Functions
ArgumentFunctions
LinearAlgebra
LogicFunctions
Spectroscopy
StandardFunctions
Trigonometric
Interfaces
Numerics
Decompositions
Others
Operators
AssigmentOperators
BinaryOperators
DotOperators
LogicOperators
UnaryOperators
Values
YAMP.csproj.user
YAMP.pidb
YAMP_Source.zip
YAMP.csproj.user
YAMP.pidb
using System;
using YAMP;
using System.Runtime.InteropServices;

namespace YAMP.Numerics
{
    /// <summary>
    /// Eigenvalues and eigenvectors of a real matrix. 
    /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
    /// diagonal and the eigenvector matrix V is orthogonal.
    /// I.e. A = V.Multiply(D.Multiply(V.Transpose())) and 
    /// V.Multiply(V.Transpose()) equals the identity matrix.
    /// If A is not symmetric, then the eigenvalue matrix D is block diagonal
    /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
    /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
    /// columns of V represent the eigenvectors in the sense that A*V = V*D,
    /// i.e. A.Multiply(V) equals V.Multiply(D).  The matrix V may be badly
    /// conditioned, or even singular, so the validity of the equation
    /// A = V*D*Inverse(V) depends upon V.cond().
    /// </summary>
    public class Eigenvalues
    {
        #region	 Members

        /// <summary>
        /// Row and column dimension (square matrix).
        /// </summary>
        int n;

        /// <summary>
        /// Symmetry flag.
        /// </summary>
        bool issymmetric;

        /// <summary>
        /// Arrays for internal storage of eigenvalues.
        /// </summary>
        double[] d, e;

        /// <summary>
        /// Array for internal storage of eigenvectors.
        /// </summary>
        double[][] V;

        /// <summary>
        /// Array for internal storage of nonsymmetric Hessenberg form.
        /// </summary>
        double[][] H;

        /// <summary>
        /// Working storage for nonsymmetric algorithm.
        /// </summary>
        double[] ort;

        #endregion //  Class variables

        #region Private Methods

        // Symmetric Householder reduction to tridiagonal form.
        void tred2()
        {
            //  This is derived from the Algol procedures tred2 by
            //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            for (int j = 0; j < n; j++)
                d[j] = V[n - 1][j];

            // Householder reduction to tridiagonal form.

            for (int i = n - 1; i > 0; i--)
            {
                // Scale to avoid under/overflow.

                var scale = 0.0;
                var h = 0.0;

                for (int k = 0; k < i; k++)
                    scale = scale + Math.Abs(d[k]);
                
                if (scale == 0.0)
                {
                    e[i] = d[i - 1];

                    for (int j = 0; j < i; j++)
                    {
                        d[j] = V[i - 1][j];
                        V[i][j] = 0.0;
                        V[j][i] = 0.0;
                    }
                }
                else
                {
                    // Generate Householder vector.

                    for (int k = 0; k < i; k++)
                    {
                        d[k] /= scale;
                        h += d[k] * d[k];
                    }

                    var f = d[i - 1];
                    var g = Math.Sqrt(h);

                    if (f > 0)
                        g = -g;
                    
                    e[i] = scale * g;
                    h = h - f * g;
                    d[i - 1] = f - g;

                    for (int j = 0; j < i; j++)
                        e[j] = 0.0;

                    // Apply similarity transformation to remaining columns.

                    for (int j = 0; j < i; j++)
                    {
                        f = d[j];
                        V[j][i] = f;
                        g = e[j] + V[j][j] * f;

                        for (int k = j + 1; k <= i - 1; k++)
                        {
                            g += V[k][j] * d[k];
                            e[k] += V[k][j] * f;
                        }

                        e[j] = g;
                    }

                    f = 0.0;

                    for (int j = 0; j < i; j++)
                    {
                        e[j] /= h;
                        f += e[j] * d[j];
                    }

                    var hh = f / (h + h);

                    for (int j = 0; j < i; j++)
                        e[j] -= hh * d[j];

                    for (int j = 0; j < i; j++)
                    {
                        f = d[j];
                        g = e[j];

                        for (int k = j; k <= i - 1; k++)
                            V[k][j] -= (f * e[k] + g * d[k]);

                        d[j] = V[i - 1][j];
                        V[i][j] = 0.0;
                    }
                }

                d[i] = h;
            }

            // Accumulate transformations.

            for (int i = 0; i < n - 1; i++)
            {
                V[n - 1][i] = V[i][i];
                V[i][i] = 1.0;
                var h = d[i + 1];

                if (h != 0.0)
                {
                    for (int k = 0; k <= i; k++)
                        d[k] = V[k][i + 1] / h;

                    for (int j = 0; j <= i; j++)
                    {
                        var g = 0.0;

                        for (int k = 0; k <= i; k++)
                            g += V[k][i + 1] * V[k][j];

                        for (int k = 0; k <= i; k++)
                            V[k][j] -= g * d[k];
                    }
                }

                for (int k = 0; k <= i; k++)
                    V[k][i + 1] = 0.0;
            }

            for (int j = 0; j < n; j++)
            {
                d[j] = V[n - 1][j];
                V[n - 1][j] = 0.0;
            }

            V[n - 1][n - 1] = 1.0;
            e[0] = 0.0;
        }

        // Symmetric tridiagonal QL algorithm.
        void tql2()
        {
            //  This is derived from the Algol procedures tql2, by
            //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            for (int i = 1; i < n; i++)
                e[i - 1] = e[i];

            e[n - 1] = 0.0;
            var f = 0.0;
            var tst1 = 0.0;
            var eps = Math.Pow(2.0, -52.0);

            for (int l = 0; l < n; l++)
            {
                // Find small subdiagonal element

                tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));
                var m = l;

                while (m < n)
                {
                    if (System.Math.Abs(e[m]) <= eps * tst1)
                        break;

                    m++;
                }

                // If m == l, d[l] is an eigenvalue,
                // otherwise, iterate.

                if (m > l)
                {
                    var iter = 0;

                    do
                    {
                        iter = iter + 1; // (Could check iteration count here.)

                        // Compute implicit shift

                        var g = d[l];
                        var p = (d[l + 1] - g) / (2.0 * e[l]);
                        var r = NumericHelpers.Hypot(p, 1.0);

                        if (p < 0)
                            r = -r;
                        
                        d[l] = e[l] / (p + r);
                        d[l + 1] = e[l] * (p + r);
                        var dl1 = d[l + 1];
                        var h = g - d[l];

                        for (int i = l + 2; i < n; i++)
                            d[i] -= h;
                        
                        f = f + h;

                        // Implicit QL transformation.

                        p = d[m];
                        var c = 1.0;
                        var c2 = c;
                        var c3 = c;
                        var el1 = e[l + 1];
                        var s = 0.0;
                        var s2 = 0.0;

                        for (int i = m - 1; i >= l; i--)
                        {
                            c3 = c2;
                            c2 = c;
                            s2 = s;
                            g = c * e[i];
                            h = c * p;
                            r = NumericHelpers.Hypot(p, e[i]);
                            e[i + 1] = s * r;
                            s = e[i] / r;
                            c = p / r;
                            p = c * d[i] - s * g;
                            d[i + 1] = h + s * (c * g + s * d[i]);

                            // Accumulate transformation.

                            for (int k = 0; k < n; k++)
                            {
                                h = V[k][i + 1];
                                V[k][i + 1] = s * V[k][i] + c * h;
                                V[k][i] = c * V[k][i] - s * h;
                            }
                        }

                        p = (-s) * s2 * c3 * el1 * e[l] / dl1;
                        e[l] = s * p;
                        d[l] = c * p;

                        // Check for convergence.
                    }
                    while (Math.Abs(e[l]) > eps * tst1);
                }

                d[l] = d[l] + f;
                e[l] = 0.0;
            }

            // Sort eigenvalues and corresponding vectors.

            for (int i = 0; i < n - 1; i++)
            {
                var k = i;
                var p = d[i];

                for (int j = i + 1; j < n; j++)
                {
                    if (d[j] < p)
                    {
                        k = j;
                        p = d[j];
                    }
                }

                if (k != i)
                {
                    d[k] = d[i];
                    d[i] = p;

                    for (int j = 0; j < n; j++)
                    {
                        p = V[j][i];
                        V[j][i] = V[j][k];
                        V[j][k] = p;
                    }
                }
            }
        }

        // Nonsymmetric reduction to Hessenberg form.
        void orthes()
        {
            //  This is derived from the Algol procedures orthes and ortran,
            //  by Martin and Wilkinson, Handbook for Auto. Comp.,
            //  Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutines in EISPACK.

            var low = 0;
            var high = n - 1;

            for (int m = low + 1; m <= high - 1; m++)
            {

                // Scale column.

                var scale = 0.0;

                for (int i = m; i <= high; i++)
                    scale = scale + Math.Abs(H[i][m - 1]);

                if (scale != 0.0)
                {
                    // Compute Householder transformation.

                    double h = 0.0;

                    for (int i = high; i >= m; i--)
                    {
                        ort[i] = H[i][m - 1] / scale;
                        h += ort[i] * ort[i];
                    }

                    double g = Math.Sqrt(h);

                    if (ort[m] > 0)
                        g = -g;
                    
                    h = h - ort[m] * g;
                    ort[m] = ort[m] - g;

                    // Apply Householder similarity transformation
                    // H = (I-u*u'/h)*H*(I-u*u')/h)

                    for (int j = m; j < n; j++)
                    {
                        double f = 0.0;

                        for (int i = high; i >= m; i--)
                            f += ort[i] * H[i][j];
                        
                        f = f / h;

                        for (int i = m; i <= high; i++)
                            H[i][j] -= f * ort[i];
                    }

                    for (int i = 0; i <= high; i++)
                    {
                        double f = 0.0;

                        for (int j = high; j >= m; j--)
                            f += ort[j] * H[i][j];

                        f = f / h;

                        for (int j = m; j <= high; j++)
                            H[i][j] -= f * ort[j];
                    }

                    ort[m] = scale * ort[m];
                    H[m][m - 1] = scale * g;
                }
            }

            // Accumulate transformations (Algol's ortran).

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    V[i][j] = (i == j ? 1.0 : 0.0);
            }

            for (int m = high - 1; m >= low + 1; m--)
            {
                if (H[m][m - 1] != 0.0)
                {
                    for (int i = m + 1; i <= high; i++)
                        ort[i] = H[i][m - 1];
                    
                    for (int j = m; j <= high; j++)
                    {
                        var g = 0.0;

                        for (int i = m; i <= high; i++)
                            g += ort[i] * V[i][j];
                        
                        // Double division avoids possible underflow
                        g = (g / ort[m]) / H[m][m - 1];

                        for (int i = m; i <= high; i++)
                            V[i][j] += g * ort[i];
                    }
                }
            }
        }


        // Complex scalar division.
        double cdivr, cdivi;
        void cdiv(double xr, double xi, double yr, double yi)
        {
            double r, d;

            if (Math.Abs(yr) > Math.Abs(yi))
            {
                r = yi / yr;
                d = yr + r * yi;
                cdivr = (xr + r * xi) / d;
                cdivi = (xi - r * xr) / d;
            }
            else
            {
                r = yr / yi;
                d = yi + r * yr;
                cdivr = (r * xr + xi) / d;
                cdivi = (r * xi - xr) / d;
            }
        }


        // Nonsymmetric reduction from Hessenberg to real Schur form.
        void hqr2()
        {
            //  This is derived from the Algol procedure hqr2,
            //  by Martin and Wilkinson, Handbook for Auto. Comp.,
            //  Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            // Initialize

            var nn = this.n;
            var n = nn - 1;
            var low = 0;
            var high = nn - 1;
            var eps = Math.Pow(2.0, -52.0);
            var exshift = 0.0;
            double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;

            // Store roots isolated by balanc and compute matrix norm

            var norm = 0.0;

            for (int i = 0; i < nn; i++)
            {
                if (i < low | i > high)
                {
                    d[i] = H[i][i];
                    e[i] = 0.0;
                }

                for (int j = Math.Max(i - 1, 0); j < nn; j++)
                    norm = norm + Math.Abs(H[i][j]);
            }

            // Outer loop over eigenvalue index

            var iter = 0;

            while (n >= low)
            {
                // Look for single small sub-diagonal element

                int l = n;

                while (l > low)
                {
                    s = Math.Abs(H[l - 1][l - 1]) + Math.Abs(H[l][l]);

                    if (s == 0.0)
                        s = norm;

                    if (Math.Abs(H[l][l - 1]) < eps * s)
                        break;

                    l--;
                }

                // Check for convergence
                // One root found

                if (l == n)
                {
                    H[n][n] = H[n][n] + exshift;
                    d[n] = H[n][n];
                    e[n] = 0.0;
                    n--;
                    iter = 0;

                    // Two roots found
                }
                else if (l == n - 1)
                {
                    w = H[n][n - 1] * H[n - 1][n];
                    p = (H[n - 1][n - 1] - H[n][n]) / 2.0;
                    q = p * p + w;
                    z = Math.Sqrt(Math.Abs(q));
                    H[n][n] = H[n][n] + exshift;
                    H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
                    x = H[n][n];

                    // Real pair

                    if (q >= 0)
                    {
                        if (p >= 0)
                            z = p + z;
                        else
                            z = p - z;

                        d[n - 1] = x + z;
                        d[n] = d[n - 1];

                        if (z != 0.0)
                            d[n] = x - w / z;
                        
                        e[n - 1] = 0.0;
                        e[n] = 0.0;
                        x = H[n][n - 1];
                        s = Math.Abs(x) + Math.Abs(z);
                        p = x / s;
                        q = z / s;
                        r = Math.Sqrt(p * p + q * q);
                        p = p / r;
                        q = q / r;

                        // Row modification

                        for (int j = n - 1; j < nn; j++)
                        {
                            z = H[n - 1][j];
                            H[n - 1][j] = q * z + p * H[n][j];
                            H[n][j] = q * H[n][j] - p * z;
                        }

                        // Column modification

                        for (int i = 0; i <= n; i++)
                        {
                            z = H[i][n - 1];
                            H[i][n - 1] = q * z + p * H[i][n];
                            H[i][n] = q * H[i][n] - p * z;
                        }

                        // Accumulate transformations

                        for (int i = low; i <= high; i++)
                        {
                            z = V[i][n - 1];
                            V[i][n - 1] = q * z + p * V[i][n];
                            V[i][n] = q * V[i][n] - p * z;
                        }

                        // Complex pair
                    }
                    else
                    {
                        d[n - 1] = x + p;
                        d[n] = x + p;
                        e[n - 1] = z;
                        e[n] = -z;
                    }

                    n = n - 2;
                    iter = 0;

                    // No convergence yet
                }
                else
                {
                    // Form shift

                    x = H[n][n];
                    y = 0.0;
                    w = 0.0;

                    if (l < n)
                    {
                        y = H[n - 1][n - 1];
                        w = H[n][n - 1] * H[n - 1][n];
                    }

                    // Wilkinson's original ad hoc shift

                    if (iter == 10)
                    {
                        exshift += x;

                        for (int i = low; i <= n; i++)
                            H[i][i] -= x;
                        
                        s = Math.Abs(H[n][n - 1]) + Math.Abs(H[n - 1][n - 2]);
                        x = y = 0.75 * s;
                        w = (-0.4375) * s * s;
                    }

                    // MATLAB's new ad hoc shift

                    if (iter == 30)
                    {
                        s = (y - x) / 2.0;
                        s = s * s + w;

                        if (s > 0)
                        {
                            s = Math.Sqrt(s);

                            if (y < x)
                                s = -s;

                            s = x - w / ((y - x) / 2.0 + s);

                            for (int i = low; i <= n; i++)
                                H[i][i] -= s;
                            
                            exshift += s;
                            x = y = w = 0.964;
                        }
                    }

                    iter = iter + 1; // (Could check iteration count here.)

                    // Look for two consecutive small sub-diagonal elements

                    int m = n - 2;

                    while (m >= l)
                    {
                        z = H[m][m];
                        r = x - z;
                        s = y - z;
                        p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
                        q = H[m + 1][m + 1] - z - r - s;
                        r = H[m + 2][m + 1];
                        s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
                        p = p / s;
                        q = q / s;
                        r = r / s;

                        if (m == l)
                            break;
                        
                        if (Math.Abs(H[m][m - 1]) * (Math.Abs(q) + Math.Abs(r)) < eps * (Math.Abs(p) * (Math.Abs(H[m - 1][m - 1]) + Math.Abs(z) + Math.Abs(H[m + 1][m + 1]))))
                            break;
                        
                        m--;
                    }

                    for (int i = m + 2; i <= n; i++)
                    {
                        H[i][i - 2] = 0.0;

                        if (i > m + 2)
                            H[i][i - 3] = 0.0;
                    }

                    // Double QR step involving rows l:n and columns m:n

                    for (int k = m; k <= n - 1; k++)
                    {
                        bool notlast = (k != n - 1);

                        if (k != m)
                        {
                            p = H[k][k - 1];
                            q = H[k + 1][k - 1];
                            r = (notlast ? H[k + 2][k - 1] : 0.0);
                            x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);

                            if (x != 0.0)
                            {
                                p = p / x;
                                q = q / x;
                                r = r / x;
                            }
                        }

                        if (x == 0.0)
                            break;

                        s = Math.Sqrt(p * p + q * q + r * r);

                        if (p < 0)
                            s = -s;
                        
                        if (s != 0)
                        {
                            if (k != m)
                                H[k][k - 1] = (-s) * x;
                            else if (l != m)
                                H[k][k - 1] = -H[k][k - 1];

                            p = p + s;
                            x = p / s;
                            y = q / s;
                            z = r / s;
                            q = q / p;
                            r = r / p;

                            // Row modification

                            for (int j = k; j < nn; j++)
                            {
                                p = H[k][j] + q * H[k + 1][j];

                                if (notlast)
                                {
                                    p = p + r * H[k + 2][j];
                                    H[k + 2][j] = H[k + 2][j] - p * z;
                                }

                                H[k][j] = H[k][j] - p * x;
                                H[k + 1][j] = H[k + 1][j] - p * y;
                            }

                            // Column modification

                            for (int i = 0; i <= Math.Min(n, k + 3); i++)
                            {
                                p = x * H[i][k] + y * H[i][k + 1];

                                if (notlast)
                                {
                                    p = p + z * H[i][k + 2];
                                    H[i][k + 2] = H[i][k + 2] - p * r;
                                }

                                H[i][k] = H[i][k] - p;
                                H[i][k + 1] = H[i][k + 1] - p * q;
                            }

                            // Accumulate transformations

                            for (int i = low; i <= high; i++)
                            {
                                p = x * V[i][k] + y * V[i][k + 1];

                                if (notlast)
                                {
                                    p = p + z * V[i][k + 2];
                                    V[i][k + 2] = V[i][k + 2] - p * r;
                                }

                                V[i][k] = V[i][k] - p;
                                V[i][k + 1] = V[i][k + 1] - p * q;
                            }
                        } // (s != 0)
                    } // k loop
                } // check convergence
            } // while (n >= low)

            // Backsubstitute to find vectors of upper triangular form

            if (norm == 0.0)
                return;

            for (n = nn - 1; n >= 0; n--)
            {
                p = d[n];
                q = e[n];

                // Real vector

                if (q == 0)
                {
                    int l = n;
                    H[n][n] = 1.0;

                    for (int i = n - 1; i >= 0; i--)
                    {
                        w = H[i][i] - p;
                        r = 0.0;

                        for (int j = l; j <= n; j++)
                            r = r + H[i][j] * H[j][n];

                        if (e[i] < 0.0)
                        {
                            z = w;
                            s = r;
                        }
                        else
                        {
                            l = i;

                            if (e[i] == 0.0)
                            {
                                if (w != 0.0)
                                    H[i][n] = (-r) / w;
                                else
                                    H[i][n] = (-r) / (eps * norm);

                                // Solve real equations
                            }
                            else
                            {
                                x = H[i][i + 1];
                                y = H[i + 1][i];
                                q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                                t = (x * s - z * r) / q;
                                H[i][n] = t;

                                if (Math.Abs(x) > Math.Abs(z))
                                    H[i + 1][n] = (-r - w * t) / x;
                                else
                                    H[i + 1][n] = (-s - y * t) / z;
                            }

                            // Overflow control

                            t = Math.Abs(H[i][n]);

                            if ((eps * t) * t > 1)
                            {
                                for (int j = i; j <= n; j++)
                                    H[j][n] = H[j][n] / t;
                            }
                        }
                    }

                    // Complex vector
                }
                else if (q < 0)
                {
                    int l = n - 1;

                    // Last vector component imaginary so matrix is triangular

                    if (Math.Abs(H[n][n - 1]) > Math.Abs(H[n - 1][n]))
                    {
                        H[n - 1][n - 1] = q / H[n][n - 1];
                        H[n - 1][n] = (-(H[n][n] - p)) / H[n][n - 1];
                    }
                    else
                    {
                        cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
                        H[n - 1][n - 1] = cdivr;
                        H[n - 1][n] = cdivi;
                    }

                    H[n][n - 1] = 0.0;
                    H[n][n] = 1.0;

                    for (int i = n - 2; i >= 0; i--)
                    {
                        double ra, sa, vr, vi;
                        ra = 0.0;
                        sa = 0.0;

                        for (int j = l; j <= n; j++)
                        {
                            ra = ra + H[i][j] * H[j][n - 1];
                            sa = sa + H[i][j] * H[j][n];
                        }

                        w = H[i][i] - p;

                        if (e[i] < 0.0)
                        {
                            z = w;
                            r = ra;
                            s = sa;
                        }
                        else
                        {
                            l = i;

                            if (e[i] == 0)
                            {
                                cdiv(-ra, -sa, w, q);
                                H[i][n - 1] = cdivr;
                                H[i][n] = cdivi;
                            }
                            else
                            {

                                // Solve complex equations

                                x = H[i][i + 1];
                                y = H[i + 1][i];
                                vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                                vi = (d[i] - p) * 2.0 * q;

                                if (vr == 0.0 & vi == 0.0)
                                    vr = eps * norm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));

                                cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
                                H[i][n - 1] = cdivr;
                                H[i][n] = cdivi;

                                if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q)))
                                {
                                    H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
                                    H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
                                }
                                else
                                {
                                    cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q);
                                    H[i + 1][n - 1] = cdivr;
                                    H[i + 1][n] = cdivi;
                                }
                            }

                            // Overflow control

                            t = Math.Max(Math.Abs(H[i][n - 1]), Math.Abs(H[i][n]));

                            if ((eps * t) * t > 1)
                            {
                                for (int j = i; j <= n; j++)
                                {
                                    H[j][n - 1] = H[j][n - 1] / t;
                                    H[j][n] = H[j][n] / t;
                                }
                            }
                        }
                    }
                }
            }

            // Vectors of isolated roots

            for (int i = 0; i < nn; i++)
            {
                if (i < low | i > high)
                {
                    for (int j = i; j < nn; j++)
                    {
                        V[i][j] = H[i][j];
                    }
                }
            }

            // Back transformation to get eigenvectors of original matrix

            for (int j = nn - 1; j >= low; j--)
            {
                for (int i = low; i <= high; i++)
                {
                    z = 0.0;
                    for (int k = low; k <= System.Math.Min(j, high); k++)
                    {
                        z = z + V[i][k] * H[k][j];
                    }
                    V[i][j] = z;
                }
            }
        }

        #endregion //  Private Methods

        #region Constructor

        /// <summary>
        /// Check for symmetry, then construct the eigenvalue decomposition
        /// </summary>
        /// <param name="Arg">Square matrix</param>
        /// <returns>Structure to access D and V.</returns>
        public Eigenvalues(MatrixValue Arg)
        {
            var A = Arg.GetRealArray();
            n = Arg.DimensionX;
            V = new double[n][];

            for (int i = 0; i < n; i++)
                V[i] = new double[n];
            
            d = new double[n];
            e = new double[n];

            issymmetric = true;

            for (int j = 0; (j < n) && issymmetric; j++)
            {
                for (int i = 0; (i < n) && issymmetric; i++)
                    issymmetric = (A[i][j] == A[j][i]);
            }

            if (issymmetric)
            {
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                        V[i][j] = A[i][j];
                }

                // Tridiagonalize.
                tred2();

                // Diagonalize.
                tql2();
            }
            else
            {
                H = new double[n][];

                for (int i2 = 0; i2 < n; i2++)
                    H[i2] = new double[n];

                ort = new double[n];

                for (int j = 0; j < n; j++)
                {
                    for (int i = 0; i < n; i++)
                        H[i][j] = A[i][j];
                }

                // Reduce to Hessenberg form.
                orthes();

                // Reduce Hessenberg to real Schur form.
                hqr2();
            }
        }

        #endregion //  Constructor

        #region Public Properties

        /// <summary>
        /// Return the real parts of the eigenvalues
        /// </summary>
        /// <returns>real(diag(D))</returns>
        virtual public double[] RealEigenvalues
        {
            get
            {
                return d;
            }
        }

        /// <summary>
        /// Return the imaginary parts of the eigenvalues
        /// </summary>
        /// <returns>imag(diag(D))</returns>
        virtual public double[] ImagEigenvalues
        {
            get
            {
                return e;
            }
        }

        /// <summary>
        /// Return the block diagonal eigenvalue matrix
        /// </summary>
        /// <returns>D</returns>
        virtual public MatrixValue D
        {
            get
            {
                var X = new MatrixValue(n, n);

                for (int i = 0; i < n; i++)
                    X[i, i] = new ScalarValue(d[i], e[i]);

                return X;
            }
        }

        #endregion //  Public Properties

        #region Public Methods

        /// <summary>
        /// Return the eigenvector matrix
        /// </summary>
        /// <returns>V</returns>
        public virtual MatrixValue GetV()
        {
            return new MatrixValue(V, n, n);
        }

        #endregion //  Public Methods
    }
}

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

About the Author

Florian Rappl
Chief Technology Officer
Germany Germany
Florian is from Regensburg, Germany. He started his programming career with Perl. After programming C/C++ for some years he discovered his favorite programming language C#. He did work at Siemens as a programmer until he decided to study Physics. During his studies he worked as an IT consultant for various companies.
 
Florian is also giving lectures in C#, HTML5 with CSS3 and JavaScript, and other topics. Having graduated from University with a Master's degree in theoretical physics he is currently busy doing his PhD in the field of High Performance Computing.
Follow on   Google+

| Advertise | Privacy | Mobile
Web04 | 2.8.140721.1 | Last Updated 30 Sep 2012
Article Copyright 2012 by Florian Rappl
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid