Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Return Q/R factors in SPQR #394

Closed
jagomezr opened this issue Jan 14, 2017 · 16 comments · Fixed by JuliaLang/julia#20464
Closed

Return Q/R factors in SPQR #394

jagomezr opened this issue Jan 14, 2017 · 16 comments · Fixed by JuliaLang/julia#20464
Labels
help wanted Extra attention is needed sparse Sparse arrays

Comments

@jagomezr
Copy link

Hi,

Has someone added the functionality to extract the Q and R factors of the SPQR factorization? So far, the result of using qrfact on a sparse matrix is the following:

Base.SparseArrays.SPQR.Factorization{Float64}(3,3,Ptr{Base.SparseArrays.SPQR.C_Factorization{Float64}} @0x0000000009f9d4a0)

I see that this issue has been reported in the past, but it seems no one has added this functionality to Julia.

Thanks!

@tkelman tkelman added linear algebra sparse Sparse arrays help wanted Extra attention is needed labels Jan 14, 2017
@jagomezr
Copy link
Author

In general, I have a system of the form:
A x = b

where A is a rectangular matrix. I want to obtain an equivalent system that is full row rank. The way I'm doing it is calculating the rank of A, and then using the QR = AE factorization and the rank r to find an equivalent system made of the first r rows of:

RE' x = Q'b.

To do this I'm using the following code:

F = qrfact(model.S);
Sfull = full(model.S);
r = rank(Sfull);
bnew = Base.SparseArrays.SPQR.qmult(Base.SparseArrays.SPQR.QTX, F, Base.SparseArrays.CHOLMOD.Dense(model.b));
Anew = sparse(Base.SparseArrays.SPQR.qmult(Base.SparseArrays.SPQR.QTX, F, Base.SparseArrays.CHOLMOD.Dense(Sfull)));
r1 = rank(full(Anew));
if r1 != r
error("The full rank code is changing the original rank of the matrix");
end
model.S = Anew[1:r,:];
model.b = bnew[1:r];

However, I'm getting the error below. Any ideas?

Thanks!

*** Error in `julia': double free or corruption (!prev): 0x00000000050961d0 ***

signal (6): Aborted
while loading no file, in expression starting on line 0
raise at /build/buildd/eglibc-2.19/signal/../nptl/sysdeps/unix/sysv/linux/raise.c:56
abort at /build/buildd/eglibc-2.19/stdlib/abort.c:89
__libc_message at /build/buildd/eglibc-2.19/libio/../sysdeps/posix/libc_fatal.c:175
malloc_printerr at /build/buildd/eglibc-2.19/malloc/malloc.c:4996 [inlined]
_int_free at /build/buildd/eglibc-2.19/malloc/malloc.c:3840
jl_gc_counted_free at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:1958
SuiteSparse_free at /home/jagomezr/Julia050/bin/../lib/julia/libsuitesparseconfig.so (unknown line)
cholmod_l_free at /home/jagomezr/Julia050/bin/../lib/julia/libcholmod.so (unknown line)
cholmod_l_free_sparse at /home/jagomezr/Julia050/bin/../lib/julia/libcholmod.so (unknown line)
free_sparse! at ./sparse/cholmod.jl:465
free! at ./sparse/cholmod.jl:1069
unknown function (ip: 0x7f1f4a0c66b2)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
jl_apply at /home/centos/buildbot/slave/package_tarball64/build/src/julia.h:1392 [inlined]
run_finalizer at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:108
jl_gc_run_finalizers_in_list at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:190 [inlined]
run_finalizers at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:211
jl_gc_enable_finalizers at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:222
jl_mutex_unlock at /home/centos/buildbot/slave/package_tarball64/build/src/./julia_threads.h:476 [inlined]
jl_compile_for_dispatch at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1314
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:184 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./REPL.jl:132
unknown function (ip: 0x7f1f4a0d09e6)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./REPL.jl:135
unknown function (ip: 0x7f1f4a0d06b6)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./multimedia.jl:143
unknown function (ip: 0x7f1f4a0d0502)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_response at ./REPL.jl:154
unknown function (ip: 0x7f1f4a090eb8)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_response at ./REPL.jl:139
unknown function (ip: 0x7f1f4a090938)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
JuliaLang/julia#22 at ./REPL.jl:652
unknown function (ip: 0x7f1f4a08f931)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
JuliaLang/julia#38 at ./REPL.jl:867
JuliaLang/julia#13 at ./LineEdit.jl:736
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
prompt! at ./LineEdit.jl:1605
run_interface at ./LineEdit.jl:1574
unknown function (ip: 0x7f21645390bf)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
run_frontend at ./REPL.jl:903
run_repl at ./REPL.jl:188
unknown function (ip: 0x7f1f4a089532)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
_start at ./client.jl:360
unknown function (ip: 0x7f21645542e8)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
unknown function (ip: 0x4018ed)
unknown function (ip: 0x4013b6)
__libc_start_main at /build/buildd/eglibc-2.19/csu/libc-start.c:287
unknown function (ip: 0x4013fc)
Allocations: 5872265 (Pool: 5871060; Big: 1205); GC: 17
Aborted (core dumped)

@tkelman
Copy link

tkelman commented Jan 15, 2017

That looks like a bug, but separate from your feature request here. Can you provide self contained input data that reproducibly causes that? And platform and version info.

@jagomezr
Copy link
Author

jagomezr commented Jan 15, 2017

You can try the following code. This is in Julia 0.5.0 in Ubuntu 14.04.1 LTS and a 64 bit machine.

test.txt

@andreasnoack
Copy link
Member

andreasnoack commented Jan 16, 2017

You are not supposed to free the SQPR object. Julia will do that for you so Julia segfaults when you have already freed the object. If I uncomment the free! line then your code works fine.

@jagomezr
Copy link
Author

If I comment that line (I do not use free), I get the following output:

*** Error in `julia': double free or corruption (!prev): 0x000000000392f1b0 ***

