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

算法名称: 带状对角矩阵的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.01.00.00.00.00.00.0 },
            
...4.01.05.00.00.00.00.0 },
            
...9.02.06.05.00.00.00.0 },
            
...0.03.05.08.09.00.00.0 },
            
...0.00.07.09.03.02.00.0 },
            
...0.00.00.03.08.04.06.0 },
            
...0.00.00.00.02.04.04.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();
    }


}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值