c_cpp 计算存储在一维点阵列中的矩阵的特征值和特征向量

Posted 2021-05-14

/// compute eigenvalues and eigen vectors of matrix stored in 1-D point array using Jacobi method.
/// </summary>
/// <param name="a">symmetric matrix with n*n size, and also used to return eigenvalues in diagonal place.</param>
/// <param name="n">size of matrix.</param>
/// <param name="v">eigen vector matrix, stored in column format.</param>
/// <param name="eps">precision parameter.</param>
/// <param name="jt">max iteration times.</param>
/// <returns>false: can't satisfy precision, and hit the max iteration times
///          true: success.</returns>
bool PcaFusion::eejcb( double a[], int n, double v[], double eps, int jt )
{
    int i, j, p, q, u, w, t, s, l;
    double fm, cn, sn, omega, x, y, d;

    l = 1;
    //初始化特征向量矩阵使其全为0
    for( i = 0; i <= n - 1; i++ )
    {
        v[i * n + i] = 1.0;
        for( j = 0; j <= n - 1; j++ )
        {
            if( i != j )
            {
                v[i * n + j] = 0.0;
            }
        }
    }

    while( true ) //循环
    {
        fm = 0.0;
        for( i = 0; i <= n - 1; i++ )   // 出, 矩阵a( 特征值 ), 中除对角线外其他元素的最大绝对值
        {
            //这个最大值是位于a[p][q] ,等于fm
            for( j = 0; j <= n - 1; j++ )
            {
                d = fabs( a[i * n + j] );

                if( ( i != j ) && ( d > fm ) )
                {
                    fm = d;
                    p = i;
                    q = j;
                }
            }
        }

        if( fm < eps )   //精度复合要求
        {
            return true;    //正常返回
        }

        if( l > jt )     //迭代次数太多
        {
            return false;    //失败返回
        }

        l ++;       //   迭代计数器
        u = p * n + q;
        w = p * n + p;
        t = q * n + p;
        s = q * n + q;
        x = -a[u];
        y = ( a[s] - a[w] ) / 2.0;		//x y的求法不同
        omega = x / sqrt( x * x + y * y );	//sin2θ

        //tan2θ=x/y = -2.0*a[u]/(a[s]-a[w])
        if( y < 0.0 )
        {
            omega = -omega;
        }

        sn = 1.0 + sqrt( 1.0 - omega * omega );
        sn = omega / sqrt( 2.0 * sn );		//sinθ
        cn = sqrt( 1.0 - sn * sn );			//cosθ

        fm = a[w];   //   变换前的a[w]   a[p][p]
        a[w] = fm * cn * cn + a[s] * sn * sn + a[u] * omega;
        a[s] = fm * sn * sn + a[s] * cn * cn - a[u] * omega;
        a[u] = 0.0;
        a[t] = 0.0;

        //   以下是旋转矩阵,旋转了了p行,q行,p列,q列
        //   但是四个特殊点没有旋转(这四个点在上述语句中发生了变化)
        //   其他不在这些行和列的点也没变
        //   旋转矩阵,旋转p行和q行
        for( j = 0; j <= n - 1; j++ )
        {
            if( ( j != p ) && ( j != q ) )
            {
                u = p * n + j;
                w = q * n + j;
                fm = a[u];
                a[u] = a[w] * sn + fm * cn;
                a[w] = a[w] * cn - fm * sn;
            }
        }

        //旋转矩阵,旋转p列和q列
        for( i = 0; i <= n - 1; i++ )
        {
            if( ( i != p ) && ( i != q ) )
            {
                u = i * n + p;
                w = i * n + q;
                fm = a[u];
                a[u] = a[w] * sn + fm * cn;
                a[w] = a[w] * cn - fm * sn;
            }
        }

        //记录旋转矩阵特征向量
        for( i = 0; i <= n - 1; i++ )
        {
            u = i * n + p;
            w = i * n + q;
            fm = v[u];
            v[u] = v[w] * sn + fm * cn;
            v[w] = v[w] * cn - fm * sn;
        }
    }

    return true;
}