Problem with UMFPACK and large matrices

[vstembera]

Hi all,

I use Sparsuite library via the CSparse.Ineropt C# project. All test example in different solvers (Cholmod, CSsparse, UMFPACK) work well. The problem is that large matrices in UMFPACK cause OutOfMemory:

Error CSparse.Interop.SuiteSparse.Umfpack.UmfpackException: ERROR: out of memory at CSparse.Interop.SuiteSparse.Umfpack.UmfpackContext1.Factorize() in C:...\CSparse.Interop\Interop\SuiteSparse\Umfpack\UmfpackContext.cs:line 60 at CSparse.Interop.SuiteSparse.Umfpack.UmfpackContext1.Solve(UmfpackSolve sys, T[] input, T[] result) in C:...\CSparse.Interop\Interop\SuiteSparse\Umfpack\UmfpackContext.cs:line 100

due to internal calling of int32 versions of UMFPACK. Is there a way how to switch to long(int64) version of UMFPACK within CSparse.Interopt? I use the following code to invert the matrix:

SparseMatrix A = new SparseMatrix(m_nDOF, m_nDOF, m_a, m_ja, m_ia); Umfpack umfpack = new Umfpack(A); umfpack.Solve(m_b, m_solution.P);

Thank you!

[wo80]

There are actually no plans to add support for the long version to this project, but it shouldn't be hard to get this working for yourself. A lot of copy and pasting of the existing code.

Before going any further, could you please debug the problem and step through the CSparse.Double.Factorization.SuiteSparse.Umfpack.DoFactorize() method. Check if the call to umfpack_di_symbolic succeeds and then report the contents of UmfpackInfo.Raw array.

[vstembera]

returned status = 0 info.Raw = [0] 0
[1] 1380959 [2] 131933819
[3] 8
[4] 4
[5] 8
[6] 8
[7] 8
[8] 0
[9] 0
[10] 0
[11] 0
[12] 0
[13] 241664959
[14] 3567871 [15] 0
[16] 1380959 [17] 0
[18] 3
[19] 1
[20] 37179628993 [21] 37889039179 [22] 528606457112171 [23] 13112955746 [24] 23308000592 [25] 350538215
[26] 37869453358 [27] 37170652715 [28] 478430757
[29] 16023
[30] 29919
[31] 1
[32] 1
[33] 1
[34] 130552860
[35] 1380959 [36] 765000299
[37] 933963677210
[38] 0
[39] 4605
[40] -1
[41] -1
[42] -1
[43] -1
[44] -1
[45] -1
[46] -1
[47] -1
[48] -1
[49] -1
[50] -1
[51] -1
[52] -1
[53] -1
[54] -1
[55] -1
[56] 0
[57] 0
[58] 1380959 [59] 1
[60] -1
[61] -1
[62] -1
[63] -1
[64] -1
[65] -1
[66] -1
[67] -1
[68] -1
[69] -1
[70] -1
[71] -1
[72] -1
[73] -1
[74] -1
[75] -1
[76] -1
[77] -1
[78] -1
[79] -1
[80] -1
[81] -1
[82] -1
[83] -1
[84] -1
[85] -1
[86] -1
[87] -1
[88] -1
[89] -1

[wo80]

Thanks. Lets look at the estimates computed by umfpack_di_symbolic. Number 20 is the final size estimated for Numeric->Memory (in bytes AFAIK): 37179628993 bytes > 34 GB. Are you sure this isn't an actual out of memory exception?

[wo80]

I've added the minimum required interop to support 64bit indices https://github.com/wo80/csparse-interop/commit/b7c055d518398d7b957fd5a5fcfcd90b9506b0fb

The code isn't tested, so you get no warranty at all. Let me know if this works!

[wo80]

Copypastas fixed in https://github.com/wo80/csparse-interop/commit/5a1fe6e98cb78e2f17c11e06cdd3b85c36d80760

[vstembera]

Thank you for new code, I am testing it! It seems however that new CompressedColumnStorage with long type of arrays is missing in the nuget CSparse.dll.

[wo80]

It seems however that new CompressedColumnStorage with long type of arrays is missing in the nuget CSparse.dll

Yes, there is no CompressedColumnStorage class with Int64 arrays and there won't be.

The problem here is not the 32-bit indexing of the input, but the fact that the number of non-zeros in the factorization exceeds the 32-bit range. To enable interop, I'm copying the Int32 storage arrays of the CompressedColumnStorage to Int64 arrays in the constructor of Umfpack_Long.

If you look at positions 23 and 24 of the UmfpackInfo array, the estimated numbers of non-zeros in L and U are 13,112,955,746 and 23,308,000,592. Those would exceed the max possible array length in .NET anyway, see docs.microsoft.com gcallowverylargeobjects.

[vstembera]

I have changed the calling to

Umfpack_Long umfpack = new Umfpack_Long(A); umfpack.Solve(m_b, m_solution.P);

However, the code still cannot find the corresponding method in the current dll:

