C#，数值计算——计算实对称矩阵所有特征值与特征向量的三角分解与QL迭代法源程序

1 文本格式

using System;

namespace Legalsoft.Truffer
{
/// <summary>
/// Computes all eigenvalues and eigenvectors of a real symmetric matrix by
/// reduction to tridiagonal form followed by QL iteration.
/// </summary>
public class Symmeig
{
public int n { get; set; }
public double[,] z;
public double[] d;
public double[] e;
public bool yesvecs { get; set; }

/// <summary>
/// Computes all eigenvalues and eigenvectors of a real symmetric matrix
/// a[0..n - 1][0..n - 1] by reduction to tridiagonal form followed by QL
/// iteration.On output, d[0..n - 1] contains the eigenvalues of a sorted into
/// descending order, while z[0..n - 1][0..n - 1] is a matrix whose columns contain
/// the corresponding normalized eigenvectors.If yesvecs is input as true (the
/// default), then the eigenvectors are computed.If yesvecs is input as false,
/// only the eigenvalues are computed.
/// </summary>
/// <param name="a"></param>
/// <param name="yesvec"></param>
public Symmeig(double[,] a, bool yesvec = true)
{
this.n = a.GetLength(0);
this.z = a;
this.d = new double[n];
this.e = new double[n];
this.yesvecs = yesvec;

tred2();
tqli();
sort();
}

/// <summary>
/// Computes all eigenvalues and(optionally) eigenvectors of a real,
/// symmetric, tridiagonal matrix by QL iteration.On input, dd[0..n - 1]
/// contains the diagonal elements of the tridi- agonal matrix.The vector
/// ee[0..n-1] inputs the subdiagonal elements of the tridiagonal matrix, with
/// ee[0] arbitrary.Output is the same as the constructor above.
/// </summary>
/// <param name="dd"></param>
/// <param name="ee"></param>
/// <param name="yesvec"></param>
public Symmeig(double[] dd, double[] ee, bool yesvec = true)
{
this.n = dd.Length;
this.d = dd;
this.e = ee;
this.z = new double[n, n];
this.yesvecs = yesvec;
for (int i = 0; i < n; i++)
{
z[i, i] = 1.0;
}

tqli();
sort();
}

public void sort()
{
if (yesvecs)
{
Jacobi.eigsrt( d, z);
}
else
{
Jacobi.eigsrt( d, z);
}
}

/// <summary>
/// Householder reduction of a real symmetric matrix z[0..n - 1][0..n-1]. (The
/// input matrix A to Symmeig is stored in z.) On output, z is replaced by the
/// orthogonal matrix Q effecting the transformation.d[0..n - 1] contains the
/// diagonal elements of the tridiagonal matrix and e[0..n - 1] the off-diagonal
/// elements, with e[0]=0. If yesvecs is false, so that only eigenvalues will
/// subsequently be determined, several statements are omitted, in which case z
/// contains no useful information on output.
/// </summary>
public void tred2()
{
for (int i = n - 1; i > 0; i--)
{
int l = i - 1;
double h = 0.0;
double scale = 0.0;
if (l > 0)
{
for (int k = 0; k < i; k++)
{
scale += Math.Abs(z[i, k]);
}
//if (scale == 0.0)
if (Math.Abs(scale) <= float.Epsilon)
{
e[i] = z[i, l];
}
else
{
for (int k = 0; k < i; k++)
{
z[i, k] /= scale;
h += z[i, k] * z[i, k];
}
double f = z[i, l];
double g = (f >= 0.0 ? -Math.Sqrt(h) : Math.Sqrt(h));
e[i] = scale * g;
h -= f * g;
z[i, l] = f - g;
f = 0.0;
for (int j = 0; j < i; j++)
{
if (yesvecs)
{
z[j, i] = z[i, j] / h;
}
g = 0.0;
for (int k = 0; k < j + 1; k++)
{
g += z[j, k] * z[i, k];
}
for (int k = j + 1; k < i; k++)
{
g += z[k, j] * z[i, k];
}
e[j] = g / h;
f += e[j] * z[i, j];
}
double hh = f / (h + h);
for (int j = 0; j < i; j++)
{
f = z[i, j];
e[j] = g = e[j] - hh * f;
for (int k = 0; k < j + 1; k++)
{
z[j, k] -= (f * e[k] + g * z[i, k]);
}
}
}
}
else
{
e[i] = z[i, l];
}
d[i] = h;
}
if (yesvecs)
{
d[0] = 0.0;
}
e[0] = 0.0;
for (int i = 0; i < n; i++)
{
if (yesvecs)
{
if (d[i] != 0.0)
{
for (int j = 0; j < i; j++)
{
double g = 0.0;
for (int k = 0; k < i; k++)
{
g += z[i, k] * z[k, j];
}
for (int k = 0; k < i; k++)
{
z[k, j] -= g * z[k, i];
}
}
}
d[i] = z[i, i];
z[i, i] = 1.0;
for (int j = 0; j < i; j++)
{
z[j, i] = z[i, j] = 0.0;
}
}
else
{
d[i] = z[i, i];
}
}
}

/// <summary>
/// QL algorithm with implicit shifts to determine the eigenvalues and
/// (optionally) the eigenvectors of a real, symmetric, tridiagonal matrix, or
/// of a real symmetric matrix previously reduced by tred2. On input,
/// d[0..n-1] contains the diagonal elements of the tridiagonal matrix. On
/// output, it returns the eigenvalues. The vector e[0..n - 1] inputs the
/// subdiagonal elements of the tridiagonal matrix, with e[0] arbitrary. On
/// output e is destroyed.If the eigenvectors of a tridiagonal matrix are
/// desired, the matrix z[0..n - 1][0..n - 1] is input as the identity matrix.If
/// the eigenvectors of a matrix that has been reduced by tred2 are required,
/// then z is input as the matrix output by tred2.In either case, column k of
/// z returns the normalized eigenvector corresponding to d[k]
/// </summary>
/// <exception cref="Exception"></exception>
public void tqli()
{
const double EPS = float.Epsilon;
for (int i = 1; i < n; i++)
{
e[i - 1] = e[i];
}
e[n - 1] = 0.0;
for (int l = 0; l < n; l++)
{
int iter = 0;
int m;
do
{
for (m = l; m < n - 1; m++)
{
double dd = Math.Abs(d[m]) + Math.Abs(d[m + 1]);
if (Math.Abs(e[m]) <= EPS * dd)
{
break;
}
}
if (m != l)
{
if (iter++ == 30)
{
throw new Exception("Too many iterations in tqli");
}
double g = (d[l + 1] - d[l]) / (2.0 * e[l]);
double r = pythag(g, 1.0);
g = d[m] - d[l] + e[l] / (g + Globals.SIGN(r, g));
double s = 1.0;
double c = 1.0;
double p = 0.0;
int i = m - 1;
for (; i >= l; i--)
{
double f = s * e[i];
double b = c * e[i];
e[i + 1] = (r = pythag(f, g));
//if (r == 0.0)
if (Math.Abs(r) <= float.Epsilon)
{
d[i + 1] -= p;
e[m] = 0.0;
break;
}
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = (d[i] - g) * s + 2.0 * c * b;
d[i + 1] = g + (p = s * r);
g = c * r - b;
if (yesvecs)
{
for (int k = 0; k < n; k++)
{
f = z[k, i + 1];
z[k, i + 1] = s * z[k, i] + c * f;
z[k, i] = c * z[k, i] - s * f;
}
}
}
//if (r == 0.0 && i >= l)
if (Math.Abs(r) <= float.Epsilon && i >= l)
{
continue;
}
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while (m != l);
}
}

public double pythag(double a, double b)
{
double absa = Math.Abs(a);
double absb = Math.Abs(b);
//return (absa > absb ? absa * Math.Sqrt(1.0 + Globals.SQR(absb / absa)) : (absb == 0.0 ? 0.0 : absb * Math.Sqrt(1.0 + Globals.SQR(absa / absb))));
return (absa > absb ? absa * Math.Sqrt(1.0 + Globals.SQR(absb / absa)) : ((absb <= float.Epsilon) ? 0.0 : absb * Math.Sqrt(1.0 + Globals.SQR(absa / absb))));
}
}
}

2 代码格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Computes all eigenvalues and eigenvectors of a real symmetric matrix by
    /// reduction to tridiagonal form followed by QL iteration.
    /// </summary>
    public class Symmeig
    {
        public int n { get; set; }
        public double[,] z;
        public double[] d;
        public double[] e;
        public bool yesvecs { get; set; }

