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];

        // 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 = SparseCholesky.Create(A, order);

            // Solve Ax = b (overwrite x).
            chol.Solve(b, x);
        }
        else if (m > n)
        {
            var At = A.Transpose();
            var AtA = At.Multiply(A);

            var c = Vector.Create(n, 0.0);

            At.Multiply(b, c);

            var chol = SparseCholesky.Create(AtA, order);

            // Solve normal equation A'Ax = A'b (overwrite x).
            chol.Solve(c, x);
        }
        else
        {
            return false;
        }

        // Compute residual b - Ax (overwrite b).
        A.Multiply(-1.0, x, 1.0, b);

        return true;
    }
}