JuliaSmoothOptimizers / QRMumps.jl

Interface to multicore QR factorization qr_mumps
GNU Lesser General Public License v3.0
16 stars 6 forks source link

Arguments check in COO input #104

Open dpo opened 2 months ago

dpo commented 2 months ago

The following code works

julia> rows = [1, 2, 3, 4]; cols = [1, 2, 3, 4]; vals = ones(4);

julia> qrm_spmat_init(2, 2, rows, cols, vals)
Sparse matrix -- qrm_spmat of size (2, 2) with 4 nonzeros.

but shouldn’t.

Where should the argument check happen? In Julia or in Fortran?

@abuttari

dpo commented 2 months ago

Also is it correct and intentional that duplicate COO entries are not accepted by QRMumps?

abuttari commented 2 months ago

The following code works

julia> rows = [1, 2, 3, 4]; cols = [1, 2, 3, 4]; vals = ones(4);

julia> qrm_spmat_init(2, 2, rows, cols, vals)
Sparse matrix -- qrm_spmat of size (2, 2) with 4 nonzeros.

but shouldn’t.

Where should the argument check happen? In Julia or in Fortran?

@abuttari

this variant of qrm_spmat_init does not correspond to any qr_mumps routine but is defined

https://github.com/JuliaSmoothOptimizers/QRMumps.jl/blob/e45811442a0f7e555fd5def595b0d22b59adab2d/src/wrapper/qr_mumps_api.jl#L68

However I am not sure this should be fixed. In this case qr_mumps should consider the last two coefficients as out-of-bounds and ignore them.

abuttari commented 2 months ago

Also is it correct and intentional that duplicate COO entries are not accepted by QRMumps?

@dpo yes, they are accepted. They are added when frontal matrices are assembled during the factorization. Are you getting some inconsistent behavior?

amontoison commented 2 months ago

@abuttari https://github.com/JuliaSmoothOptimizers/QRMumps.jl/issues/101

abuttari commented 2 months ago

@abuttari #101

sorry I overlooked this issue. I am not sure what's happening but if I read this matrix directly in the timing/qrm_spgeqr.F90 benchmark of qr_mumps (i.e., not passing through the Julia interface) everything work fine

amontoison commented 2 months ago

Can you try with the C interface but passing through the Julia interface?

abuttari commented 2 months ago

yes, I can reproduce the error. I am supposed to see messages printed within the qr_mumps code? for example I have a write statement inside qrm_spfct_init (where no error occurs) but it doesn't get printed in julia

dpo commented 2 months ago

Do you mean that an error is printed related to out-of-bound indices?

abuttari commented 2 months ago

no, I mean that I added print statements inside qr_mumps to see how far the execution gets but the messages are not printed when I run qr_mumps from within julia. Actually I am not 100% sure I am really using the qr_mumps I have manually built. All I have to do is

export JULIA_QRMUMPS_LIBRARY_PATH=...

before launching Julia, right?

amontoison commented 2 months ago

It's right Alfredo, you can check what is installed with QRMumps.QRMUMPS_INSTALLATION. If needed, you can force a recompilation to use your local version with:

force_recompile(package_name::String) = Base.compilecache(Base.identify_package(package_name))
force_recompile("QRMumps")
using QRMumps
QRMumps.QRMUMPS_INSTALLATION
dpo commented 2 months ago

I think so. You can also check that QRMumps.libqrm_common points to the right library.

In my view, it would be more useful to users if out-of-bound indices raised an error instead of being silently ignore. But I could easily add that to the Julia interface if you think it doesn’t belong in the Fortran library.

abuttari commented 2 months ago

ok, I forced the recompilation and verified that QRMumps.libqrm_common points to the good library. Then I do

using QRMumps, DelimitedFiles
m = 361
n = 802

rows = readdlm("rows.txt",Int32)[:]
cols = readdlm("cols.txt",Int32)[:]
vals = readdlm("vals.txt")[:]

spmat = qrm_spmat_init(m,n,rows,cols,vals; sym=false)
spfct = qrm_spfct_init(spmat)
qrm_analyse!(spmat, spfct; transp='t')

and it goes to completion without any error. With the YGGDRASIL install, instead, I have

free(): corrupted unsorted chunks

With the custom installation though I have another strange behavior: as soon as I execute qrm_init() Julia goes up to 100% CPU utilization and doesn't come back down .

amontoison commented 2 months ago

The call to qrm_init() is mandatory before running any code with the custom installation.

I do it automatically with the Yggdrasil installation because It's the sequential version and I know the number of CPU (one) and GPU (zero).