        /// <summary>
        /// Computes all eigenvalues and eigenvectors of a real symmetric matrix
        /// a[0..n - 1][0..n - 1] by reduction to tridiagonal form followed by QL
        /// iteration.On output, d[0..n - 1] contains the eigenvalues of a sorted into
        /// descending order, while z[0..n - 1][0..n - 1] is a matrix whose columns contain
        /// the corresponding normalized eigenvectors.If yesvecs is input as true (the
        /// default), then the eigenvectors are computed.If yesvecs is input as false,
        /// only the eigenvalues are computed.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="yesvec"></param>
        public Symmeig(double[,] a, bool yesvec = true)
        {
            this.n = a.GetLength(0);
            this.z = a;
            this.d = new double[n];
            this.e = new double[n];
            this.yesvecs = yesvec;

            tred2();
            tqli();
            sort();
        }

        /// <summary>
        /// Computes all eigenvalues and(optionally) eigenvectors of a real,
        /// symmetric, tridiagonal matrix by QL iteration.On input, dd[0..n - 1]
        /// contains the diagonal elements of the tridi- agonal matrix.The vector
        /// ee[0..n-1] inputs the subdiagonal elements of the tridiagonal matrix, with
        /// ee[0] arbitrary.Output is the same as the constructor above.
        /// </summary>
        /// <param name="dd"></param>
        /// <param name="ee"></param>
        /// <param name="yesvec"></param>
        public Symmeig(double[] dd, double[] ee, bool yesvec = true)
        {
            this.n = dd.Length;
            this.d = dd;
            this.e = ee;
            this.z = new double[n, n];
            this.yesvecs = yesvec;
            for (int i = 0; i < n; i++)
            {
                z[i, i] = 1.0;
            }

            tqli();
            sort();
        }

