我的线性代数库（LU分解求矩阵的逆矩阵，行列式的值）

#include "stdafx.h"

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <string.h>

// N可以随意指定
#define N 4

void main()
{
    double a[N][N];
    double L[N][N], U[N][N];
    double r[N][N], u[N][N];
    double out[N][N], out1[N][N], E[N][N];
    memset(a, 0, sizeof(a));    // 初始化
    memset(L, 0, sizeof(L));
    memset(U, 0, sizeof(U));
    memset(r, 0, sizeof(r));
    memset(u, 0, sizeof(u));
    int n = N;
    int k,i,j;
    double s,t;
    srand(time(0));

    // 输入矩阵
    printf("input A=\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            //scanf("%f", &a[i][j]);
            do{
                a[i][j] = rand() % 100;
            }while(a[i][j] == 0);
        }
    }
    
    printf("输入矩阵:\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf(" %lf", a[i][j]);
        }
        printf("\n");
    }

    // U矩阵的第一行不变
    //for(j=0;j<n;j++){
    //    a[0][j] = a[0][j];            // 计算U矩阵的第一行
    //}
    for(i=1;i<n;i++){
        a[i][0] = a[i][0]/a[0][0];    // 计算U矩阵的第一列
    }
    
    for(k=1;k<n;k++){
        for(j=k;j<n;j++){
            s = 0;
            for(i=0;i<k;i++){
                s += a[k][i] * a[i][j];    //累加
            }
            a[k][j] = a[k][j] - s;    // 计算U矩阵的其他元素
        }
        for(i=k+1;i<n;i++){
            t = 0;
            for(j=0;j<k;j++){
                t += a[i][j] * a[j][k];    //累加
            }
            a[i][k] = (a[i][k] - t) / a[k][k];    // 计算L矩阵的其他元素
        }
    }

    for(i=0; i<n; i++){
        for(j=0; j<n; j++){
            if(i > j){
                // 如果i>j，说明行大于列，计算矩阵的下三角部分，得出L的值，U的为0
                L[i][j] = a[i][j];
                U[i][j] = 0;
            }else if(i == j){
                U[i][j] = a[i][j];
                L[i][j] = 1;
            }else{
                // 否则 i<j 说明行小于列，计算矩阵的上三角部分，得出U值，L的为0
                U[i][j] = a[i][j];
                L[i][j] = 0;
            }
        }
    }

    double temp = 1.0;
    for(i=0;i<n;i++){
        temp *= U[i][i];
    }
    if(temp == 0){
        printf("\n矩阵不可逆，行列式值为0");
        return;
    }


    //求L和U矩阵的逆

    //求U矩阵的逆
    for(i=0;i<n;i++){
        u[i][i] = 1 / U[i][i];    // 对角元素的值直接取倒数
        for(k=i-1;k>=0;k--){
            s = 0;
            for(j=k+1;j<=i;j++){
                s += U[k][j] * u[j][i];
            }
            u[k][i] = -s / U[k][k];    // 迭代计算，按列倒序一次得到每一个值
        }
    }

    //求L矩阵的逆
    for(i=0;i<n;i++){
        r[i][i] = 1;    // 对角元素的值直接取倒数，这里为1
        for(k=i+1;k<n;k++){
            for(j=i;j<=k-1;j++){
                r[k][i] -= L[k][j]*r[j][i];    //迭代计算，按列顺序一次得到每一个值
            }
        }
    }


    // 打印L,U
    printf("\nL矩阵:");
    for(i=0;i<n;i++){
        printf("\n");
        for(j=0;j<n;j++){
            printf(" %lf", L[i][j]);
        }
    }
    printf("\nU矩阵:");
    for(i=0;i<n;i++){
        printf("\n");
        for(j=0;j<n;j++){
            printf(" %lf", U[i][j]);
        }
    }
    printf("\n");

    // 打印L,U的逆矩阵
    printf("\nL矩阵的逆矩阵:");
    for(i=0;i<n;i++){
        printf("\n");
        for(j=0;j<n;j++){
            printf(" %lf", r[i][j]);
        }
    }
    printf("\nU矩阵的逆矩阵:");
    for(i=0;i<n;i++){
        printf("\n");
        for(j=0;j<n;j++){
            printf(" %lf", u[i][j]);
        }
    }
    printf("\n");

    //验证L和U相乘,得到原矩阵
    printf("\nL矩阵和U矩阵乘积\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            out[i][j] = 0;
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            for(k=0;k<n;k++){
                out[i][j] += L[i][k]*U[k][j];
            }
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf(" %lf", out[i][j]);
        }
        printf("\n");
    }

    // 将r和u相乘,得到逆矩阵
    printf("\n原矩阵的逆矩阵\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            out1[i][j] = 0;
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            for(k=0;k<n;k++){
                out1[i][j] += u[i][k]*r[k][j];
            }
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf(" %lf", out1[i][j]);
        }
        printf("\n");
    }

    // 验证
    printf("\n原矩阵和其逆矩阵乘积\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            E[i][j] = 0;
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            for(k=0;k<n;k++){
                E[i][j] += out[i][k]*out1[k][j];
            }
        }
    }
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf(" %lf", E[i][j]);
        }
        printf("\n");
    }

    double det = 1;
    for(i=0;i<n;i++){
        det *= U[i][i];
    }
    printf("\nA的行列式值=%f\n",det);
}