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 b7c055d

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

[wo80]

Copypastas fixed in 5a1fe6e

[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\Umfpack_Long.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 wo80/csparse-extensions@185dd0d . 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.