        public void sort()
        {
            if (yesvecs)
            {
                Jacobi.eigsrt( d,  z);
            }
            else
            {
                Jacobi.eigsrt( d,  z);
            }
        }

        /// <summary>
        /// Householder reduction of a real symmetric matrix z[0..n - 1][0..n-1]. (The
        /// input matrix A to Symmeig is stored in z.) On output, z is replaced by the
        /// orthogonal matrix Q effecting the transformation.d[0..n - 1] contains the
        /// diagonal elements of the tridiagonal matrix and e[0..n - 1] the off-diagonal
        /// elements, with e[0]=0. If yesvecs is false, so that only eigenvalues will
        /// subsequently be determined, several statements are omitted, in which case z
        /// contains no useful information on output.
        /// </summary>
        public void tred2()
        {
            for (int i = n - 1; i > 0; i--)
            {
                int l = i - 1;
                double h = 0.0;
                double scale = 0.0;
                if (l > 0)
                {
                    for (int k = 0; k < i; k++)
                    {
                        scale += Math.Abs(z[i, k]);
                    }
                    //if (scale == 0.0)
                    if (Math.Abs(scale) <= float.Epsilon)
                    {
                        e[i] = z[i, l];
                    }
                    else
                    {
                        for (int k = 0; k < i; k++)
                        {
                            z[i, k] /= scale;
                            h += z[i, k] * z[i, k];
                        }
                        double f = z[i, l];
                        double g = (f >= 0.0 ? -Math.Sqrt(h) : Math.Sqrt(h));
                        e[i] = scale * g;
                        h -= f * g;
                        z[i, l] = f - g;
                        f = 0.0;
                        for (int j = 0; j < i; j++)
                        {
                            if (yesvecs)
                            {
                                z[j, i] = z[i, j] / h;
                            }
                            g = 0.0;
                            for (int k = 0; k < j + 1; k++)
                            {
                                g += z[j, k] * z[i, k];
                            }
                            for (int k = j + 1; k < i; k++)
                            {
                                g += z[k, j] * z[i, k];
                            }
                            e[j] = g / h;
                            f += e[j] * z[i, j];
                        }
                        double hh = f / (h + h);
                        for (int j = 0; j < i; j++)
                        {
                            f = z[i, j];
                            e[j] = g = e[j] - hh * f;
                            for (int k = 0; k < j + 1; k++)
                            {
                                z[j, k] -= (f * e[k] + g * z[i, k]);
                            }
                        }
                    }
                }
                else
                {
                    e[i] = z[i, l];
                }
                d[i] = h;
            }
            if (yesvecs)
            {
                d[0] = 0.0;
            }
            e[0] = 0.0;
            for (int i = 0; i < n; i++)
            {
                if (yesvecs)
                {
                    if (d[i] != 0.0)
                    {
                        for (int j = 0; j < i; j++)
                        {
                            double g = 0.0;
                            for (int k = 0; k < i; k++)
                            {
                                g += z[i, k] * z[k, j];
                            }
                            for (int k = 0; k < i; k++)
                            {
                                z[k, j] -= g * z[k, i];
                            }
                        }
                    }
                    d[i] = z[i, i];
                    z[i, i] = 1.0;
                    for (int j = 0; j < i; j++)
                    {
                        z[j, i] = z[i, j] = 0.0;
                    }
                }
                else
                {
                    d[i] = z[i, i];
                }
            }
        }

