Mercurial > octave

changeset 32287:dde186dc7b4c stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix getting number of non-zero elements in SuiteSparse matrices. Based on changes by @DrTimothyAldenDavis on GitHub: https://github.com/DrTimothyAldenDavis/SuiteSparse/issues/221#issuecomment-1710586946 * liboctave/numeric/sparse-qr.cc (rcs2ros, ccs2cos): Add cholmod_common as input argument. Use "cholmod_l_nnz" to get number of non-zero matrix elements. (sparse_qr<T>::sparse_qr_rep::V): Pass cholmod_common to conversion function. Get correct number of non-zero matrix elements. (sparse_qr<T>::sparse_qr_rep::R): Use "cholmod_l_nnz" to get number of non-zero matrix elements. Get correct number of non-zero matrix elements. (sparse_qr<T>::sparse_qr_rep::tall_solve): Use "cholmod_l_nnz" to get number of non-zero matrix elements. (sparse_qr<T>::min2norm_solve): Pass cholmod_common to conversion function. * liboctave/numeric/sparse-qr.cc (parse_chol<chol_type>::L): Use "cholmod_[l_]nnz" to get number of non-zero matrix elements. * liboctave/array/CSparse.cc, dSparse.cc (fsolve): Use "cholmod_[l_]nnz" to get number of non-zero matrix elements.
author Markus Mützel <markus.muetzel@gmx.de>
date Thu, 07 Sep 2023 21:29:43 +0200
parents 6fa106924c19
children 197dc8992606 b262966c853c
files liboctave/array/CSparse.cc liboctave/array/dSparse.cc liboctave/numeric/sparse-chol.cc liboctave/numeric/sparse-qr.cc
diffstat 4 files changed, 50 insertions(+), 27 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/array/CSparse.cc	Mon Aug 28 11:44:31 2023 +0200
+++ b/liboctave/array/CSparse.cc	Thu Sep 07 21:29:43 2023 +0200
@@ -6042,12 +6042,14 @@
               retval = SparseComplexMatrix
                        (static_cast<octave_idx_type> (X->nrow),
                         static_cast<octave_idx_type> (X->ncol),
-                        static_cast<octave_idx_type> (X->nzmax));
+                        octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
+                                                       (X, cm)));
               for (octave_idx_type j = 0;
                    j <= static_cast<octave_idx_type> (X->ncol); j++)
                 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
               for (octave_idx_type j = 0;
-                   j < static_cast<octave_idx_type> (X->nzmax); j++)
+                   j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
+                   j++)
                 {
                   retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
                   retval.xdata (j) = static_cast<Complex *>(X->x)[j];
@@ -6550,12 +6552,14 @@
               retval = SparseComplexMatrix
                        (static_cast<octave_idx_type> (X->nrow),
                         static_cast<octave_idx_type> (X->ncol),
-                        static_cast<octave_idx_type> (X->nzmax));
+                        octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
+                                                       (X, cm)));
               for (octave_idx_type j = 0;
                    j <= static_cast<octave_idx_type> (X->ncol); j++)
                 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
               for (octave_idx_type j = 0;
-                   j < static_cast<octave_idx_type> (X->nzmax); j++)
+                   j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
+                   j++)
                 {
                   retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
                   retval.xdata (j) = static_cast<Complex *>(X->x)[j];
--- a/liboctave/array/dSparse.cc	Mon Aug 28 11:44:31 2023 +0200
+++ b/liboctave/array/dSparse.cc	Thu Sep 07 21:29:43 2023 +0200
@@ -6053,14 +6053,17 @@
 
               cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
 
-              retval = SparseMatrix (static_cast<octave_idx_type> (X->nrow),
-                                     static_cast<octave_idx_type> (X->ncol),
-                                     static_cast<octave_idx_type> (X->nzmax));
+              retval = SparseMatrix
+                       (static_cast<octave_idx_type> (X->nrow),
+                        static_cast<octave_idx_type> (X->ncol),
+                        octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
+                                                       (X, cm)));
               for (octave_idx_type j = 0;
                    j <= static_cast<octave_idx_type> (X->ncol); j++)
                 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
               for (octave_idx_type j = 0;
-                   j < static_cast<octave_idx_type> (X->nzmax); j++)
+                   j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
+                   j++)
                 {
                   retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
                   retval.xdata (j) = static_cast<double *>(X->x)[j];
@@ -6548,12 +6551,14 @@
               retval = SparseComplexMatrix
                        (static_cast<octave_idx_type> (X->nrow),
                         static_cast<octave_idx_type> (X->ncol),
-                        static_cast<octave_idx_type> (X->nzmax));
+                        octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
+                                                       (X, cm)));
               for (octave_idx_type j = 0;
                    j <= static_cast<octave_idx_type> (X->ncol); j++)
                 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
               for (octave_idx_type j = 0;
-                   j < static_cast<octave_idx_type> (X->nzmax); j++)
+                   j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
+                   j++)
                 {
                   retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
                   retval.xdata (j) = static_cast<Complex *>(X->x)[j];
--- a/liboctave/numeric/sparse-chol.cc	Mon Aug 28 11:44:31 2023 +0200
+++ b/liboctave/numeric/sparse-chol.cc	Thu Sep 07 21:29:43 2023 +0200
@@ -97,6 +97,13 @@
   {
     return m_L;
   }
+
+  cholmod_common *cc (void) const
+  {
+    cholmod_common *ret = const_cast<cholmod_common *> (&m_common);
+    return ret;
+  }
+
 #endif
 
   octave_idx_type P (void) const
@@ -421,7 +428,8 @@
   cholmod_sparse *m = m_rep->L ();
 
   octave_idx_type nc = m->ncol;
-  octave_idx_type nnz = m->nzmax;
+  octave_idx_type nnz
+    = from_suitesparse_long (CHOLMOD_NAME(nnz) (m, m_rep->cc ()));
 
   chol_type ret (m->nrow, nc, nnz);
 
--- a/liboctave/numeric/sparse-qr.cc	Mon Aug 28 11:44:31 2023 +0200
+++ b/liboctave/numeric/sparse-qr.cc	Thu Sep 07 21:29:43 2023 +0200
@@ -371,11 +371,13 @@
 // Convert real sparse cholmod matrix to real sparse octave matrix.
 // Returns a "shallow" copy of y.
 static SparseMatrix
-rcs2ros (const cholmod_sparse *y)
+rcs2ros (cholmod_sparse *y, const cholmod_common *cc1)
 {
   octave_idx_type nrow = from_size_t (y->nrow);
   octave_idx_type ncol = from_size_t (y->ncol);
-  octave_idx_type nz = from_size_t (y->nzmax);
+  octave_idx_type nz
+    = from_suitesparse_long
+      (cholmod_l_nnz (y, const_cast<cholmod_common *> (cc1)));
   SparseMatrix ret (nrow, ncol, nz);
 
   SuiteSparse_long *y_p = reinterpret_cast<SuiteSparse_long *> (y->p);
@@ -396,11 +398,11 @@
 // Convert complex sparse cholmod matrix to complex sparse octave matrix.
 // Returns a "deep" copy of a.
 static SparseComplexMatrix
-ccs2cos (const cholmod_sparse *a)
+ccs2cos (cholmod_sparse *a, cholmod_common *cc1)
 {
   octave_idx_type nrow = from_size_t (a->nrow);
   octave_idx_type ncol = from_size_t (a->ncol);
-  octave_idx_type nz = from_size_t (a->nzmax);
+  octave_idx_type nz = from_suitesparse_long (cholmod_l_nnz (a, cc1));
   SparseComplexMatrix ret (nrow, ncol, nz);
 
   SuiteSparse_long *a_p = reinterpret_cast<SuiteSparse_long *> (a->p);
@@ -595,7 +597,7 @@
 {
 #if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
 
-  return rcs2ros (m_H);
+  return rcs2ros (m_H, &m_cc);
 
 #elif defined (HAVE_CXSPARSE)
 
@@ -609,7 +611,7 @@
   CXSPARSE_DNAME (_spfree) (D);
 
   octave_idx_type nc = N->L->n;
-  octave_idx_type nz = N->L->nzmax;
+  octave_idx_type nz = N->L->p[nc];
   SparseMatrix ret (N->L->m, nc, nz);
 
   for (octave_idx_type j = 0; j < nc+1; j++)
@@ -638,7 +640,9 @@
 
   octave_idx_type nr = static_cast<octave_idx_type> (m_R->nrow);
   octave_idx_type nc = static_cast<octave_idx_type> (m_R->ncol);
-  octave_idx_type nz = static_cast<octave_idx_type> (m_R->nzmax);
+  octave_idx_type nz
+    = from_suitesparse_long
+      (cholmod_l_nnz (m_R, const_cast<cholmod_common *> (&m_cc)));
 
   // FIXME: Does this work if econ = true?
   SparseMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz);
@@ -668,7 +672,7 @@
   CXSPARSE_DNAME (_spfree) (D);
 
   octave_idx_type nc = N->U->n;
-  octave_idx_type nz = N->U->nzmax;
+  octave_idx_type nz = N->U->p[nc];
 
   SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
 
@@ -942,7 +946,7 @@
       for (octave_idx_type i = 0; i < ncols+1; i++)
         R2_p[i] = suitesparse_long_to_suitesparse_integer (R_p[i]);
       R2.p = R2_p;
-      octave_idx_type nnz = m_R->nzmax;
+      octave_idx_type nnz = from_suitesparse_long (cholmod_l_nnz (m_R, &m_cc));
       R2_i = new suitesparse_integer[nnz];
       SuiteSparse_long *R_i = reinterpret_cast<SuiteSparse_long *> (m_R->i);
       for (octave_idx_type i = 0; i < nnz; i++)
@@ -1542,7 +1546,7 @@
   CXSPARSE_ZNAME (_spfree) (D);
 
   octave_idx_type nc = N->L->n;
-  octave_idx_type nz = N->L->nzmax;
+  octave_idx_type nz = N->L->p[nc];
   SparseComplexMatrix ret (N->L->m, nc, nz);
 
   for (octave_idx_type j = 0; j < nc+1; j++)
@@ -1571,7 +1575,9 @@
 
   octave_idx_type nr = from_size_t (m_R->nrow);
   octave_idx_type nc = from_size_t (m_R->ncol);
-  octave_idx_type nz = from_size_t (m_R->nzmax);
+  octave_idx_type nz
+    = from_suitesparse_long
+      (cholmod_l_nnz (m_R, const_cast<cholmod_common *> (&m_cc)));
 
   // FIXME: Does this work if econ = true?
   SparseComplexMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz);
@@ -1601,7 +1607,7 @@
   CXSPARSE_ZNAME (_spfree) (D);
 
   octave_idx_type nc = N->U->n;
-  octave_idx_type nz = N->U->nzmax;
+  octave_idx_type nz = N->U->p[nc];
 
   SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows),
                            nc, nz);
@@ -2822,7 +2828,7 @@
       delete [] reinterpret_cast<SuiteSparse_long *> (B.i);
     }
 
-  x = rcs2ros (X);
+  x = rcs2ros (X, &cc);
   cholmod_l_finish (&cc);
   info = 0;
 
@@ -2897,7 +2903,7 @@
     }
   cholmod_l_finish (&cc);
 
-  SparseComplexMatrix ret = ccs2cos(X);
+  SparseComplexMatrix ret = ccs2cos (X, &cc);
 
   info = 0;
 
@@ -3022,7 +3028,7 @@
 
   info = 0;
 
-  return ccs2cos (X);
+  return ccs2cos (X, &cc);
 
 }
 
@@ -3047,7 +3053,7 @@
   X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, &A, B, &cc);
   spqr_error_handler (&cc);
 
-  SparseComplexMatrix ret = ccs2cos(X);
+  SparseComplexMatrix ret = ccs2cos (X, &cc);
 
   if (sizeof (octave_idx_type) != sizeof (SuiteSparse_long))
     {