LUP分解求解线性方程组及求逆矩阵 java

public class LU
{
/**矩阵相乘*/
public static float[][] multiply(float[][]A,float[][]B){
int n = A[0].length,i = 0,j = 0,k = 0;
float[][] C = new float

;
for(i = 0;i < n;i++){
for(j = 0;j < n;j++){
for(k = 0;k < n;k++){
C[i][j] += A[i][k]*B[k][j];
}
}
}
return C;
}

/**矩阵LU分解*/
public static void LU_Decomposition(float[][] A,float[][] L,float[][] U){
int n = A[0].length;
float u = 0f,l = 0f;
int i=0,j=0,k=0;
//initial L and U
for(i = 0;i < n-1;i++){
u = A[i][i];
//第j行
for(j = i+1;j < n;j++){
l = A[j][i]/u;
A[j][i] = l;
//update A'
for(k = i+1;k < n;k++){
A[j][k] = A[j][k] - l*A[i][k];
}
}
}//for
//make L and U
for(i=0;i < n;i++){
for(j = 0;j <= i;j++){
if(i!=j)
L[i][j] = A[i][j];
else L[i][j] = 1;
}
for(k = i;k < n;k++){
U[i][k] = A[i][k];
}
}
}
/**矩阵LUP分解
*	相比LU分解，LUP分解可以防止出现除以0的情况和除数非常小导致计算不准确的情况
*/
public static void LUP_Decomposition(float[][] A,float[][] L,float[][] U,int[] P){
int n = A[0].length;
int i = 0,j = 0,k = 0,row = 0;
for(i = 0;i < n;i++) P[i] = i;
for(i = 0;i < n-1;i++){
float p = 0f;
for(j = i;j < n;j++){
if(Math.abs(A[j][i]) > p){
p = Math.abs(A[j][i]);
row = j;
}
}
if(p == 0) {System.err.println("singular matrix");return;}
//exchange P[i] with P[row]
int tmp = P[i];
P[i] = P[row];
P[row] = tmp;

float tmp2 = 0f;
for(j = 0;j < n;j++){
//exchange A[i][j] with A[row][j]
tmp2 = A[i][j];
A[i][j] = A[row][j];
A[row][j] = tmp2;
}
//下面与LU分解一样了
float u = A[i][i],l = 0f;
for(j = i+1;j < n;j++){
l = A[j][i]/u;
A[j][i] = l;
//update A'
for(k = i+1;k < n;k++){
A[j][k] = A[j][k] - A[i][k]*l;
}
}
}//for
//make L and U
for(i=0;i < n;i++){
for(j = 0;j <= i;j++){
if(i!=j)
L[i][j] = A[i][j];
else L[i][j] = 1;
}
for(k = i;k < n;k++){
U[i][k] = A[i][k];
}
}
}

/**求解线性方程组
*	Ax = b ==> PAx = Pb ==>LUx = Pb ==> Ly = Pb ==> Ux = y
*/
public static float[] LUP_Solve(float[][] L,float[][] U,int[]P,float[] b){
int n = L[0].length,i = 0,j = 0;
float[] x = new float
;
float[] y = new float
;

//正向替换
for(i = 0;i < n;i++){
y[i] = b[P[i]];
for(j = 0;j < i;j++){
y[i] = y[i] - L[i][j]*y[j];
}
}
//反向替换
for(i = n-1;i >= 0; i--){
x[i] = y[i];
for(j = n-1;j > i;j--){
x[i] = x[i] - U[i][j]*x[j];
}
x[i] /= U[i][i];
}
return x;
}
/**矩阵置换*/
public static float[][] transposition(float[][] A){
int i=0,j=0,n = A[0].length;
float tmp = 0f;
for(i=0;i < n;i++){
for(j=0;j<i;j++){
tmp = A[i][j];
A[i][j] = A[j][i];
A[j][i] = tmp;
}
}
return A;
}
/**矩阵求逆
*	通过解线程方程组的方式求逆矩阵。b分别取单位阵的各个列向量
*	所得到的解向量x就是逆矩阵的各个列向量，拼成逆矩阵即可
*/
public static float[][] inverse(float[][] A){
int i=0,j=0,n = A[0].length;
float[][] inv = new float

;
float[][] B = new float

;
//init identity matrix
for(i = 0;i < n;i++){
B[i][i] = 1;
}
//each column of the inverse matrix
for(i = 0;i < n;i++){
float[][] L = new float

;
float[][] U = new float

;
int[] P = new int
;
float[][] copyOfA = copy(A);
LUP_Decomposition(copyOfA,L,U,P);
inv[i] = LUP_Solve(L,U,P,B[i]);
}
transposition(inv);
return inv;
}

public static void main(String[] args){

float[][] A = {{2,3,1,5},{6,13,5,19},{2,19,10,23},{4,10,11,31}};
float[][] L = new float[A.length][A.length];
float[][] U = new float[A.length][A.length];
/*
LU_Decomposition(A,L,U);
print(L);
System.out.println("=======================");
print(U);
System.out.println("**********************");

float[][] B = {{2,0,2,0.6f},{3,3,4,-2},{5,5,4,2},{-1,-2,3.4f,-1}};
L = new float[B.length][B.length];
U = new float[B.length][B.length];
int[] P = new int[A.length];
LUP_Decomposition(B,L,U,P);
print(L);
System.out.println("=======================");
print(U);
System.out.println("**********************");
*/
float[][] C = {{1,5,4},{2,0,3},{5,8,2}};
/*
L = new float[C.length][C.length];
U = new float[C.length][C.length];
int[] P = new int[C.length];
float[] b = {12,9,5};
float[] x;
LUP_Decomposition(C,L,U,P);
x = LUP_Solve(L,U,P,b);
System.out.println("解为:");
print(x);
*/
print(multiply(C,inverse(C)));
}

public static float[][] copy(float[][] A){
int len = A[0].length;
float[][] copy = new float[len][len];
for(int i = 0;i < len;i++){
for(int j=0;j < len; j++){
copy[i][j] = A[i][j];
}
}
return copy;
}
public static void print(float[][] a){
int len = a[0].length;
for(int i=0;i < len;i++){
for(int j=0;j < len;j++){
System.out.print(a[i][j]+"\t");
}
System.out.println();
}
}
public static void print(float[] a){
for(int i=0;i < a.length;i++){
System.out.print(a[i]+" ");
}
System.out.println();
}
}