        /// <summary>
        /// QL algorithm with implicit shifts to determine the eigenvalues and
        /// (optionally) the eigenvectors of a real, symmetric, tridiagonal matrix, or
        /// of a real symmetric matrix previously reduced by tred2. On input,
        /// d[0..n-1] contains the diagonal elements of the tridiagonal matrix. On
        /// output, it returns the eigenvalues. The vector e[0..n - 1] inputs the
        /// subdiagonal elements of the tridiagonal matrix, with e[0] arbitrary. On
        /// output e is destroyed.If the eigenvectors of a tridiagonal matrix are
        /// desired, the matrix z[0..n - 1][0..n - 1] is input as the identity matrix.If
        /// the eigenvectors of a matrix that has been reduced by tred2 are required,
        /// then z is input as the matrix output by tred2.In either case, column k of
        /// z returns the normalized eigenvector corresponding to d[k]
        /// </summary>
        /// <exception cref="Exception"></exception>
        public void tqli()
        {
            const double EPS = float.Epsilon;
            for (int i = 1; i < n; i++)
            {
                e[i - 1] = e[i];
            }
            e[n - 1] = 0.0;
            for (int l = 0; l < n; l++)
            {
                int iter = 0;
                int m;
                do
                {
                    for (m = l; m < n - 1; m++)
                    {
                        double dd = Math.Abs(d[m]) + Math.Abs(d[m + 1]);
                        if (Math.Abs(e[m]) <= EPS * dd)
                        {
                            break;
                        }
                    }
                    if (m != l)
                    {
                        if (iter++ == 30)
                        {
                            throw new Exception("Too many iterations in tqli");
                        }
                        double g = (d[l + 1] - d[l]) / (2.0 * e[l]);
                        double r = pythag(g, 1.0);
                        g = d[m] - d[l] + e[l] / (g + Globals.SIGN(r, g));
                        double s = 1.0;
                        double c = 1.0;
                        double p = 0.0;
                        int i = m - 1;
                        for (; i >= l; i--)
                        {
                            double f = s * e[i];
                            double b = c * e[i];
                            e[i + 1] = (r = pythag(f, g));
                            //if (r == 0.0)
                            if (Math.Abs(r) <= float.Epsilon)
                            {
                                d[i + 1] -= p;
                                e[m] = 0.0;
                                break;
                            }
                            s = f / r;
                            c = g / r;
                            g = d[i + 1] - p;
                            r = (d[i] - g) * s + 2.0 * c * b;
                            d[i + 1] = g + (p = s * r);
                            g = c * r - b;
                            if (yesvecs)
                            {
                                for (int k = 0; k < n; k++)
                                {
                                    f = z[k, i + 1];
                                    z[k, i + 1] = s * z[k, i] + c * f;
                                    z[k, i] = c * z[k, i] - s * f;
                                }
                            }
                        }
                        //if (r == 0.0 && i >= l)
                        if (Math.Abs(r) <= float.Epsilon && i >= l)
                        {
                            continue;
                        }
                        d[l] -= p;
                        e[l] = g;
                        e[m] = 0.0;
                    }
                } while (m != l);
            }
        }

        public double pythag(double a, double b)
        {
            double absa = Math.Abs(a);
            double absb = Math.Abs(b);
            //return (absa > absb ? absa * Math.Sqrt(1.0 + Globals.SQR(absb / absa)) : (absb == 0.0 ? 0.0 : absb * Math.Sqrt(1.0 + Globals.SQR(absa / absb))));
            return (absa > absb ? absa * Math.Sqrt(1.0 + Globals.SQR(absb / absa)) : ((absb <= float.Epsilon) ? 0.0 : absb * Math.Sqrt(1.0 + Globals.SQR(absa / absb))));
        }
    }
}