1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Computes all eigenvalues and eigenvectors of
    /// a real symmetric matrix by Jacobi's method.
    /// </summary>
    public class Jacobi
    {
        private int n { get; set; }
        private double[,] a;
        private double[,] v;
        private double[] d;
        private int nrot { get; set; }
        private double EPS { get; set; }

        /// <summary>
        /// Computes all eigenvalues and eigenvectors of a real symmetric matrix
        /// a[0..n - 1][0..n - 1]. On output, d[0..n - 1] contains the eigenvalues of a
        /// sorted into descending order, while v[0..n - 1][0..n - 1] is a matrix whose
        /// columns contain the corresponding normalized eigenvectors.nrot contains
        /// the number of Jacobi rotations that were required.Only the upper triangle
        /// of a is accessed.
        /// </summary>
        /// <param name="aa"></param>
        /// <exception cref="Exception"></exception>
        public Jacobi(double[,] aa)
        {
            this.n = aa.GetLength(0);
            this.a = aa;
            this.v = new double[n, n];
            this.d = new double[n];
            this.nrot = 0;
            this.EPS = float.Epsilon;

            double[] b = new double[n];
            double[] z = new double[n];

            for (int ip = 0; ip < n; ip++)
            {
                for (int iq = 0; iq < n; iq++)
                {
                    v[ip, iq] = 0.0;
                }
                v[ip, ip] = 1.0;
            }
            for (int ip = 0; ip < n; ip++)
            {
                b[ip] = d[ip] = a[ip, ip];
                z[ip] = 0.0;
            }
            for (int i = 1; i <= 50; i++)
            {
                double sm = 0.0;
                for (int ip = 0; ip < n - 1; ip++)
                {
                    for (int iq = ip + 1; iq < n; iq++)
                    {
                        sm += Math.Abs(a[ip, iq]);
                    }
                }
                //if (sm == 0.0)
                if (Math.Abs(sm) <= float.Epsilon)
                {
                    eigsrt( d,  v);
                    return;
                }
                double tresh;
                if (i < 4)
                {
                    tresh = 0.2 * sm / (n * n);
                }
                else
                {
                    tresh = 0.0;
                }
                for (int ip = 0; ip < n - 1; ip++)
                {
                    for (int iq = ip + 1; iq < n; iq++)
                    {
                        double g = 100.0 * Math.Abs(a[ip, iq]);
                        if (i > 4 && g <= EPS * Math.Abs(d[ip]) && g <= EPS * Math.Abs(d[iq]))
                        {
                            a[ip, iq] = 0.0;
                        }
                        else if (Math.Abs(a[ip, iq]) > tresh)
                        {
                            double h = d[iq] - d[ip];
                            double t;
                            if (g <= EPS * Math.Abs(h))
                            {
                                t = (a[ip, iq]) / h;
                            }
                            else
                            {
                                double theta = 0.5 * h / (a[ip, iq]);
                                t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0 + theta * theta));
                                if (theta < 0.0)
                                {
                                    t = -t;
                                }
                            }
                            double c = 1.0 / Math.Sqrt(1 + t * t);
                            double s = t * c;
                            double tau = s / (1.0 + c);
                            h = t * a[ip, iq];
                            z[ip] -= h;
                            z[iq] += h;
                            d[ip] -= h;
                            d[iq] += h;
                            a[ip, iq] = 0.0;
                            for (int j = 0; j < ip; j++)
                            {
                                rot( a, s, tau, j, ip, j, iq);
                            }
                            for (int j = ip + 1; j < iq; j++)
                            {
                                rot( a, s, tau, ip, j, j, iq);
                            }
                            for (int j = iq + 1; j < n; j++)
                            {
                                rot( a, s, tau, ip, j, iq, j);
                            }
                            for (int j = 0; j < n; j++)
                            {
                                rot( v, s, tau, j, ip, j, iq);
                            }
                            ++nrot;
                        }
                    }
                }
                for (int ip = 0; ip < n; ip++)
                {
                    b[ip] += z[ip];
                    d[ip] = b[ip];
                    z[ip] = 0.0;
                }
            }
            throw new Exception("Too many iterations in routine jacobi");
        }

        public void rot(double[,] a, double s, double tau, int i, int j, int k, int l)
        {
            double g = a[i, j];
            double h = a[k, l];
            a[i, j] = g - s * (h + g * tau);
            a[k, l] = h + s * (g - h * tau);
        }

        /// <summary>
        /// Given the eigenvalues d[0..n - 1] and(optionally) the eigenvectors
        /// v[0..n - 1][0..n - 1] as determined by Jacobi or tqli, this routine sorts the
        /// eigenvalues into descending order and rearranges the columns of v
        /// correspondingly.The method is straight insertion.
        /// </summary>
        /// <param name="d"></param>
        /// <param name="v"></param>
        public static void eigsrt(double[] d, double[,] v)
        {
            int k;
            int n = d.Length;
            for (int i = 0; i < n - 1; i++)
            {
                double p = d[k = i];
                for (int j = i; j < n; j++)
                {
                    if (d[j] >= p)
                    {
                        p = d[k = j];
                    }
                }
                if (k != i)
                {
                    d[k] = d[i];
                    d[i] = p;
                    if (v != null)
                    {
                        for (int j = 0; j < n; j++)
                        {
                            p = v[j, i];
                            v[j, i] = v[j, k];
                            v[j, k] = p;
                        }
                    }
                }
            }
        }
    }
}
 

