Example: LU Factorization of a Sparse Matrix

The LU Factorization of the sparse 6 \times 6 matrix

 A=\begin{pmatrix} 10.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 10.0 & -3.0 & -1.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 15.0 & 0.0 & 0.0 & 0.0 \\ -2.0 & 0.0 & 0.0 & 10.0 & -1.0 & 0.0 \\ -1.0 & 0.0 & 0.0 & -5.0 & 1.0 & -3.0 \\ -1.0 & -2.0 & 0.0 & 0.0 & 0.0 & 6.0 \end{pmatrix}
is computed. The sparse coordinate form for A is given by a series of row, column, value triplets:

row column value
0 0 10.0
1 1 10.0
1 2 -3.0
1 3 -1.0
2 2 15.0
3 0 -2.0
3 3 10.0
3 4 -1.0
4 0 -1.0
4 3 -5.0
4 4 1.0
4 5 -3.0
5 0 -1.0
5 1 -2.0
5 5 6.0

Let

y^T = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
so that b_1:=Ay = {(10.0, 7.0, 45.0, 33.0, -34.0, 31.0)}^T, and
b_2:=A^Ty = {(-9.0, 8.0, 39.0, 13.0, 1.0, 21.0)}^T
.

The LU factorization of A is used to solve the sparse linear systems Ax=b_1 and A^Tx=b_2 with iterative refinement. The reciprocal pivot growth factor and the reciprocal condition number are also computed.


import com.imsl.math.*;

public class SuperLUEx1 {

    public static void main(String args[]) throws Exception {
        int m;
        SuperLU slu;
        double conditionNumber, recip_pivot_growth;
        double[] sol = null;
        double Ferr, Berr;

        double[] b1 = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
        double[] b2 = {-9.0, 8.0, 39.0, 13.0, 1.0, 21.0};

        // Initialize matrix A.
        m = 6;
        SparseMatrix a = new SparseMatrix(m, m);

        a.set(0, 0, 10.0);
        a.set(1, 1, 10.0);
        a.set(1, 2, -3.0);
        a.set(1, 3, -1.0);
        a.set(2, 2, 15.0);
        a.set(3, 0, -2.0);
        a.set(3, 3, 10.0);
        a.set(3, 4, -1.0);
        a.set(4, 0, -1.0);
        a.set(4, 3, -5.0);
        a.set(4, 4, 1.0);
        a.set(4, 5, -3.0);
        a.set(5, 0, -1.0);
        a.set(5, 1, -2.0);
        a.set(5, 5, 6.0);

        // Compute the sparse LU factorization of a
        slu = new SuperLU(a);

        slu.setEquilibrate(false);
        slu.setColumnPermutationMethod(SuperLU.NATURAL_ORDERING);
        slu.setPivotGrowth(true);

        //slu.setPerformanceTuningParameters(1,2);
        //slu.setPerformanceTuningParameters(2,1);
        // Set option of iterative refinement
        slu.setIterativeRefinement(true);

        // Solve sparse system A*x = b1
        System.out.println();
        System.out.println("Solve sparse System Ax=b1");
        System.out.println("=========================");
        System.out.println();

        sol = slu.solve(b1);
        new PrintMatrix("Solution").print(sol);

        Ferr = slu.getForwardErrorBound();
        Berr = slu.getRelativeBackwardError();

        System.out.println();
        System.out.println("Forward error bound: " + Ferr);
        System.out.println("Relative backward error: " + Berr);
        System.out.println();
        System.out.println();

        // Solve sparse system A^T*x = b2
        System.out.println();
        System.out.println("Solve sparse System A^Tx=b2");
        System.out.println("===========================");
        System.out.println();

        sol = slu.solveTranspose(b2);
        new PrintMatrix("Solution").print(sol);

        Ferr = slu.getForwardErrorBound();
        Berr = slu.getRelativeBackwardError();

        System.out.println();
        System.out.println("Forward error bound: " + Ferr);
        System.out.println("Relative backward error: " + Berr);
        System.out.println();
        System.out.println();

        // Compute reciprocal pivot growth factor and condition number
        recip_pivot_growth = slu.getReciprocalPivotGrowthFactor();
        conditionNumber = slu.getConditionNumber();

        System.out.println("Pivot growth factor and condition number");
        System.out.println("========================================");
        System.out.println();

        System.out.println("Reciprocal pivot growth factor: "
                + recip_pivot_growth);
        System.out.println("Reciprocal condition number: " + conditionNumber);
        System.out.println();
    }
}

Output


Solve sparse System Ax=b1
=========================

Solution
   0  
0  1  
1  2  
2  3  
3  4  
4  5  
5  6  


Forward error bound: 2.7448942589780064E-14
Relative backward error: 0.0



Solve sparse System A^Tx=b2
===========================

Solution
   0  
0  1  
1  2  
2  3  
3  4  
4  5  
5  6  


Forward error bound: 2.413791011361905E-15
Relative backward error: 0.0


Pivot growth factor and condition number
========================================

Reciprocal pivot growth factor: 1.0
Reciprocal condition number: 0.02445109780439122

Link to Java source.