[Octave-bug-tracker] [bug #57033] Replace CXSPARSE with SPQR

[Markus Mützel]

Follow-up Comment #41, bug #57033 (project octave):

I tried to implement something like outlined in comment #40.

The member m_E is passed as an argument to functions from SPQR and from
CXSparse.
For SPQR:
https://hg.savannah.gnu.org/hgweb/octave/file/d8cb0f06ffb1/liboctave/numeric/sparse-qr.cc#l263

      SuiteSparseQR<double> (order, SPQR_DEFAULT_TOL, A.nrow,
                             const_cast<cholmod_sparse*> (&A), &m_R, &m_E,
&m_H,
                             &m_HPinv, &m_Htau, &m_cc);

With the signature:

template <typename Entry> SuiteSparse_long SuiteSparseQR
(
    // inputs, not modified
    int ordering,           // all, except 3:given treated as 0:fixed
    double tol,             // only accept singletons above tol
    SuiteSparse_long econ,  // number of rows of C and R to return
    cholmod_sparse *A,      // m-by-n sparse matrix
    // outputs
    cholmod_sparse **R,     // the R factor
    SuiteSparse_long **E,   // permutation of 0:n-1, NULL if identity
    cholmod_sparse **H,     // the Householder vectors (m-by-nh)
    SuiteSparse_long **HPinv,// size m; row permutation for H
    cholmod_dense **HTau,   // size nh, Householder coefficients
    cholmod_common *cc      // workspace and parameters
) ;

(Some `static_cast`s seem to be missing here in our code.)

For CXSparse:
https://hg.savannah.gnu.org/hgweb/octave/file/d8cb0f06ffb1/liboctave/numeric/sparse-qr.cc#l700

          CXSPARSE_DNAME (_ipvec)
            (m_E, &(static_cast<double*> (QTB->x)[j * b_nr]), &x_vec[j *
ncols],
             ncols);

With the signature (for systems with OCTAVE_ENABLE_64 *not* defined):

int cs_di_ipvec (const int *p, const double *b, double *x, int n) ;


So `m_E` must be of type `SuiteSparse_long*` (i.e. convertible to
`ptrdiff_t*`) for the one function and `int*` for the other. On 64bit systems,
those types are incompatible.
The only solution I can think of at the moment is copying the content of m_E
to a new array of the correct type.
Is there a better solution?

If OCTAVE_ENABLE_64 *is* defined, the signature of the used function is:

cs_long_t cs_dl_ipvec (const cs_long_t *p, const double *b, double *x,
cs_long_t n) ;

where `cs_long_t` is `SuiteSparse_long`.
So the types coincidently match.


Also, thinking about it, I'm no longer sure if our definition of
`octave::suitesparse_integer` is correct. At least, if the intent was to
follow CHOLMOD's BLAS_INT (which looks eerily similar), the conditions
shouldn't be OCTAVE_ENABLE_64, but rather whether OCTAVE_SIZEOF_F77_INT_TYPE
is 4 or 8...

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?57033>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/