工程计算-DoolittleのLU分解

import java.util.Scanner;

/**
 * @author Janus
 * @creator 2020-02-21-18:43
 */
public class Main {
    public static void main(String[] Args) {
        Scanner scanner=new Scanner(System.in);
        int r=scanner.nextInt();
        int c=scanner.nextInt();
        double[][] a = new double[r][c];
        for(int i=0;i<r;++i){
            for(int j=0;j<c;++j){
                a[i][j]=scanner.nextDouble();
            }
        }
        Output(a);
        double[] b=new double[r];
        for(int i=0;i<r;i++){
            b[i]=scanner.nextDouble();
        }
        double[][] l = LUDecomposition(a, "l");
        double[][] u = LUDecomposition(a, "u");
        double[] y=SolveY(l,b);
        for(double temp:y){
            System.out.printf("%+3.6f\t", temp);
        }
        System.out.print('\n');
        double[] x=SolveX(u,y);
        for(double temp:x){
            System.out.printf("%+3.6f\t", temp);
        }
        System.out.print('\n');
    }

    public static double[] SolveX(double[][] u,double[] y){
        double[]x=new double[y.length];
        double temp;
        for(int i=y.length-1;i>=0;--i){
            temp=y[i];
            for(int j=y.length-1;j>i;--j){
                temp-=x[j]*u[i][j];
            }
            x[i]=temp/u[i][i];
        }
        return x;
    }

    public static double[] SolveY(double[][] l,double[] b){
        double[] y=new double[b.length];
        double temp;
        for(int i=0;i<b.length;++i){
            temp=b[i];
            for(int j=0;j<i;++j){
                temp-=l[i][j]*y[j];
            }
            y[i]=temp/l[i][i];
        }
        return y;
    }

    public static double[][] LUDecomposition(double a[][], String X) {
        int row = a.length;//row为矩阵行数
        int line = a[0].length;//line为矩阵列数
        double temp = 0;
        if (row != line) {
            System.out.println("矩阵行列数不等！！");
            return a;
        }//如果矩阵行列数不等，返回原矩阵。
        double[][] l = new double[row][line];//定义L矩阵
        double[][] u = new double[row][line];//定义U矩阵
        int k = 0;
        for (int j = 0; j < line; j++) {//首先计算U矩阵第一行各元素
            u[k][j] = a[k][j];
        }
        for (int i = 0; i < line; i++) {//计算L矩阵第一列各元素
            l[i][k] = a[i][k] / u[k][k];
        }
        for (k = 1; k < row; k++){//计算U矩阵第k行元素
            for (int j = k; j < row; j++) {
                temp = 0;
                for (int p = 0; p < k; p++) {
                    temp += l[k][p] * u[p][j];
                }
                u[k][j] = a[k][j] - temp;
            }//求出第k行的u
            for (int i = k; i < line; i++){ //计算L矩阵第k列元素
                temp = 0;
                for (int p = 0; p < k; p++) {
                    temp += l[i][p] * u[p][k];
                }
                l[i][k] = (a[i][k] - temp) / u[k][k];
            }//求出第k列的j
        }
        System.out.println("L矩阵为：");
        Output(l);//打印矩阵l,此为自定义函数
        System.out.println("U矩阵为：");
        Output(u);//打印矩阵l,此为自定义函数
        if (X == "U" || X == "u")
            return u;//返回U矩阵
        else if (X == "L" || X == "l")
            return l;//返回L矩阵
        else
            return a;//返回原矩阵
    }

    public static final void Output(double[][] matrix) {
        int rowNum=matrix.length;
        int columnNum=matrix[0].length;
        for (int i = 0; i < rowNum; ++i) {
            for (int j = 0; j < columnNum; ++j) {
                System.out.printf("%+3.6f\t", matrix[i][j]);
            }
            System.out.print("\n");
        }
        System.out.print("\n");
    }


}