Skip to content

Sparse Cholesky example

wo80 edited this page Jun 8, 2015 · 5 revisions

This example shows how to create a sparse Cholesky factorization:

using CSparse;
using CSparse.Double;
using CSparse.Double.Factorization;
using CSparse.IO;

public static class Example
{
    public static bool Solve(string filePath)
    {
        // Load matrix from a file.
        var A = MatrixMarketReader.ReadMatrix<double>(filePath);

        int m = A.RowCount;
        int n = A.ColumnCount;

        // Create test data.
        var x = Vector.Create(n, 1.0);
        var b = new double[m];
        var r = new double[m];

        // Compute right hand side vector b.
        A.Multiply(1.0, x, 0.0, b);

        // Apply column ordering to A to reduce fill-in.
        var order = ColumnOrdering.MinimumDegreeAtPlusA;

        if (m == n)
        {
            var chol = new SparseCholesky(A, order);

            // Solve Ax = b.
            Vector.Copy(b, x);
            chol.Solve(x);
        }
        else if (m > n)
        {
            // WARNING: hope you know what you're doing!
            var At = A.Transpose() as CompressedColumnStorage;
            var AtA = At.Multiply(A);

            var b2 = Vector.Create(b.Length, 0.0);
            var x2 = Vector.Create(b.Length, 0.0);

            At.Multiply(b, b2);

            var chol = new SparseCholesky(AtA, order);

            // Solving normal equation A'Ax = A'b
            chol.Solve(b2, x2);
            Vector.Copy(x2, x, n);
        }
        else
        {
            return false;
        }

        // Compute residual r = b - Ax.
        Vector.Copy(b, r);
        A.Multiply(-1.0, x, 1.0, r);

        return true;
    }
}