wo80 / CSparse.NET

A concise library for solving sparse linear systems with direct methods.
GNU Lesser General Public License v2.1
58 stars 25 forks source link

is it possible to add UMFPACK solver ? #11

Closed bigworld12 closed 6 years ago

bigworld12 commented 6 years ago

although we already have 4 good solvers, i think we still need UMFPACK, which has the best performance for large square non-symmetric matrices, e.g. solving a 1 000 000 x 1 000 000 sparse matrix with ~900 000 non-zero elements using Cholesky with A+A` column ordering takes about 3 minutes and about 2 gb of ram for just the factorization. but with umfpack it takes only ~50 seconds

bigworld12 commented 6 years ago

here is more information about the solver https://people.sc.fsu.edu/~jburkardt/cpp_src/umfpack/umfpack.html also the source from SuiteSparse https://github.com/PetterS/SuiteSparse/tree/master/UMFPACK

wo80 commented 6 years ago

I already have bindings for some of the SuiteSparse solvers (Cholmod, UMFPACK, SPQR) and I still plan to publish them here on Github (already mentioned in https://github.com/wo80/CSparse.NET/issues/7#issuecomment-315426892).

It will be a separate project (CSparse.Interop). I don't have a lot of time to work on it, but if you're interested (and maybe want to contribute), I could publish the code I have so far.

bigworld12 commented 6 years ago

that would be awesome ! i am already working on implementing such algorithms into c# as the current available solutions are very complicated to build and were last updated 3~6 years ago, so i will gladly try my best to help you

wo80 commented 6 years ago

I've added a repository that contains a Visual Studio solution to compile SuiteSparse: https://github.com/wo80/vs-suitesparse

Let me know if this works for you (the only tricky part should be setting up BLAS/LAPACK). I hope I will find some time to publish the C# bindings next week.

wo80 commented 6 years ago

I've pushed some code to https://github.com/wo80/csparse-interop. Have a nice weekend!

A simple test program could look like

namespace ConsoleApp
{
    using CSparse.Double;
    using CSparse.Double.Factorization;
    using CSparse.Interop.Umfpack;
    using System;

    class TestUmfpack
    {
        public static void Test(SparseMatrix A)
        {
            int n = A.RowCount;

            var x = Vector.Create(n, 1.0);
            var b = Vector.Create(n, 0.0);

            // Save true solution.
            var s = Vector.Clone(x);

            A.Multiply(x, b);

            Vector.Clear(x);

            using (var solver = new Umfpack(A))
            {
                solver.Solve(b, x);

                var info = solver.Info;

                Console.WriteLine("ordering: {0}", (UmfpackOrdering)info.ORDERING_USED);
                Console.WriteLine("strategy: {0}", (UmfpackStrategy)info.STRATEGY_USED);
                Console.WriteLine("  LU nnz: {0}", info.LU_ENTRIES);
                Console.WriteLine("   rcond: {0}", info.RCOND);

                Console.WriteLine("L2-Error: {0:0.00e00}", ComputeL2Error(x, s));
            }
        }

        public static double ComputeL2Error(double[] actual, double[] expected, bool relativeError = true)
        {
            var e = Vector.Clone(actual);

            Vector.Axpy(-1.0, expected, e);

            if (relativeError)
            {
                return Vector.Norm(e) / Vector.Norm(expected);
            }

            return Vector.Norm(e);
        }
    }
}