反幂法求矩阵特征值的C语言实现

#include <stdio.h>
#include <math.h>

int n = 5;

double * LUdecomposition(double A[n][n], double Y[n])
{
    static double result[4];//n
    double U[n][n];
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {           
            U[i][j] = 0;       
        } 
    }
    
    double L[n][n];
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (i < j)
            {
                L[i][j] = 0;
            }
            else if (i == j)
            {
                L[i][j] = 1;
            } 
            else
            {
                L[i][j] = 0;
            }    
        } 
    }

    for (int k = 0; k < n; k++)
    {
        for (int j = k; j < n; j++)
        {
            if (k == 0)
            {
                U[k][j] = A[k][j];
            }
            else
            {
                double sum = 0;
                for (int r = 0; r < k; r++)
                {
                    sum = sum +L[k][r]*U[r][j];
                }    
                U[k][j] = A[k][j] - sum;
            }
        }
        
        for (int i = k+1; i < n; i++)
        {
            if (k == 0)
            {
                L[i][k] = A[i][k]/U[k][k];
            }
            else
            {
                double sum = 0;
                for (int r = 0; r < k; r++)
                {
                    sum = sum + L[i][r]*U[r][k];
                }
                L[i][k] = (A[i][k] - sum)/U[k][k];
            }
        }
    }
    
    double Z[n];
    for (int i = 0; i < n; i++)
    {
        if (i == 0)
        {
            Z[i] = Y[i];
        }
        else
        {
            double sum = 0;
            for (int j = 0; j < i; j++)
            {
                sum = sum + L[i][j]*Z[j];
            }
            Z[i] = Y[i] - sum;
        }
    }
    
    for (int i = n-1; i >= 0; i--)
    {
        if (i == n-1)
        {
            result[i] = Z[i]/U[i][i];
        }
        else
        {
            double sum = 0;
            for (int j = i+1; j < n; j++)
            {
                sum = sum + U[i][j]*result[j];
            }
            
            result[i] = (Z[i] - sum)/U[i][i];
        }
    }
    
    return result;
}

double max(double X[n])
{
    double m = 0;
    for (int i = 0; i < n; i++)
    {
        if (fabs(X[i]) > m)
        {
            m = fabs(X[i]);
        }    
    }
    return m;
}

int max_order(double X[n])
{
    double m = 0;
    int mo = 0;
    for (int i = 0; i < n; i++)
    {
        if (fabs(X[i]) > m)
        {
            m = fabs(X[i]);
            mo = i;
        }    
    }
    return mo;
}

int main()
{   
    /*double A[5][5]= {{1.0/9, 1.0/8, 1.0/7, 1.0/6, 1.0/5}, 
                     {1.0/8, 1.0/7, 1.0/6, 1.0/5, 1.0/4}, 
                     {1.0/7, 1.0/6, 1.0/5, 1.0/4, 1.0/3}, 
                     {1.0/6, 1.0/5, 1.0/4, 1.0/3, 1.0/2}, 
                     {1.0/5, 1.0/4, 1.0/3, 1.0/2, 1.0/1}};*/

    double A[n][n];
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            A[i][j] = 1.0/(10 - i- j - 1);
        }  
    }
    
    double Y[n];
    for (int i = 0; i < n; i++)
    {
        Y[i] = 1;
    }
    
    printf(" k \t\t\t X_(k+1) \t\t\t\t\t\t\t\t Y_(k) \t\t\t\t 1/λ_(k) \n");
    printf("-------------------------------------------------------------------------------------------------------------------------------------------\n");

    double value = 1;
    int num = 0;
    double λ = 1; 
    double xi;
    
    while (value >= pow(10, -5))
    {
        double *X = LUdecomposition(A, Y);
        double m = max(X);
        int mo = max_order(X);
         
        value = fabs(1.0/λ - 1.0/(X[mo]/Y[mo]));
        λ = X[mo]/Y[mo]; 

        printf(" %d ", num);
        printf("\t(");
        for (int i = 0; i < n; i++)
        {
            printf("%.5f ", X[i]);
        }
        printf(")\t");
        printf("\t(");
        for (int i = 0; i < n; i++)
        {
            printf("%.5f ", Y[i]);
        }
        printf(")\t");
        printf("%.12f \n", 1.0/λ);
        
        for (int i = 0; i < n; i++)
        {
            Y[i] = X[i]/m;
        }

        if (num > 99)
        {
            break;
        }

        num = num +1;
    }
    printf("Modulo Minimum Eigenvalue = %.12f\n", 1.0/λ);
    printf("Modulo Minimum Eigenvector = ");
    printf("(");
    for (int i = 0; i < n; i++)
    {
        printf("%.5f ", Y[i]);
    }
    printf(")\n");

    return 0;
}