signal (6): Aborted
while loading no file, in expression starting on line 0
raise at /build/buildd/eglibc-2.19/signal/../nptl/sysdeps/unix/sysv/linux/raise.c:56
abort at /build/buildd/eglibc-2.19/stdlib/abort.c:89
__libc_message at /build/buildd/eglibc-2.19/libio/../sysdeps/posix/libc_fatal.c:175
malloc_printerr at /build/buildd/eglibc-2.19/malloc/malloc.c:4996 [inlined]
_int_free at /build/buildd/eglibc-2.19/malloc/malloc.c:3840
jl_gc_counted_free at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:1958
SuiteSparse_free at /home/jagomezr/Julia050/bin/../lib/julia/libsuitesparseconfig.so (unknown line)
cholmod_l_free at /home/jagomezr/Julia050/bin/../lib/julia/libcholmod.so (unknown line)
cholmod_l_free_sparse at /home/jagomezr/Julia050/bin/../lib/julia/libcholmod.so (unknown line)
free_sparse! at ./sparse/cholmod.jl:465
free! at ./sparse/cholmod.jl:1069
unknown function (ip: 0x7f2cf7bf28b2)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
jl_apply at /home/centos/buildbot/slave/package_tarball64/build/src/julia.h:1392 [inlined]
run_finalizer at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:108
jl_gc_run_finalizers_in_list at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:190 [inlined]
run_finalizers at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:211
jl_gc_enable_finalizers at /home/centos/buildbot/slave/package_tarball64/build/src/gc.c:222
jl_mutex_unlock at /home/centos/buildbot/slave/package_tarball64/build/src/./julia_threads.h:476 [inlined]
jl_compile_for_dispatch at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1314
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:184 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./multimedia.jl:143
unknown function (ip: 0x7f2cf7bfbb82)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_response at ./REPL.jl:154
unknown function (ip: 0x7f2cf7bc1e08)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_response at ./REPL.jl:139
unknown function (ip: 0x7f2cf7bc1888)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
JuliaLang/julia#22 at ./REPL.jl:652
unknown function (ip: 0x7f2cf7bbd931)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
JuliaLang/julia#38 at ./REPL.jl:867
JuliaLang/julia#13 at ./LineEdit.jl:736
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
prompt! at ./LineEdit.jl:1605
run_interface at ./LineEdit.jl:1574
unknown function (ip: 0x7f2f120670bf)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
run_frontend at ./REPL.jl:903
run_repl at ./REPL.jl:188
unknown function (ip: 0x7f2cf7bb7532)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
_start at ./client.jl:360
unknown function (ip: 0x7f2f120822e8)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
unknown function (ip: 0x4018ed)
unknown function (ip: 0x4013b6)
__libc_start_main at /build/buildd/eglibc-2.19/csu/libc-start.c:287
unknown function (ip: 0x4013fc)
Allocations: 5324334 (Pool: 5323189; Big: 1145); GC: 14
Aborted (core dumped)

@andreasnoack
Copy link
Member

In this file line 75 has a free! which causes a segfault on my machine. If you are getting a segfault with a different program then please post a link to that.

@jagomezr
Copy link
Author

The test.jl in the attached folder fails with
*** Error in `julia': double free or corruption (!prev): 0x0000000004a2fc70 ***

SPQR failure.zip

@andreasnoack
Copy link
Member

Thanks. Just tried it on my Mac. It works on both master and 0.5. What is your versioninfo()?

@jagomezr
Copy link
Author

Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
System: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.7.1 (ORCJIT, ivybridge)

@andreasnoack
Copy link
Member

Just tried on a Linux machine with 0.5 and now I can reproduce the segfault. I'll look into it.

@jagomezr
Copy link
Author

Thanks for looking into this!

For the record, the same code in Julia 0.5.0 in Windows 10 also fails. No error is returned, the terminal just shuts down.

julia> versioninfo()
Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
System: NT (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.7.1 (ORCJIT, ivybridge)

@andreasnoack
Copy link
Member

I've spent some time on this but so far I haven't been able to figure out what the problem is including why this is not showing up on Mac.

@RalphAS
Copy link

RalphAS commented Feb 4, 2017

The failure in the example "SPQR failure.zip" posted above is due to a malformed input: one of the matrices returned by matread has more indices than values, which leads to heap corruption after SuiteSparse frees its working arrays. Whether this heap corruption causes a crash seems to be platform dependent. I suspect a bug in MAT.jl or a flaw in the data file.

With regard to the request for Q and R matrices, they are not immediately accessible from the C API which Julia uses because they are wasteful for most purposes. One could of course apply qmult to the identity.

@andreasnoack
Copy link
Member

See the discussion in JuliaLang/julia#20464. It appears that we actually made some wrong assumptions about the buffer sizes. The PR fixes the issue on the Linux machine where I was able to reproduce the problem.

@jagomezr
Copy link
Author

What does PR mean? I see you committed some changes that you are waiting to merge with master. How can I adopt the changes? Should I wait until the changes are pulled into the master?

Thanks!

@andreasnoack
Copy link
Member

PR means pull request. There was a request for an extra test which I haven't added yet. It should happen soon.

@KristofferC KristofferC transferred this issue from JuliaLang/julia Nov 26, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
help wanted Extra attention is needed sparse Sparse arrays
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants