bugGNU Octave - Bugs: bug #57033, Replace CXSPARSE with SPQR

 
 

bug #57033: Replace CXSPARSE with SPQR

Submitted by:  Rik <rik5>
Submitted on:  Thu 10 Oct 2019 04:57:19 PM UTC  
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Feature Request
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Sun 29 Mar 2020 10:41:06 AM UTC, comment #16: 

The test suite is crashing with a segfault on this test:

%!testif HAVE_CXSPARSE
%! n = 20;  d = 0.2;
%! ## initialize generators to make behavior reproducible
%! rand ("state", 42);
%! randn ("state", 42);
%! a = sprandn (n,n,d) + speye (n,n);
%! [c,r] = qr (a, ones (n,1));
%! assert (r\c, full (a)\ones (n,1), 10e-10);

Running the test alone doesn't crash Octave for me.
Minimal steps to reproduce for me are:

a = [0, 2, 1; 2, 1, 2];
[q, r] = qr (a);
n = 20;  d = 0.2;
a = sprandn (n,n,d) + speye (n,n);
[c,r] = qr (a, ones (n,1));

Digging a little bit with gdb, the segfault occurs on the attempt to free "cc->Flag".
It could be that something isn't correctly cleaned up for the dense case. Or some initialization is missing for the sparse case.
I'm no suitesparse expert though.

Wrt the testif syntax: Are you sure the complex syntax with "||" works before the semicolon? I only found it used after the semicolon.
I could only find examples where HAVE_* feature tests are combined with a comma, which seems to imply an "&&".
Testing the syntax, any symbol can be used to separate the HAVE_* features, and they all are interpreted as "&&" (including "||").

Markus Mützel <mmuetzel>
Project Member
Sat 28 Mar 2020 08:24:57 PM UTC, comment #15: 

I think from Octave's perspective the underlying library shourdn't matter.  Hence, the BIST tests for qr() or whatever function should be executed if either CXSPARSE or SPQR is available.

The %!test syntax supports complex expressions such as

%!testif HAVE_CXSPARSE || HAVE_SPQR

But, first thing is that the code shouldn't segfault.

Rik <rik5>
Project Administrator
Sat 28 Mar 2020 12:47:03 PM UTC, comment #14: 

Thanks for the update.

The attached patch includes some minor style review (mostly white spaces) and restores the lost commit message.

It compiles without errors. But crashes during "make check".
Running the test suite with gdb ("make check RUN_OCTAVE_OPTIONS=-g"), I see the following top of the back trace when the crash occurs:

(gdb) bt
#0  0x00007ffff518a8a0 in __GI___libc_free (mem=0xffff00001fbf) at malloc.c:3102
#1  0x00007ffff2e2d733 in SuiteSparse_free () at /usr/lib/x86_64-linux-gnu/libsuitesparseconfig.so.5
#2  0x00007ffff2ef78ea in cholmod_l_free () at /usr/lib/x86_64-linux-gnu/libcholmod.so.3
#3  0x00007ffff2ef281f in cholmod_l_free_work () at /usr/lib/x86_64-linux-gnu/libcholmod.so.3
#4  0x00007ffff5e742da in octave::math::sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep() (this=0x7fffcd0a6090, __in_chrg=<optimized out>) at ../liboctave/numeric/sparse-qr.cc:346
#5  0x00007ffff5e7ffb1 in octave::math::sparse_qr<SparseMatrix>::~sparse_qr() (this=<optimized out>, __in_chrg=<optimized out>) at ../liboctave/numeric/sparse-qr.cc:2782
#6  0x00007ffff5e7ffb1 in octave::math::sparse_qr<SparseMatrix>::~sparse_qr() (this=<optimized out>, __in_chrg=<optimized out>) at ../liboctave/numeric/sparse-qr.cc:2782
#7  0x00007ffff774c2ab in Fqr(octave_value_list const&, int) (args=..., nargout=<optimized out>) at ../liboctave/array/MSparse.h:74
#8  0x00007ffff70592b0 in octave_builtin::call(octave::tree_evaluator&, int, octave_value_list const&) (this=<optimized out>, tw=..., nargout=2, args=...)
    at ../libinterp/octave-value/ov-builtin.cc:65