_00:00:03.59, Error System.EntryPointNotFoundException: Unable to find an entry point named 'umfpack_dl_defaults' in DLL 'libumfpack'. at CSparse.Interop.SuiteSparse.Umfpack.NativeMethods.umfpack_dl_defaults(Double[] Control) at CSparse.Double.Factorization.SuiteSparse.Umfpack_Long.DoInitialize() in C:\D_A_T_A\PROJEKTY\A_femCalc\jobreader\CSparse.Interop\Double\Factorization\SuiteSparse\Umfpack_Long.cs:line 47 at CSparse.Interop.SuiteSparse.Umfpack.UmfpackContext1..ctor(CompressedColumnStorage1 matrix) in C:\D_A_T_A\PROJEKTY\A_femCalc\jobreader\CSparse.Interop\Interop\SuiteSparse\Umfpack\UmfpackContext.cs:line 43 at CSparse.Double.Factorization.SuiteSparse.Umfpack_Long..ctor(SparseMatrix matrix) in C:\D_A_T_A\PROJEKTY\A_femCalc\jobreader\CSparse.Interop\Double\Factorization\SuiteSparse\UmfpackLong.cs:line 24

I have checked libumpfpack.dll (from http://wo80.bplaced.net/math/packages.html) by Dependancy Walker and dl variants of the functions are really missing, ony di and zi variants are present.

[wo80]

dl variants of the functions are really missing, ony di and zi variants are present.

Because the dlls are intended to match the features of CSparse.Interop, as clearly stated:

Please note, that the packages were built to work with CSparse.Interop and that the list of exported functions might not be complete.

I'm not sure when I will find some time to recompile the libraries, maybe next weekend. All dependencies will have to include the long versions, like AMD. METIS reordering will definitely not work, since one has to make the choice of IDXTYPEWIDTH = 32 or 64 at compile time.

[vstembera]

Thanks very much for help! I am in the situation, where the matrices, which I get, are big and not positive definite at the same time. Therefore CSparse, Cholmod and other Cholesky solvers fail. Only one, which works, is Pardiso, which however is very bad with respect to round-off errors (it is well known, but I also have the same numerical experience with Pardiso). I know that UMFPACK works far better in these cases according to tests of my colleague, who however works in Linux and different software.

[wo80]

Have you considered using iterative solvers? MINRES if your problem is symmetric, otherwise GMRES or BiCGStab.

[vstembera]

Yes, it is also a possibility. Until now I tried to avoid it due to known problems with preconditioning, where as far as I know lots of heuristics takes place. However, do you know any good library with iterative solvers for large sparse matrices, which would be not difficult to connect to C# project? I think I should also try this way. (My matrices are symmetric.)

[wo80]

Until now I tried to avoid it due to known problems with preconditioning, where as far as I know lots of heuristics takes place.

Yes, the main problem is choosing a good preconditioner. You could try MathNet. It has a couple of iterative solvers and two preconditioners, diagonal and ilu0.

I've pushed MINRES to https://github.com/wo80/csparse-extensions/commit/185dd0d351cc4a4970bf551e3c0b184699e77911 . I haven't tested this in practice, but you could give it a try: MinResTest.cs

[vstembera]

I have just implemented MathNet.Numerics in my project. It was a good idea! It has 4 iterative solvers and 4 preconditioners, however only 2 combinations work well (GpBiCg+diagonal preconditioner and BiCgStab+diagonal preconditioner), other returns garbage or diverge at all. Interesting will be comparison with UMFPACK.

[wo80]

It's been some time since I used MathNet. The two preconditioners I mentioned that make sense for large matrices are Diagonal and MILU0. The other two - ILU0 and ILUTP - aren't optimized for sparse matrices.

I you find the time, take a look at MINRES. I'd really like to know, if it works. The setup is similar to MathNet.

[wo80]

I've added BiCgStab and MILU(0) with https://github.com/wo80/csparse-extensions/commit/e57929ecab383ec878dd27df9b77875aa5a2e062

MINRES expects a symmetric preconditioner, so don't try to use MILU(0) with it.

[vstembera]

Thank you. However I am now again very sceptic about iterative methods. From 16 possible combinations of MathNet.Numerics only 2 work on my test matrix with 2.5k dofs. On my original matrix with 1.38M dofs none of these 2 work and return only garbage. Moreover, it takes so long, from 30 minutes up to 6 hours, Pardiso solves the same system in 30 seconds.

[wo80]

I've added the 64-bit integer version to vs-suitesparse. Binaries available at http://wo80.bplaced.net/math/packages.html

Regarding iterative solvers, to get an acceptable running time, it's crucial to choose a good preconditioner. Diagonal preconditioning usually doesn't help much, so the only one making sense here is MILU(0). In general, I'd choose an iterative over a direct solver only if the factorization exceeds your hardware capabilities.

[vstembera]

The dll works! Perfect. I check also the MINRES, when I find some time and I let you know.

P.S. For really large matrices in UMFPACK (I have tried 28797598 unknowns) I got an overflow problem, but that is probably the real boundary of UMFPACK as such.
Error System.OverflowException: Arithmetic operation resulted in an overflow. at MathLib.SparseMatrixStorageCalculus.Create(Int32 nDOF, Int32[][] incidenceMatrix, Int32[,] DOFMatrix, Double[,] LagrangeMultiplierMatrix, DirectLinearSolverType solverType) in C:\D_A_T_A\PROJEKTY\A_femCalc\jobreader\Solver\MathLib\SparseMatrixStorageCalculus.cs:line 182

[wo80]

Probably not UMFPACK, but the hard limit for array sizes in .NET, see https://github.com/wo80/csparse-interop/issues/11#issuecomment-1061120822

I think you can forget about MINRES at the moment, since there's no useful preconditioner. I'm currently adding HYPRE interop to this project, so you might want to check that out when it's ready: https://github.com/wo80/csparse-interop/pull/12

I'll also publish the native libraries, which will contain the sequential interface (no MPI or OpenMP).