矩阵的LU分解和求逆

package Algorithms_算法.Class._Matrix;

import java.util.Arrays;

public class LU {
    double[][] matrix;
    double[][] expand_Matrix;
    double[][] reverse_matrix;
    double[][] l;
    double[][] u;
    public void initialize(double[][] matrix){
        if( matrix[0] == null || matrix.length != matrix[0].length){
            System.out.println("不是方阵或者矩阵是空的");
        }
        this.matrix = matrix;
        l = new double[matrix.length][matrix.length];
        u = new double[matrix.length][matrix.length];
        expand_Matrix = new double[matrix.length][matrix.length * 2];
        reverse_matrix = new double[matrix.length][matrix.length];
    }
    public void solve(){
        int n = matrix.length;
//        L = {*            U = {****
//             **                 ***
//             ***                 **
//             ****}                *}
//        因为A=LU，所以A的第一行等于L的第一行第一个元素乘上U的每一列，
//         A的第一列等于L的每一行乘上U的第一行第一个元素
//        U[1][i]=A[1][i],L[i][1]=A[i][1]/U[1][1]
        for(int i=0;i<n;i++){
            u[0][i] = matrix[0][i];//upper的第一行和原矩阵相同，因为原矩阵等于LU相乘，而L的第一行只有一个元素
            //                                                                                              u【0】【0】不就等于matrix【0】【0】的平方了吗
            //确实，此外的u的每一列第一个元素等于原来位置的元素
            l[i][0] = matrix[i][0] / u[0][0];
            //L的第一个元素就等于源数据的倒数
        }

        //计算U的第r行，L的第r列元素
//        a[r][i] = Σ(k=0->n) l[r][k] * u[k][i],这是最基本的规律
//        但是当k大于r的时候，后面的乘积恒为0(LU矩阵性质),所以k <= r,i
//        a[r][i] = Σ(k=0->r) l[r][k] * u[k][i]
//        u[r][i] = a[r][i] - Σ(k=1 -> r-1) l[r][k] * u[k][i]    (i=r,r+1,...,n)
//        l[i][r] = a[i][r] - Σ(k=1 -> r-1) l[i][k] * u[k][r]    (i=r+1,r+2,...,n且r!=n)
        for (int r = 1; r < n; r++) {
            for (int i = r; i < n; i++) {
                u[r][i] = matrix[r][i] - sumLrkUki(r, i);
//                                                                                      if(i==r) l[r][r]=1;
//                                                                                      else if(r==n) l[n][n]=1;
                l[i][r] = (matrix[i][r] - sumLikUkr(r, i)) / u[r][r];
            }
        }
    }
    private double sumLrkUki(int r, int i) {
        double re = 0.0;
        for (int k = 0; k < r; k++) {
            re += l[r][k] * u[k][i];
        }
        return re;
    }
    private double sumLikUkr(int r, int i) {
        double re = 0.0;
        for (int k = 0; k < r; k++) {
            re += l[i][k] * u[k][r];
        }
        return re;
    }
    public boolean reverse(){
        int i,j,k=0;
        int num = matrix.length;
        double[] temp = new double[num*2];
        double coe=1;
        for(i=0;i<num;i++) {
            for(j=0;j<num;j++)
                expand_Matrix[i][j]=matrix[i][j];
        }//扩充矩阵初始化
        for(i=0;i<num;i++) {
            for(j=num;j<2*num;j++){
                expand_Matrix[i][j] = 0;
                expand_Matrix[i][i+num] = 1;
            }
        }
        //放入单位矩阵


        for(i=0;i<num;i++){
            //对矩阵进行重排，确保矩阵主对角线上元素不为零
            if(expand_Matrix[i][i]==0){
                if(i==num-1)
                    return false;
                //判断矩阵是否满秩
                for(j=i;j<num;j++){
                    if(expand_Matrix[j][i] != 0){
                        k = j;
                        break;
                    }
                }
                temp = expand_Matrix[k];
                expand_Matrix[k] = expand_Matrix[i];
                expand_Matrix[i] = temp;
            }
            //初等行变换
            for(j=0;j<num;j++){
                if(j!=i){
                    if(expand_Matrix[j][i]!= 0){
                        coe = expand_Matrix[j][i]/expand_Matrix[i][i];
                        for(k=i;k<2*num;k++){
                            expand_Matrix[j][k] -= coe*expand_Matrix[i][k];}
                    }
                }
            }
            //对主对角线元素对应的每一列，用该数除以主对角线的数得到系数，
            //用该行每一个数减去主对角线所在的行对应列的数乘以该系数，
            //确保该列每一个数都是零
            coe = expand_Matrix[i][i];
            for(j=i;j<2*num;j++)
                if(coe!=0)
                    expand_Matrix[i][j]/=coe;
        }//对主对角线的数，用该行每一个数（包括自身）除以对角线上的数，确保主对角线上的数为一，得到左侧单位阵

        for(i=0;i<num;i++){
            for(j=0;j<num;j++){
                reverse_matrix[i][j] = expand_Matrix[i][j+num];
            }
        }
        return true;
    }
    public static void main(String[] args) {
        double[][] data= {
                {2,2,6},
                {2,5,15},
                {6,15,46},
        };
        LU demo = new LU();
        demo.initialize(data);
        double[][] l = demo.l;
        double[][] u = demo.u;
        double[][] reverse_matrix = demo.reverse_matrix;
        demo.solve();
        System.out.println("U阵：");
        for (int i = 0; i < l.length; i++) {
            System.out.println(Arrays.toString(u[i]));
        }
        System.out.println("---------------------");
        System.out.println("L阵：");
        for (int i = 0; i < u.length; i++) {
            System.out.println(Arrays.toString(l[i]));
        }

        demo.reverse();
        System.out.println("---------------------");
        System.out.println("逆矩阵：");
        for (int i = 0; i < reverse_matrix.length; i++) {
            System.out.println(Arrays.toString(reverse_matrix[i]));
        }
    }
    //U的每一行就是这一轮，轮数所在的列被清空之后，留下的元素，相当于保存了每一次处等行变换之后的A矩阵
    //L的每一列就是这一轮，从这一轮所在的行开始向下，每一列与当前基准列的，要清掉的元素的比值，相当于保存了每一次的变换矩阵
    //因为自己是自己的1倍，所以L主对角线对应每一次作比值的时候和自己的比值，一定是1
}