C#，码海拾贝（35）

using System;

namespace Zhou.CSharp.Algorithm {     ///     /// 矩阵类     /// 作者：周长发     /// 改进：深度混淆     /// https://blog.csdn.net/beijinghorn     ///     public partial class Matrix     {

        ///         /// 求实对称矩阵特征值与特征向量的雅可比法         ///         /// 源矩阵         /// 一维数组，长度为矩阵的阶数，返回时存放特征值         /// 返回时存放特征向量矩阵，其中第i列为与数组dblEigenValue中第j个特征值对应的特征向量         /// 迭代次数         /// 计算精度         /// 求解是否成功         public static bool ComputeEvJacobi(Matrix src, out double[] dblEigenValue, out Matrix mtxEigenVector, int nMaxIt = 100, double eps = 1.0E-7)         {             int p = 0, q = 0, u, w, t, s;             double fm, cn, sn, omega, x, y, d;

            int n = src.Columns;             dblEigenValue = new double[n];             mtxEigenVector = new Matrix(n, n);

            int k = 1;             for (int i = 0; i < n; i++)             {                 mtxEigenVector[i * n + i] = 1.0;                 for (int j = 0; j < n; j++)                 {                     if (i != j)                     {                         mtxEigenVector[i * n + j] = 0.0;                     }                 }             }             while (true)             {                 fm = 0.0;                 for (int i = 1; i < n; i++)                 {                     for (int j = 0; j < i; j++)                     {                         d = Math.Abs(src[i * n + j]);                         if ((i != j) && (d > fm))                         {                             fm = d;                             p = i;                             q = j;                         }                     }                 }

                if (fm < eps)                 {                     for (int i = 0; i < n; ++i)                     {                         dblEigenValue[i] = src[i, i];                     }                     return true;                 }

                if (k > nMaxIt)                 {                     return false;                 }                 k = k + 1;                 u = p * n + q;                 w = p * n + p;                 t = q * n + p;                 s = q * n + q;                 x = -src[u];                 y = (src[s] - src[w]) / 2.0;                 if (Math.Abs(x) < float.Epsilon || Math.Abs(y) < float.Epsilon)                 {                     return false;                 }                 omega = x / Math.Sqrt(x * x + y * y);

                if (y < 0.0)                 {                     omega = -omega;                 }                 if (Math.Abs(omega - 1.0) < float.Epsilon)                 {                     return false;                 }                 sn = 1.0 + Math.Sqrt(1.0 - omega * omega);                 if (Math.Abs(sn) < float.Epsilon)                 {                     return false;                 }                 sn = omega / Math.Sqrt(2.0 * sn);                 if (Math.Abs(1.0 - sn) < float.Epsilon)                 {                     return false;                 }                 cn = Math.Sqrt(1.0 - sn * sn);                 fm = src[w];                 src[w] = fm * cn * cn + src[s] * sn * sn + src[u] * omega;                 src[s] = fm * sn * sn + src[s] * cn * cn - src[u] * omega;                 src[u] = 0.0;                 src[t] = 0.0;                 for (int j = 0; j < n; j++)                 {                     if ((j != p) && (j != q))                     {                         u = p * n + j; w = q * n + j;                         fm = src[u];                         src[u] = fm * cn + src[w] * sn;                         src[w] = -fm * sn + src[w] * cn;                     }                 }

                for (int i = 0; i < n; i++)                 {                     if ((i != p) && (i != q))                     {                         u = i * n + p;                         w = i * n + q;                         fm = src[u];                         src[u] = fm * cn + src[w] * sn;                         src[w] = -fm * sn + src[w] * cn;                     }                 }

                for (int i = 0; i < n; i++)                 {                     u = i * n + p;                     w = i * n + q;                     fm = mtxEigenVector[u];                     mtxEigenVector[u] = fm * cn + mtxEigenVector[w] * sn;                     mtxEigenVector[w] = -fm * sn + mtxEigenVector[w] * cn;                 }             }         }     } }