#9  0x00007ffff7248141 in octave::tree_index_expression::evaluate_n(octave::tree_evaluator&, int) (this=0x7fffcd0bf9d0, tw=..., nargout=2) at /usr/include/c++/9/ext/new_allocator.h:89

This is with Ubuntu 19.10 and suitesparse 5.4.0.

One more doubt I have: Should the BISTs that depend on CXSPARSE be replicated for SPQR so that they run on systems with SPQR but without CXSPARSE?

(file #48691)

Markus Mützel <mmuetzel>
Project Member
Fri 13 Mar 2020 10:53:10 PM UTC, comment #13: 

added documentation of new features to doc string

(file #48596)

Simon <sim79>
Sun 23 Feb 2020 02:25:51 PM UTC, comment #12: 

"qr.cc" was moved from dldfcn to corefcn. That caused the patch to no longer apply cleanly.
The attached patch is rebased on the current tip of the default branch. It also cleans up some minor issues (like reverted copyright blocks).

It compiles for me and all BISTs pass.

Could you please add some documentation of the new features to the doc string of qr?

(file #48481)

Markus Mützel <mmuetzel>
Project Member
Fri 10 Jan 2020 03:39:18 PM UTC, comment #11: 

This is good.  I don't think I can fit it in to the 6.1 release, but immediately afterwards.

Rik <rik5>
Project Administrator
Fri 10 Jan 2020 12:22:23 AM UTC, comment #10: 

I appended a new patch (it contains rev2-patch).

Changes:
Economy versions of C, Q and R are returned if qr(...,0) is called.
Column permutation of A is returned as (vector|matrix) if qr(...,("vector"|"matrix")) is called.
Additional tests if HAVE_SPQR is set.

(file #48196)

Simon <sim79>
Mon 30 Dec 2019 01:47:45 PM UTC, comment #9: 

I appended bug57033_spqr_rev2.patch. Tests in qr.cc-tst and sparse.tst are passed, but some warnings are issued.

Changes:
1) Removed fill-reducing ordering for calls [R] = qr(A), [Q,R] = qr(A) and [C,R] = qr(A,B).
2) If last argument of qr is "vector", a permutation vector is returned. Otherwise a permutation matrix is returned.

(file #48139, file #48140)

Simon <sim79>
Sat 14 Dec 2019 01:29:10 AM UTC, comment #8: 

It would be nice to have a review soon.  Otherwise, it won't be able to make it in to the 6.1 release.

Rik <rik5>
Project Administrator
Sun 01 Dec 2019 04:03:54 PM UTC, comment #7: 

I cleaned up your patch (mostly white space changes) and wanted to push it.
To silence the shadowing warnings, I renamed the local variable "S" and removed the property "econ" that doesn't seem to be used.

Unfortunately with the patch, some of the tests for "qr" fail. I attached a log file with the failing tests. Could you please have a look?
Could you also please base any changes on the attached patch? It took a considerable amount of time to get there.

Do you want to be added to the list of contributors at the beginning of the manual?

(file #47957, file #47958)

Markus Mützel <mmuetzel>
Project Member
Fri 18 Oct 2019 09:16:30 PM UTC, comment #6: 

Hello,
here are some additional notes on the test script.

Let A be a sparse m x n matrix. The test script is divided into four sections.
In the first section the QR-factorization [Q,R,E] = qr(A) is calculated. E is an array of size n that represents a (fill-reducing) column permutation of A. The script creates a permutation matrix P from E (with Q*R = A*P) and checks whether Q*R*P' - A is small. 
In the second section some overdetermined systems A*x = b (m > n) are solved by calling x = qr(A,b). The script prints norm (A*x - b,2). Usually the values of norm(A*x - b,2) are between 1.0 and 6.0. The n columns of A span a subspace of the m-dimensional vector space. If b lies outside this subspace, there is no x with A * x = b. Only an approximation A*x of b can be calculated.
In the third section some underdetermined systems A*x = b (m < n) are solved by calling x = qr(A,b). The script prints norm(A*x - b,2). Usually the values of norm(A*x - b,2) have the order of magnitude -16. If m < n, there is a high chance that a column-subset of A spans the complete m-dimensional vector space. In this case, there exists x with A*x = b and norm(A*x - b,2) is very small. If no such subset exists, the value of norm (A*x - b,2) is large (> 0.1).
The fourth section checks whether the QR-factorization for a dense A still works.

Simon

Simon <sim79>
Sun 13 Oct 2019 03:14:43 PM UTC, comment #5: 

Thanks, Rik.

The attached patch adds a test for -lspqr to configure. However, I don't know how to check for a function from a template. Instead I checked for the C function "SuiteSparseQR_C" which happens to live in the same library.
Also some more logic could probably be added to make sure that we don't link to CXSparse in case SPQR can be used. (If this is what we want to do.)
There is a comment in configure.ac that the order of the suitesparse flags matters on some systems. The current order works on my system (Ubuntu). I don't know what the "correct" order is though.

The patch applies on top of the patch from comment #0. Feel free to merge the patches if you want to.

During compilation, I see a few shadowed declaration warnings:

../liboctave/numeric/sparse-qr.cc: In member function ‘SPARSE_T octave::math::sparse_qr<SPARSE_T>::sparse_qr_rep::R(bool) const [with SPARSE_T = SparseMatrix]’:
../liboctave/numeric/sparse-qr.cc:371:59: warning: declaration of ‘econ’ shadows a member of ‘octave::math::sparse_qr<SparseMatrix>::sparse_qr_rep’ [-Wshadow]
     sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const
                                                           ^~~~~
../liboctave/numeric/sparse-qr.cc:130:23: note: shadowed declaration is here
       octave_idx_type econ;
                       ^~~~
../liboctave/numeric/sparse-qr.cc: In member function ‘RET_T octave::math::sparse_qr<SPARSE_T>::sparse_qr_rep::wide_solve(const RHS_T&, octave_idx_type&) const [with RHS_T = SparseMatrix; RET_T = SparseMatrix; SPARSE_T = SparseMatrix; octave_idx_type = long int]’:
../liboctave/numeric/sparse-qr.cc:1059:23: warning: declaration of ‘S’ shadows a member of ‘octave::math::sparse_qr<SparseMatrix>::sparse_qr_rep’ [-Wshadow]
       cholmod_sparse *S;
                       ^
../liboctave/numeric/sparse-qr.cc:132:57: note: shadowed declaration is here
       typename cxsparse_types<SPARSE_T>::symbolic_type *S;
                                                         ^
../liboctave/numeric/sparse-qr.cc: In member function ‘SPARSE_T octave::math::sparse_qr<SPARSE_T>::sparse_qr_rep::R(bool) const [with SPARSE_T = SparseComplexMatrix]’:
../liboctave/numeric/sparse-qr.cc:1542:66: warning: declaration of ‘econ’ shadows a member of ‘octave::math::sparse_qr<SparseComplexMatrix>::sparse_qr_rep’ [-Wshadow]
     sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
                                                                  ^~~~~
../liboctave/numeric/sparse-qr.cc:130:23: note: shadowed declaration is here
       octave_idx_type econ;
                       ^~~~

The tests from the .zip file also produce some results and don't crash for me. I'm also not sure how to interpret the results though.

I'll try to cross-build for Windows next.

(file #47661)

Markus Mützel <mmuetzel>
Project Member
Sat 12 Oct 2019 11:48:21 PM UTC, comment #4: 

So I can get it to compile by adding "-lspqr" to the environment variable LDFLAGS and then re-configuring and compiling.

This suggests only that some more work has to be done with the configuration and build system to check for the existence of libspqr during configuration, and then to add "-lspqr" to the LIBS variable.

With the compiled version I was able to run test_spqr.m in the zip file.  I don't know what output it is supposed to produce, but it doesn't segfault or do anything bad.

Rik <rik5>
Project Administrator
Sat 12 Oct 2019 11:18:47 PM UTC, comment #3: 

I also fail to compile on Linux.  This is with suitesparse-dev-5.1.2.  It is failing at the link stage.  It may be that the SQPR routines are in a different package such as libspqr2 (but I have that installed), or that some different link flags need to be set.

Rik <rik5>
Project Administrator
Sat 12 Oct 2019 10:47:42 AM UTC, comment #2: 

I tried to compile with this patch (to eventually try with a cross-build on Windows). However, the linker fails to resolve some references:

  CXXLD    src/octave-cli
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_dense_struct* SuiteSparseQR_qmult<double>(int, cholmod_sparse_struct*, cholmod_dense_struct*, long*, cholmod_dense_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_sparse_struct* SuiteSparseQR_min2norm<double>(int, double, cholmod_sparse_struct*, cholmod_sparse_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `long SuiteSparseQR<double>(int, double, long, cholmod_sparse_struct*, cholmod_sparse_struct**, long**, cholmod_sparse_struct**, long**, cholmod_dense_struct**, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `long SuiteSparseQR<std::complex<double> >(int, double, long, cholmod_sparse_struct*, cholmod_sparse_struct**, long**, cholmod_sparse_struct**, long**, cholmod_dense_struct**, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_sparse_struct* SuiteSparseQR_qmult<double>(int, cholmod_sparse_struct*, cholmod_dense_struct*, long*, cholmod_sparse_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_dense_struct* SuiteSparseQR_qmult<std::complex<double> >(int, cholmod_sparse_struct*, cholmod_dense_struct*, long*, cholmod_dense_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_dense_struct* SuiteSparseQR_min2norm<std::complex<double> >(int, double, cholmod_sparse_struct*, cholmod_dense_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_sparse_struct* SuiteSparseQR_min2norm<std::complex<double> >(int, double, cholmod_sparse_struct*, cholmod_sparse_struct*, cholmod_common_struct*)'
/usr/bin/ld: liboctave/.libs/liboctave.so: undefined reference to `cholmod_dense_struct* SuiteSparseQR_min2norm<double>(int, double, cholmod_sparse_struct*, cholmod_dense_struct*, cholmod_common_struct*)'
collect2: error: ld returned 1 exit status
make[2]: *** [Makefile:15663: src/octave-cli] Error 1
make[2]: Leaving directory '/home/osboxes/Repositories/Octave/octave-1/.build'
make[1]: *** [Makefile:27400: all-recursive] Error 1
make[1]: Leaving directory '/home/osboxes/Repositories/Octave/octave-1/.build'
make: *** [Makefile:11032: all] Error 2

My steps were:
1) Update octave repo to hg id 3fe26656e73c
2) Apply the patch from comment #0.
3) Run "make maintainer-clean"
4) Complete "bootstrap"-"configure"-"make all" cycle.

This was on Ubuntu 19.04 with libsuitesparse-dev 5.4.0.
I understand that this is an incomplete patch. Please, let me know if it is too early to try it on a cross-build for Windows.
Otherwise: What do I have to do to get this to compile and link (on Ubuntu-Linux as a first step)?

Markus Mützel <mmuetzel>
Project Member
Thu 10 Oct 2019 04:58:44 PM UTC, comment #1: 

I filed this issue report to keep track of a promising patch to implement SPQR for sparse operations in Octave.  I've added the author of the patch (Simon) to the CC list for this bug, as well as Marco and Kai who I believe are interested in this topic.

Rik <rik5>
Project Administrator
Thu 10 Oct 2019 04:57:19 PM UTC, original submission:  

I appended a test patch and a zip file. The zip file contains the script test_spqr.m with test cases.  Two test cases use SuiteSparse test matrices in Matrix Market Format. So I added the test matrices and MatrixMarketUtils to the zip file. I hope the patch will work on your system. 

The current status of the patch is
1) QR factorization can be calculated for sparse (real|complex) matrices.
2) least square and min2norm solutions can be calculated for sparse (real|complex) matrices.

Remarks
1) I will improve code formatting. I didn't apply the GNU coding standards consistently.
2) Probably the error handling and validation of parameters can be improved.
3) The solve methods for least square and min2norm are moved from class sparse_qr_rep to class sparse_qr. I can change this if necessary.
4) Macros HAVE_SUITESPARSE_SPQR_H and HAVE_SPQR are defined in oct-sparse.h. I don't know how to put them into config.h.
5) I didn't do any tests on Windows or Mac.
6) SPQR stores the results of the QR factorization either in a collection of data structures (H, R, HPinv, HTau) or in one "big" structure of type SuiteSparseQR_factorization. Unfortunately, there is no user-callable SPQR function to extract H, R, HPinv and HTau from a SuiteSparseQR_factorization-object. At least I didn't find one.

The implementation of the R-,C-,Q- and V-methods in class sparse_qr_rep is based on H, R, HPinv and HTau (calculated in the constructor). For the implementation of the solve methods the function SuiteSparseQR_min2norm is very convenient. This function uses the SuiteSparseQR_factorization and does not reuse H, R, HPinv and HTau. So for a matrix A I do one QR factorization in order to get R, Q etc. and a further factorization to solve Ax = b.

In the old implementation of sparse_qr_rep there was only one QR factorization in the constructor. The solve methods and the calculation of R, C, Q and V based on this single factorization. Therefore I think that the old approach is more efficient than my implementation.

I am thankful for tips or comments how to proceed.

Simon

Rik <rik5>
Project Administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #48691:  bug57033_spqr_rev6.patch added by mmuetzel (72KiB - application/octet-stream)
file #48596:  bug57033_spqr_rev5.patch added by sim79 (72KiB - text/x-patch)
file #48481:  bug57033_spqr_rev4.patch added by mmuetzel (66KiB - application/octet-stream)
file #48196:  bug57033_spqr_rev3.patch added by sim79 (77KiB - text/x-patch)
file #48139:  bug57033_spqr_rev2.patch added by sim79 (54KiB - text/x-patch)
file #48140:  bug57033_tests_rev2.log added by sim79 (4KiB - text/x-log)
file #47957:  bug57033_failing_tests.log added by mmuetzel (124KiB - application/octet-stream)
file #47958:  bug57033_spqr_rev1.patch added by mmuetzel (63KiB - application/octet-stream)
file #47661:  bug57033_spqr_configure.patch added by mmuetzel (4KiB - application/octet-stream)
file #47640:  test_spqr_20191006.patch added by rik5 (61KiB - text/x-patch)
file #47641:  test_spqr_20191006.zip added by rik5 (1MiB - application/zip)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by sim79 (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Submitted the item)
  • -email is unavailable- added by rik5 (New patch to reimplement sparse operations with SPQR)
  • -email is unavailable- added by rik5 (New patch to reimplement sparse operations with SPQR)
  • -email is unavailable- added by rik5 (New patch to reimplement sparse operations with SPQR)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-03-28 mmuetzel Attached File- => Added bug57033_spqr_rev6.patch, #48691
    2020-03-13 sim79 Attached File- => Added bug57033_spqr_rev5.patch, #48596
    2020-02-23 mmuetzel Attached File- => Added bug57033_spqr_rev4.patch, #48481
    2020-01-10 sim79 Attached File- => Added bug57033_spqr_rev3.patch, #48196
    2019-12-30 sim79 Attached File- => Added bug57033_spqr_rev2.patch, #48139
        Attached File- => Added bug57033_tests_rev2.log, #48140
    2019-12-01 mmuetzel Attached File- => Added bug57033_failing_tests.log, #47957
        Attached File- => Added bug57033_spqr_rev1.patch, #47958
    2019-10-13 mmuetzel Dependencies- => bugs #51950 is dependent
    2019-10-13 mmuetzel Attached File- => Added bug57033_spqr_configure.patch, #47661
    2019-10-10 rik5 Attached File- => Added test_spqr_20191006.patch, #47640
        Attached File- => Added test_spqr_20191006.zip, #47641
        Carbon-Copy- => Added marco.caliari@univr.it
        Carbon-Copy- => Added k.ohlhus@gmail.com
        Carbon-Copy- => Added simon.hau79@gmail.com

    Back to the top


    Powered by Savane 3.5