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

算法名称： 带状对角矩阵的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();
    }

}