带状对角矩阵的LU分解及回代求解算法实现

来源：互联网发布：淘宝金冠店铺排名编辑：程序博客网时间：2024/04/28 00:27

算法名称：带状对角矩阵的LU分解及回代求解

算法描述：

分解主要是使用笔者前面几篇文章提到过的Crout方法。因为不可能把一个带状对角矩阵A的LU分解也像其压缩形式本是一样紧凑的存储起来，因为分解产生了附加的非零元素填入。一种直接的存储方案是，坝上三角因子（U）返回到以前占有的相同的空间中，把下三角因子（L）返回到单独的N×m₁压缩矩阵中。U的对角线元素被存放在A的存储空间的第一列。

一旦矩阵A被分解，任意数目的右端项都可以通过重复调用回代函数求解，方法参考笔者前面的几篇文章。

运行示例：

The origin matrix:

3.0 1.0 0.0 0.0 0.0 0.0 0.0

4.0 1.0 5.0 0.0 0.0 0.0 0.0

9.0 2.0 6.0 5.0 0.0 0.0 0.0

0.0 3.0 5.0 8.0 9.0 0.0 0.0

0.0 0.0 7.0 9.0 3.0 2.0 0.0

0.0 0.0 0.0 3.0 8.0 4.0 6.0

0.0 0.0 0.0 0.0 2.0 4.0 4.0

----------------------------------

The matrix (after compression):

0.0 0.0 3.0 1.0

0.0 4.0 1.0 5.0

9.0 2.0 6.0 5.0

3.0 5.0 8.0 9.0

7.0 9.0 3.0 2.0

3.0 8.0 4.0 6.0

2.0 4.0 4.0 0.0

----------------------------------

The U matrix compressed:

9.0 2.0 6.0 5.0

3.0 5.0 8.0 9.0

7.0 9.0 3.0 2.0

-5.280423280423281 -1.2539682539682542 -0.6137566137566138 0.0

7.287575150300601 3.6513026052104207 6.0 0.0

2.9979375773408496 2.353361748934415 0.0 0.0

-0.4729407448174646 0.0 0.0 0.0

----------------------------------

The L matrix compressed:

0.4444444444444444 0.3333333333333333

0.11111111111111112 0.037037037037037056

0.3068783068783069 -0.36507936507936506

-0.13827655310621242 -0.5681362725450901

-0.010724597827581485 0.27443970851093086

0.2283067327095946 0.0

0.0 0.0

----------------------------------

The right vector:

1.0

2.0

3.0

4.0

5.0

6.0

7.0

----------------------------------

The solution vector:

-0.12296353762606659

1.3688906128782

0.22459270752521324

0.00426687354538399

-0.14041892940263748

1.905352986811482

-0.08514352211016335

----------------------------------

示例程序：

package com.nc4nr.chapter02.bandec;

public class BanDec ...{

double[][] a = ...{

...{ 3.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },

...{ 4.0, 1.0, 5.0, 0.0, 0.0, 0.0, 0.0 },

...{ 9.0, 2.0, 6.0, 5.0, 0.0, 0.0, 0.0 },

...{ 0.0, 3.0, 5.0, 8.0, 9.0, 0.0, 0.0 },

...{ 0.0, 0.0, 7.0, 9.0, 3.0, 2.0, 0.0 },

...{ 0.0, 0.0, 0.0, 3.0, 8.0, 4.0, 6.0 },

...{ 0.0, 0.0, 0.0, 0.0, 2.0, 4.0, 4.0 }

};

int anrow = 7,

m1 = 2, // 对角线下有m1行

m2 = 1; // 对角线上有m2行

double[][] ca = new double[anrow][m1+m2+1];

double[][] al = new double[anrow][m1];

double[][] b = ...{

...{1.0},

...{2.0},

...{3.0},

...{4.0},

...{5.0},

...{6.0},

...{7.0}

};

int[] indx = new int[anrow];

double parity = 1.0;

private void cmpban() ...{

int n = anrow;

for (int i = 0; i < n; i++) ...{

int k = i - m1;

int tmploop = Math.min(m1+m2+1,n - k);

for (int j = Math.max(0, -k); j < tmploop; j++) ca[i][j] = a[i][j+k];

}

private void bandec() ...{

int mm = m1 + m2 + 1,

l = m1,

n = anrow;

for (int i = 0; i < m1; i++) ...{

for (int j=m1-i; j<mm; j++) ca[i][j-l] = ca[i][j];

l--;

for (int j=mm-l-1; j<mm; j++) ca[i][j] = 0.0;

}

l = m1;

for (int k = 0; k < n; k++) ...{

double dum = ca[k][0];

int i = k;

if (l<n) l++;

for (int j = k+1; j < l; j++) ...{

if (Math.abs(ca[j][0]) > Math.abs(dum)) ...{

dum = ca[j][0];

i = j;

}

indx[k] = i+1;

if (i != k) ...{

parity = -parity;

for (int j=0; j<mm; j++) ...{

double mid = ca[k][j];

ca[k][j] = ca[i][j];

ca[i][j] = mid;

}

for (int j = k+1; j < l; j++) ...{

dum = ca[j][0]/ca[k][0];

//System.out.println("dum=" + dum);

al[k][j-k-1]=dum;

for (int ji = 1; ji<mm; ji++) ...{

// System.out.println("ca[" + k + "][" + ji + "]=" + ca[k][ji]);

// System.out.println("ca[" + j + "][" + ji + "]=" + ca[j][ji]);

ca[j][ji-1]=ca[j][ji]-dum*ca[k][ji];

// System.out.println("ca[" + j + "][" + (ji-1) + "]=" + ca[j][ji-1]);

// System.out.println("--------------------------------------");

}

ca[j][mm-1] = 0.0;

}

private void banbks() ...{

int mm = m1 + m2 + 1,

n = anrow,

l = m1;

for (int k = 0; k < n; k++) ...{

int j = indx[k] - 1;

if (j != k) ...{

double mid = b[j][0];

b[j][0] = b[k][0];

b[k][0] = mid;

}

if (l < n) l++;

for (int i = k+1; i < l; i++) b[i][0] -= al[k][i-k-1]*b[k][0];

}

l = 1;

for (int i = n-1; i >= 0; i--) ...{

double dum = b[i][0];

for (int k = 1; k < l; k++) dum -= ca[i][k]*b[k+i][0];

b[i][0] = dum/ca[i][0];

if (l < mm) l++;

}

private void output(double[][] a, int nrow, int ncol) ...{

for (int i = 0; i < nrow; i++) ...{

String str = "";

for (int j = 0; j < ncol; j++) ...{

str += a[i][j] + " ";

}

System.out.println(str);

}

System.out.println("----------------------------------");

}

public BanDec() ...{

System.out.println("The origin matrix:");

output(a,anrow,anrow);

cmpban();

System.out.println("The matrix (after compression):");

output(ca,anrow,m1+m2+1);

bandec();

System.out.println("The U matrix compressed:");

output(ca,anrow,m1+m2+1);

System.out.println("The L matrix compressed:");

output(al,anrow,m1);

System.out.println("The right vector:");

output(b,anrow,1);

banbks();

System.out.println("The solution vector:");

output(b,anrow,1);

}

public static void main(String[] args) ...{

new BanDec();

}