package com.imsl.test.example.math;

import com.imsl.math.*;

/**
 * <p>
 * SuperLU Example 1: Computes the LU factorization of a sparse matrix.</p>
 *
 *
 * <p>
 * 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} $$
 * </p>
 * <p>
 * is computed. The sparse coordinate form for \(A\) is given by a series of
 * row, column, and value triplets:
 * <center>
 * <table border="0">
 * <tr><th>row</th><th>column</th><th>value</th></tr>
 * <tr><td>\(0\)</td><td>\(0\)</td><td>\(10.0\)</td></tr>
 * <tr><td>\(1\)</td><td>\(1\)</td><td>\(10.0\)</td> </tr>
 * <tr><td>\(1\)</td><td>\(2\)</td><td>\(-3.0\)</td></tr>
 * <tr><td>\(1\)</td><td>\(3\)</td><td>\(-1.0\)</td></tr>
 * <tr><td>\(2\)</td><td>\(2\)</td><td>\(15.0\)</td></tr>
 * <tr><td>\(3\)</td><td>\(0\)</td><td>\(-2.0\)</td></tr>
 * <tr><td>\(3\)</td><td>\(3\)</td><td>\(10.0\)</td></tr>
 * <tr><td>\(3\)</td><td>\(4\)</td><td>\(-1.0\)</td></tr>
 * <tr><td>\(4\)</td><td>\(0\)</td><td>\(-1.0\)</td></tr>
 * <tr><td>\(4\)</td><td>\(3\)</td><td>\(-5.0\)</td></tr>
 * <tr><td>\(4\)</td><td>\(4\)</td><td>\(1.0\)</td></tr>
 * <tr><td>\(4\)</td><td>\(5\)</td><td>\(-3.0\)</td></tr>
 * <tr><td>\(5\)</td><td>\(0\)</td><td>\(-1.0\)</td></tr>
 * <tr><td>\(5\)</td><td>\(1\)</td><td>\(-2.0\)</td></tr>
 * <tr><td>\(5\)</td><td>\(5\)</td><td>\(6.0\)</td></tr>
 * </table>
 * </center>
 * </p>
 * <p>
 * Let \(y^T = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)\). Then we have \(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 \).
 * </p>
 * <p>
 * In the example, the The LU factorization of \(A\) is used to solve the sparse
 * linear system $$\begin{array}{cc} Ax &=& b_1 \\ A^Tx &=& b_2 \end{array}$$
 * <p>
 * with iterative refinement. The reciprocal pivot growth factor and the
 * reciprocal condition number are also computed. Note that by construction
 * \(x=y\) is the solution to this system.
 * </p>
 *
 * @see <a href="SuperLUEx1.java">Code</a>
 * @see <a href="doc-files/SuperLUEx1.out">Output</a>
 */
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();
    }
}