2 代码格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Computes all eigenvalues and eigenvectors of
    /// a real symmetric matrix by Jacobi's method.
    /// </summary>
    public class Jacobi
    {
        private int n { get; set; }
        private double[,] a;
        private double[,] v;
        private double[] d;
        private int nrot { get; set; }
        private double EPS { get; set; }


        /// <summary>
        /// Computes all eigenvalues and eigenvectors of a real symmetric matrix
        /// a[0..n - 1][0..n - 1]. On output, d[0..n - 1] contains the eigenvalues of a
        /// sorted into descending order, while v[0..n - 1][0..n - 1] is a matrix whose
        /// columns contain the corresponding normalized eigenvectors.nrot contains
        /// the number of Jacobi rotations that were required.Only the upper triangle
        /// of a is accessed.
        /// </summary>
        /// <param name="aa"></param>
        /// <exception cref="Exception"></exception>
        public Jacobi(double[,] aa)
        {
            this.n = aa.GetLength(0);
            this.a = aa;
            this.v = new double[n, n];
            this.d = new double[n];
            this.nrot = 0;
            this.EPS = float.Epsilon;

            double[] b = new double[n];
            double[] z = new double[n];

            for (int ip = 0; ip < n; ip++)
            {
                for (int iq = 0; iq < n; iq++)
                {
                    v[ip, iq] = 0.0;
                }
                v[ip, ip] = 1.0;
            }
            for (int ip = 0; ip < n; ip++)
            {
                b[ip] = d[ip] = a[ip, ip];
                z[ip] = 0.0;
            }
            for (int i = 1; i <= 50; i++)
            {
                double sm = 0.0;
                for (int ip = 0; ip < n - 1; ip++)
                {
                    for (int iq = ip + 1; iq < n; iq++)
                    {
                        sm += Math.Abs(a[ip, iq]);
                    }
                }
                //if (sm == 0.0)
                if (Math.Abs(sm) <= float.Epsilon)
                {
                    eigsrt( d,  v);
                    return;
                }
                double tresh;
                if (i < 4)
                {
                    tresh = 0.2 * sm / (n * n);
                }
                else
                {
                    tresh = 0.0;
                }
                for (int ip = 0; ip < n - 1; ip++)
                {
                    for (int iq = ip + 1; iq < n; iq++)
                    {
                        double g = 100.0 * Math.Abs(a[ip, iq]);
                        if (i > 4 && g <= EPS * Math.Abs(d[ip]) && g <= EPS * Math.Abs(d[iq]))
                        {
                            a[ip, iq] = 0.0;
                        }
                        else if (Math.Abs(a[ip, iq]) > tresh)
                        {
                            double h = d[iq] - d[ip];
                            double t;
                            if (g <= EPS * Math.Abs(h))
                            {
                                t = (a[ip, iq]) / h;
                            }
                            else
                            {
                                double theta = 0.5 * h / (a[ip, iq]);
                                t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0 + theta * theta));
                                if (theta < 0.0)
                                {
                                    t = -t;
                                }
                            }
                            double c = 1.0 / Math.Sqrt(1 + t * t);
                            double s = t * c;
                            double tau = s / (1.0 + c);
                            h = t * a[ip, iq];
                            z[ip] -= h;
                            z[iq] += h;
                            d[ip] -= h;
                            d[iq] += h;
                            a[ip, iq] = 0.0;
                            for (int j = 0; j < ip; j++)
                            {
                                rot( a, s, tau, j, ip, j, iq);
                            }
                            for (int j = ip + 1; j < iq; j++)
                            {
                                rot( a, s, tau, ip, j, j, iq);
                            }
                            for (int j = iq + 1; j < n; j++)
                            {
                                rot( a, s, tau, ip, j, iq, j);
                            }
                            for (int j = 0; j < n; j++)
                            {
                                rot( v, s, tau, j, ip, j, iq);
                            }
                            ++nrot;
                        }
                    }
                }
                for (int ip = 0; ip < n; ip++)
                {
                    b[ip] += z[ip];
                    d[ip] = b[ip];
                    z[ip] = 0.0;
                }
            }
            throw new Exception("Too many iterations in routine jacobi");
        }

        public void rot(double[,] a, double s, double tau, int i, int j, int k, int l)
        {
            double g = a[i, j];
            double h = a[k, l];
            a[i, j] = g - s * (h + g * tau);
            a[k, l] = h + s * (g - h * tau);
        }

        /// <summary>
        /// Given the eigenvalues d[0..n - 1] and(optionally) the eigenvectors
        /// v[0..n - 1][0..n - 1] as determined by Jacobi or tqli, this routine sorts the
        /// eigenvalues into descending order and rearranges the columns of v
        /// correspondingly.The method is straight insertion.
        /// </summary>
        /// <param name="d"></param>
        /// <param name="v"></param>
        public static void eigsrt(double[] d, double[,] v)
        {
            int k;
            int n = d.Length;
            for (int i = 0; i < n - 1; i++)
            {
                double p = d[k = i];
                for (int j = i; j < n; j++)
                {
                    if (d[j] >= p)
                    {
                        p = d[k = j];
                    }
                }
                if (k != i)
                {
                    d[k] = d[i];
                    d[i] = p;
                    if (v != null)
                    {
                        for (int j = 0; j < n; j++)
                        {
                            p = v[j, i];
                            v[j, i] = v[j, k];
                            v[j, k] = p;
                        }
                    }
                }
            }
        }
    }
}