Bug in UMFPack Interface

[cortner]

I had some unexplained \ behaviour that is, according to Ivar Nesje, a bug in /base/linalg/umfpack.jl#L18

Here is a short code snippet that works ok. I get the correct solution from it. Primarily the last line is relevant. ( I can produce the complete code on request. )

 # construct a finite element grid
 X, T = square(10)
 # assemble the stiffness matrix and rhs (for - \Delta u = 1)
 A, F = simple_stiffmat(X, T)
 # find the interior nodes
 Ifree = find( (minimum(X,1) .> 0) & (maximum(X, 1) .< 1) )
 # solve the linear system and plot the solution
 U = zeros(size(X,2))
 U[Ifree] = cholfact(A[Ifree,Ifree]) \ F[Ifree]

If I replace the last line with

 U[Ifree] = A[Ifree,Ifree] \ F[Ifree]

then I get the following error message:

umferror has no method matching umferror(::Int32)
while loading In[112], in expression starting on line 3

in umfpack_symbolic! at linalg/umfpack.jl:155
in umfpack_numeric! at linalg/umfpack.jl:176
in lufact at linalg/umfpack.jl:118
in A_ldiv_B! at linalg/sparse.jl:210
in \ at linalg/generic.jl:233

Same happens with lufact. I also tried copying B = copy(A[Ifree,Ifree]) first. Same result.

I have used UMFPack successfully in other cases. E.g.,

 U = (A+speye(size(X,2)) \ F

works fine.

I copy Ivar Nesje's comment from Google Groups:
This is at least one bug.

In /base/linalg/umfpack.jl#L18, we have a convenience function to lookup the status code, seems to be restricted to Int64, even though it is called with Int32 on some occations. If you compile yourself, or know how to regenerate the system image, it is pretty easy to change Int to Integer on that line. A real fix for this would require someone to actually check the return types of the C functions called, and include some more tests to ensure that all code paths are followed in the test suite.

[ivarne]

Link to google groups discussion: https://groups.google.com/forum/#!topic/julia-users/_3PLfz3vMZk

[ViralBShah]

Oops, sorry @andreasnoack - Assigned to you by mistake while trying to set the label. The github android interface is not the greatest...

[tkelman]

andreasnoack, linear algebra, same thing

[andreasnoack]

@cortner Is it possible for you to provide a reproducible code example for this problem? It makes it much easier to fix.

[cortner]

This is from a little FEM code I wrote - I will reduce it to a bare minimum and make it available as a notebook.

[cortner]

@andreasnoack I uploaded a notebook with the code on my homepage:
FEM-Minimal or open with nbviewer.

The second cell are the matrix assembly routines. Then in the subsequent cells I have some examples that work and one that does not. I hope this help.

[andreasnoack]

The "THIS DOES NOT WORK" part works for me. What is your versioninfo()

[cortner]

Julia Version 0.3.2
Commit 21d5433* (2014-10-21 20:18 UTC)
Platform Info:
System: Darwin (x86_64-apple-darwin13.3.0)
CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
WORD_SIZE: 64
BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Sandybridge)
LAPACK: libopenblas
LIBM: libopenlibm
LLVM: libLLVM-3.3

Strange, I definitely updated to 0.3.3 - I will check to see what is going on and then post again.

[cortner]

After updating IJulia to 3.3, it now does indeed work. I am very sorry for having wasted your time. (In my defense - a few days ago it was still a bug.)

Out of curiosity: what changed between versions?

Christoph

[andreasnoack]

No problem. It didn't take long time to check your example. Actually, I don't know which change that you have affected this.

[ivarne]

@cortner see JuliaLang/julia@v0.3.2...v0.3.3

[ViralBShah]

@ivarne That is pretty cool. Didn't know about that.

[cortner]

So I ran into the problem again. I updated the notebook with a reproducible example:
View
Download

I did not get the problem when the system matrix was symmetric, but I did get it when it was non-symmetric. So I suspect that Julia v.3.3 checks for symmetry and then uses a different solver.

@andreasnoack would you be willing to double-check this, and
@ViralBShah if you can confirm it, can the issue be reopened?

I would also be interested if there is a temporary workaround?

Thanks, Christoph

[ivarne]

@andreasnoack How about tests?

cc: @JuliaBackports

[tkelman]

Already backported apparently 6c995fedefea316d94f630c622e3d990dfe49355

Tests would be a very good idea

[andreasnoack]

Thanks for the -x trick. I'll write some tests. Apparently, the UMFPACK wrappers haven't been tested for 32 bit integers on 64 bit systems. The error in this issue also seems to have hidden another one. @tkelman, have you seen this one before

julia> using Base.LinAlg.UMFPACK.increment!

julia> A = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
                  increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
                  [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.], 5, 5)
5x5 sparse matrix with 12 Float64 entries:
    [1, 1]  =  2.0
    [2, 1]  =  3.0
    [1, 2]  =  3.0
    [3, 2]  =  -1.0
    [5, 2]  =  4.0
    [2, 3]  =  4.0
    [3, 3]  =  -3.0
    [4, 3]  =  1.0
    [5, 3]  =  2.0
    [3, 4]  =  2.0
    [2, 5]  =  6.0
    [5, 5]  =  1.0

julia> lufact(convert(SparseMatrixCSC{Float64,Int32}, A))
UMFPACK LU Factorization of a 5-by-5 sparse matrix
Ptr{Void} @0x00007fcccf99def0


julia> gc()

julia> gc()
julia(39834,0x7fff7a133300) malloc: *** error for object 0x2: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug

signal (6): Abort trap: 6
__pthread_kill at /usr/lib/system/libsystem_kernel.dylib (unknown line)
Abort trap: 6

[tkelman]

My goodness. Nope. Sort of like JuliaLang/julia#9185 in that it appears to have something to do with the umfpack finalizers?