对称正定矩阵LDL^T分解的C语言实现

作业不能白做

于 2023-05-29 12:58:42 发布

阅读量294

点赞数

文章标签：矩阵 c语言算法

本文链接：https://blog.csdn.net/weixin_53383166/article/details/130925991

版权

该代码实现展示了如何使用C++对对称正定矩阵进行LDLT分解，其中L是单位下三角矩阵，D是对角矩阵，LT是L的转置。程序首先初始化单位矩阵，然后计算D和L，最后解决线性方程组。

摘要由CSDN通过智能技术生成

//当A为对称正定矩阵时，有分解 ：A = LDLT，其中L为单位下三角阵，D为对角阵，LT为L的转置

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

void print_matrix(int row, int col, double matrix[row][col])
{
    for(int i = 0;i < row; i++)
    {
        for(int j = 0;j < col; j++)
        {
            printf("%lf\t",matrix[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
}

void print_vector(int n, double vector[n])
{
    for(int i = 0;i < n;i ++)
    {
        printf("%.12lf\n",vector[i]);   
    }
    printf("----------------------------------------------");
    printf("\n");
}

int main()
{
    int n = 3;
    double A[3][3] = {{1, -1, 1}, {-1, 3, -2}, {1, -2, 4.5}};//A[n][n]
    printf("A = \n");
    print_matrix(n, n, A);

    double L[n][n];
    for (int j = 0; j < n; j++)//initialize L = I
    {
        for (int k = 0; k < n; k++)
        {
            if (k == j)
            {
                L[j][k] = 1;
            }
            else
            {
                L[j][k] = 0;
            }
        }        
    }
    double D[n][n];
    for (int j = 0; j < n; j++)//initialize D = I
    {
        for (int k = 0; k < n; k++)
        {
            if (k == j)
            {
                D[j][k] = 1;
            }
            else
            {
                D[j][k] = 0;
            }
        }        
    }
    //------------------------------------------------------------------------------------------------------------------
    for (int k = 0; k < n; k++)//calculate D, L
    {
        if (k ==0)
        {
            D[k][k] = A[k][k];
        }
        else
        {
            double sum1 = 0;
            for (int r = 0; r <= k-1; r++)
            {
                sum1 = sum1 + D[r][r]*L[k][r]*L[k][r];
            }
            D[k][k] = A[k][k] - sum1;
        }
        for (int i = k+1; i < n; i++)
        {
            if (k ==0)
            {
                L[i][k] = A[i][k]/D[k][k];
            }
            else
            {
                double sum2 = 0;
                for (int r = 0; r <= k-1; r++)
                {
                    sum2 = sum2 + D[r][r]*L[i][r]*L[k][r];
                }
                L[i][k] = (A[i][k] - sum2)/D[k][k];
            }
        }        
    }
    printf("L = \n");
    print_matrix(n, n, L);
    printf("D = \n");
    print_matrix(n, n, D);
    //------------------------------------------------------------------------------------------------------------------
    //solve LZ = b
    double b[3] = {4, 6, 5};
    double Z[n];
    for (int i = 0; i < n; i++)//calculate Z
    {
        if (i == 0)
        {
            Z[i] = b[i];
        }
        else
        {
            double sum2 = 0;
            for (int j = 0; j <= i-1; j++)
            {
                sum2 = sum2 + L[i][j]*Z[j];
            }
            Z[i] = b[i] - sum2;
        }        
    }
    //solve DY = Z
    double Y[n];
    for (int i = 0; i < n; i++)//calculate Y
    {
        Y[i] = Z[i]/D[i][i];
    }
    //solve LTX = Y
    double X[n];
    for (int i = n-1; i >= 0; i--)//calculate X
    {
        if (i == n-1)
        {
            X[i] = Y[i];
        }
        else
        {
            double sum3 = 0;
            for (int j = i+1; j <n; j++)
            {
                sum3 = sum3 + L[j][i]*X[j];
            }
            X[i] = Y[i] - sum3;
        }        
    }
    //print Z, Y, X
    printf("Z = \n");
    print_vector(n, Z);
    printf("Y = \n");
    print_vector(n, Y);
    printf("X = \n");
    print_vector(n, X);

    return 0;
}