Mercurial > octave
changeset 28739:3eb2fab71028
Allocate memory for Jacobian in __ode15__ more efficiently (bug #55905).
* __ode15__.cc (set_up): Initially allocate sparse matrix with 0 non-zero
elements.
* __ode15__.cc (jacsparse_impl): Increase allocated memory to match actual
number of non-zero elements.
* m4/acinclude.m4: Check if "SUNSparseMatrix_Reallocate" can be used.
| author | Bill Greene <w.h.greene@gmail.com> |
|---|---|
| date | Fri, 20 Sep 2019 08:48:45 +0200 |
| parents | e93ad0cfd9bd |
| children | fed3b09f2efc |
| files | libinterp/dldfcn/__ode15__.cc m4/acinclude.m4 |
| diffstat | 2 files changed, 20 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/__ode15__.cc Mon Sep 14 09:18:10 2020 -0700 +++ b/libinterp/dldfcn/__ode15__.cc Fri Sep 20 08:48:45 2019 +0200 @@ -395,13 +395,17 @@ if (m_havejacsparse) { # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) - // FIXME : one should not allocate space for a full Jacobian - // when using a sparse format. Consider allocating less space - // then possibly using SUNSparseMatrixReallocate to increase it. +# if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE) + // Initially allocate memory for 0 entries. We will reallocate when we + // get the Jacobian matrix from the user and know the actual number of + // entries. + m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT); +# else // FIXME: m_num*m_num might be larger than the largest // octave_f77_int_type (integer overflow). Consider using a save // calculation method. m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT); +# endif if (! m_sunJacMatrix) error ("Unable to create sparse Jacobian for Sundials"); @@ -482,6 +486,18 @@ else jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj); +# if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE) + octave_f77_int_type nnz = to_f77_int (jac.nnz ()); + if (nnz > SUNSparseMatrix_NNZ (Jac)) + { + // Allocate required memory for sparse Jacobian defined in user function + // This will always be required once since we set the number of + // non-zero elements to zero initially. + if (SUNSparseMatrix_Reallocate (Jac, nnz)) + error("Unable to allocate sufficient memory for IDA sparse matrix"); + } +# endif + SUNMatZero_Sparse (Jac); // We have to use "sunindextype *" here but still need to check that // conversion of each element to "octave_f77_int_type" is save.
--- a/m4/acinclude.m4 Mon Sep 14 09:18:10 2020 -0700 +++ b/m4/acinclude.m4 Fri Sep 20 08:48:45 2019 +0200 @@ -1937,7 +1937,7 @@ ac_octave_save_LIBS=$LIBS LIBS="$SUNDIALS_IDA_LIBS $SUNDIALS_NVECSERIAL_LIBS $LIBS" dnl Current API functions present in SUNDIALS version 4 - AC_CHECK_FUNCS([IDASetJacFn IDASetLinearSolver SUNLinSol_Dense]) + AC_CHECK_FUNCS([IDASetJacFn IDASetLinearSolver SUNLinSol_Dense SUNSparseMatrix_Reallocate]) dnl FIXME: The purpose of the following tests is to detect the deprecated dnl API from SUNDIALS version 3, which should only be used if the current dnl API tests above failed. For now, always test for ida_direct.h.