算法名称：带状对角矩阵的LU分解及回代求解
算法描述：
分解主要是使用笔者前面几篇文章提到过的Crout方法。因为不可能把一个带状对角矩阵A的LU分解也像其压缩形式本是一样紧凑的存储起来，因为分解产生了附加的非零元素填入。一种直接的存储方案是，坝上三角因子（U）返回到以前占有的相同的空间中，把下三角因子（L）返回到单独的N×m1压缩矩阵中。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();
}
}
发表于 @ 2007年12月09日 23:04:00|评论(loading...)|编辑