LU Factorization of a Sparse Matrix

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

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 is used to solve the sparse linear systems and with iterative refinement. The reciprocal pivot growth factor and the reciprocal condition number are also computed.

using System;
using Imsl.Math;

public class SuperLUEx1
{
	public static void  Main(String[] args)
	{
		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.Equilibrate = false;
		slu.ColumnOrderingMethod = SuperLU.ColumnOrdering.Natural;
		slu.PivotGrowth = true;
		
		// Set option of iterative refinement
		slu.IterativeRefinement = true;
		
		// Solve sparse system A*x = b1
		Console.Out.WriteLine();
		Console.Out.WriteLine("Solve sparse System Ax=b1");
		Console.Out.WriteLine("=========================");
		Console.Out.WriteLine();
		
		sol = slu.Solve(b1);
		new PrintMatrix("Solution").Print(sol);
		
		Ferr = slu.ForwardErrorBound;
		Berr = slu.RelativeBackwardError;
		
		Console.Out.WriteLine();
		Console.Out.WriteLine("Forward error bound: " + Ferr);
		Console.Out.WriteLine("Relative backward error: " + Berr);
		Console.Out.WriteLine();
		Console.Out.WriteLine();
		
		
		// Solve sparse system (A^T)*x = b2
		Console.Out.WriteLine();
		Console.Out.WriteLine("Solve sparse System (A^T)*x=b2");
		Console.Out.WriteLine("==============================");
		Console.Out.WriteLine();
		
		sol = slu.SolveTranspose(b2);
		new PrintMatrix("Solution").Print(sol);
		
		Ferr = slu.ForwardErrorBound;
		Berr = slu.RelativeBackwardError;
		
		System.Console.Out.WriteLine();
		System.Console.Out.WriteLine("Forward error bound: " + Ferr);
		System.Console.Out.WriteLine("Relative backward error: " + Berr);
		System.Console.Out.WriteLine();
		System.Console.Out.WriteLine();
		
		// Compute reciprocal pivot growth factor and condition number
		
		recip_pivot_growth = slu.ReciprocalPivotGrowthFactor;
		conditionNumber = slu.ConditionNumber;
		
		Console.Out.WriteLine("Pivot growth factor and condition number");
		Console.Out.WriteLine("========================================");
		Console.Out.WriteLine();
		
		Console.Out.WriteLine("Reciprocal pivot growth factor: " +
         recip_pivot_growth);
		Console.Out.WriteLine("Reciprocal condition number: " + conditionNumber);
		Console.Out.WriteLine();
	}
}

Output


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

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


Forward error bound: 2.74489425897801E-14
Relative backward error: 0



Solve sparse System (A^T)*x=b2
==============================

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


Forward error bound: 2.41379101136191E-15
Relative backward error: 0


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

Reciprocal pivot growth factor: 1
Reciprocal condition number: 0.0244510978043912

Link to C# source.