# HG changeset patch # User Sarrah Bastawala <84874044+sarrah-basta@users.noreply.github.com> # Date 1662956177 -19800 # Mon Sep 12 09:46:17 2022 +0530 # Node ID f525ca4baa98b75af7e419b9f01cb1a709ccab7d # Parent 44cf6bbeeca91d028c74a3737b3916d5f23ba0a2 [PATCH] add support for sparse computations in ode15{i,s} using Octave classes, making the dependency on SUNDIALS interface to KLU optional ( GSoC '22) * libinterp/dldfcn/oct-sundials.h : add declarations of custom implementations to SUNDIALS NVECTOR, SUNMATRIX and SUNLINEARSOLVER APIs using Octave's classes. * libinterp/dldfcn/oct-sundials.cc: add definitions of exported functions declared in corresponding header file * libinterp/dldfcn/__ode15__.cc : add motifications and optimizations to make ode15{i,s} solvers default to using Octave's custom implementations and solver in IDA for sparse computations if KLU is not found.binterp * libinterp/dldfcn/config-module.awk : add build dependencies for __ode15__ oct-file. * configure.ac : modify warnings to let user know ode15{i,s} support sparse computations even if KLU is not found. * scripts/ode/ode15i.m : modify tests using sparse Jacobians to run even without KLU * scripts/ode/ode15s.m : modify tests using sparse Jacobians to run even without KLU --- configure.ac | 4 +- libinterp/dldfcn/__ode15__.cc | 295 ++++++-- libinterp/dldfcn/config-module.awk | 10 +- libinterp/dldfcn/oct-sundials.cc | 1134 ++++++++++++++++++++++++++++ libinterp/dldfcn/oct-sundials.h | 244 ++++++ scripts/ode/ode15i.m | 4 +- scripts/ode/ode15s.m | 10 +- 7 files changed, 1618 insertions(+), 83 deletions(-) create mode 100644 libinterp/dldfcn/oct-sundials.cc create mode 100644 libinterp/dldfcn/oct-sundials.h diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -2189,7 +2189,9 @@ ### Check for KLU library and header. OCTAVE_CHECK_LIB(klu, KLU, - [KLU library not found. This will result in some lack of functionality for sparse matrices.], + [KLU library not found. This will result in some lack of functionality for sparse matrices. + Sparse computations in ode15{i,s} solvers will still be supported but will suffer a hit in + time taken.], [klu.h suitesparse/klu.h klu/klu.h], [klu_solve], [], [don't use KLU library, disable some sparse matrix functionality]) diff --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc --- a/libinterp/dldfcn/__ode15__.cc +++ b/libinterp/dldfcn/__ode15__.cc @@ -23,16 +23,18 @@ // //////////////////////////////////////////////////////////////////////// +/* include guards so that files included + in dependencies are not included again */ +#ifndef __CONFIG_H_INCLUDE_GUARD +#define __CONFIG_H_INCLUDE_GUARD #if defined (HAVE_CONFIG_H) # include "config.h" #endif +#endif -#include "dColVector.h" -#include "dMatrix.h" -#include "dSparse.h" +#include #include "f77-fcn.h" #include "lo-utils.h" - #include "Cell.h" #include "defun-dld.h" #include "error.h" @@ -43,11 +45,14 @@ #include "pager.h" #include "parse.h" + #if defined (HAVE_SUNDIALS) -# if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H) -# include -# endif +# include "oct-sundials.h" + + #if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H) + #include + #endif # if defined (HAVE_IDA_IDA_H) # include @@ -78,7 +83,8 @@ # include # endif # include -# endif +# endif + #endif @@ -121,6 +127,25 @@ # endif static inline realtype * + nv_data_c (N_Vector& v) + { +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Disable warning from GCC about old-style casts in Sundials + // macro expansions. Do this in a function so that this + // diagnostic may still be enabled for the rest of the file. +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wold-style-cast" +# endif + + return NV_DATA_C (v); + +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Restore prevailing warning state for remainder of the file. +# pragma GCC diagnostic pop +# endif + } + + static inline realtype * nv_data_s (N_Vector& v) { # if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) @@ -139,6 +164,45 @@ # endif } + static inline _N_VectorContent_Serial* + nv_content_s (N_Vector& v) + { +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Disable warning from GCC about old-style casts in Sundials + // macro expansions. Do this in a function so that this + // diagnostic may still be enabled for the rest of the file. +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wold-style-cast" +# endif + + return NV_CONTENT_S (v); + +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Restore prevailing warning state for remainder of the file. +# pragma GCC diagnostic pop +# endif + } + + static inline ColumnVector * + nv_content_c (N_Vector& v) + { +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Disable warning from GCC about old-style casts in Sundials + // macro expansions. Do this in a function so that this + // diagnostic may still be enabled for the rest of the file. +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wold-style-cast" +# endif + + return NV_CONTENT_C (v); + +# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) + // Restore prevailing warning state for remainder of the file. +# pragma GCC diagnostic pop +# endif + } + + class IDA { public: @@ -257,9 +321,11 @@ static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n); # if defined (HAVE_SUNDIALS_SUNCONTEXT) - N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n); + N_Vector ColToNVec_Serial (const ColumnVector& data, octave_f77_int_type n); + N_Vector ColToNVec_Octave (const ColumnVector& data, octave_f77_int_type n); # else - static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n); + static N_Vector ColToNVec_Serial (const ColumnVector& data, octave_f77_int_type n); + static N_Vector ColToNVec_Octave (const ColumnVector& data, octave_f77_int_type n); # endif void @@ -292,7 +358,6 @@ jacdense_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ); -# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) static int jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp, N_Vector, SUNMatrix Jac, void *user_data, N_Vector, @@ -306,7 +371,7 @@ void jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, SUNMatrix& Jac); -# endif + void set_maxstep (realtype maxstep); @@ -399,23 +464,24 @@ ColumnVector res = (*m_fcn) (y, yp, t, m_ida_fcn); - realtype *puntrr = nv_data_s (rr); - - for (octave_idx_type i = 0; i < m_num; i++) - puntrr[i] = res(i); - } + realtype *puntrr; + if(N_VGetVectorID(rr) == SUNDIALS_NVEC_CUSTOM) + { + ColumnVector *rest = const_cast (nv_content_c(rr)); + *rest = res; + } + else + { + puntrr = nv_data_s (rr); + for (octave_idx_type i = 0; i < m_num; i++) + puntrr[i] = res(i); + } -# if defined (HAVE_SUNDIALS_SUNCONTEXT) -# define OCTAVE_SUNCONTEXT , m_sunContext -# else -# define OCTAVE_SUNCONTEXT -# endif + } void IDA::set_up (const ColumnVector& y) { - N_Vector yy = ColToNVec (y, m_num); - if (m_havejacsparse) { # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) @@ -423,38 +489,39 @@ // 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. + std::cout<<"Setting up using KLU sparse solver \n"; + N_Vector yy = ColToNVec_Serial (y, m_num); m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT OCTAVE_SUNCONTEXT); -# else - octave_f77_int_type max_elems; - if (math::int_multiply_overflow (m_num, m_num, &max_elems)) - error ("Unable to allocate memory for sparse Jacobian"); - - m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT - OCTAVE_SUNCONTEXT); -# endif +# endif if (! m_sunJacMatrix) error ("Unable to create sparse Jacobian for Sundials"); m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix OCTAVE_SUNCONTEXT); if (! m_sunLinearSolver) + error ("Unable to create KLU sparse solver"); +# else + std::cout<<"Setting up using Octave sparse solver \n"; + N_Vector yy = ColToNVec_Octave (y, m_num); + m_sunJacMatrix = OCTSparseMatrix (m_num, m_num, 0 OCTAVE_SUNCONTEXT); + if (! m_sunJacMatrix) + error ("Unable to create sparse Jacobian for Sundials"); + + m_sunLinearSolver = OCTLinSol_Gen (yy, m_sunJacMatrix + OCTAVE_SUNCONTEXT); + if (! m_sunLinearSolver) error ("Unable to create KLU sparse solver"); - +# endif if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix)) error ("Unable to set sparse linear solver"); IDASetJacFn (m_mem, IDA::jacsparse); -# else - error ("SUNDIALS SUNLINSOL KLU was unavailable or disabled when " - "Octave was built"); - -# endif - } else { + N_Vector yy = ColToNVec_Serial (y, m_num); m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT); if (! m_sunJacMatrix) @@ -479,12 +546,9 @@ N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ) { - octave_f77_int_type Neq = NV_LENGTH_S (yy); - + octave_f77_int_type Neq = N_VGetLength(yy); ColumnVector y = NVecToCol (yy, Neq); - ColumnVector yp = NVecToCol (yyp, Neq); - Matrix jac; if (m_havejacfcn) @@ -498,16 +562,14 @@ SUNDenseMatrix_Data (JJ)); } -# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) void IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, SUNMatrix& Jac) { +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) ColumnVector y = NVecToCol (yy, m_num); - ColumnVector yp = NVecToCol (yyp, m_num); - SparseMatrix jac; if (m_havejacfcn) @@ -542,26 +604,57 @@ rowvals[i] = to_f77_int (jac.ridx (i)); d[i] = jac.data (i); } +# else + ColumnVector *y = const_cast (nv_content_c(yy)); + ColumnVector *yp = const_cast (nv_content_c(yyp)); + SparseMatrix jac; + + if (m_havejacfcn) + jac = (*m_jacspfcn) (*y, *yp, t, cj, m_ida_jac); + else + jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj); + + SUNMatZero(Jac); + // We have to use "sunindextype *" here but still need to check that + // conversion of each element to "octave_f77_int_type" is save. + + SparseMatrix *jnew = new SparseMatrix(); + *jnew = jac; + SparseMatrix *content = const_cast (jnew); + Jac->content = content; +# endif } -# endif ColumnVector IDA::NVecToCol (N_Vector& v, octave_f77_int_type n) { ColumnVector data (n); - realtype *punt = nv_data_s (v); + if(N_VGetVectorID(v) == SUNDIALS_NVEC_CUSTOM) + { + ColumnVector *punt = const_cast (nv_content_c(v)); + data = *punt; + } + else + { + realtype *punt = nv_data_s (v); - for (octave_f77_int_type i = 0; i < n; i++) - data(i) = punt[i]; - + for (octave_f77_int_type i = 0; i < n; i++) + data(i) = punt[i]; + } return data; } N_Vector - IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n) + IDA::ColToNVec_Octave (const ColumnVector& data, octave_f77_int_type n) + { + N_Vector v = N_VMake_Octave (data OCTAVE_SUNCONTEXT); + return v; + } + + N_Vector + IDA::ColToNVec_Serial (const ColumnVector& data, octave_f77_int_type n) { N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT); - realtype *punt = nv_data_s (v); for (octave_f77_int_type i = 0; i < n; i++) @@ -590,13 +683,23 @@ # else m_mem = IDACreate (); # endif - - N_Vector yy = ColToNVec (m_y0, m_num); - - N_Vector yyp = ColToNVec (m_yp0, m_num); + N_Vector yy,yyp; + if (m_havejacsparse){ +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + yy = ColToNVec_Serial (m_y0, m_num); + yyp = ColToNVec_Serial (m_yp0, m_num); +# else + yy = ColToNVec_Octave (m_y0, m_num); + yyp = ColToNVec_Octave (m_yp0, m_num); +# endif + } + else + { + yy = ColToNVec_Serial (m_y0, m_num); + yyp = ColToNVec_Serial (m_yp0, m_num); + } IDA::set_userdata (); - if (IDAInit (m_mem, IDA::resfun, m_t0, yy, yyp) != 0) error ("IDA not initialized"); } @@ -604,12 +707,23 @@ void IDA::set_tolerance (ColumnVector& abstol, realtype reltol) { - N_Vector abs_tol = ColToNVec (abstol, m_num); + N_Vector abs_tol; + if (m_havejacsparse){ +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + abs_tol = ColToNVec_Serial (abstol, m_num); +# else + abs_tol = ColToNVec_Octave (abstol, m_num); +# endif + } + else + { + abs_tol = ColToNVec_Serial (abstol, m_num); + } if (IDASVtolerances (m_mem, reltol, abs_tol) != 0) error ("IDA: Tolerance not set"); - - N_VDestroy_Serial (abs_tol); + + N_VDestroy (abs_tol); } void @@ -639,10 +753,21 @@ realtype tsol = tspan(0); realtype tend = tspan(numt-1); - - N_Vector yyp = ColToNVec (yp, m_num); - - N_Vector yy = ColToNVec (y, m_num); + N_Vector yy, yyp; + if (m_havejacsparse){ +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + yy = ColToNVec_Serial (y, m_num); + yyp = ColToNVec_Serial (yp, m_num); +# else + yy = ColToNVec_Octave (y, m_num); + yyp = ColToNVec_Octave (yp, m_num); +# endif + } + else + { + yy = ColToNVec_Serial (y, m_num); + yyp = ColToNVec_Serial (yp, m_num); + } // Initialize OutputFcn if (haveoutputfcn) @@ -654,7 +779,6 @@ status = IDA::event (event_fcn, te, ye, ie, tsol, y, "init", yp, oldval, oldisterminal, olddir, cont, temp, tsol, yold, num_event_args); - if (numt > 2) { // First output value @@ -669,10 +793,8 @@ for (octave_idx_type j = 1; j < numt && status == 0; j++) { // IDANORMAL already interpolates tspan(j) - if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0) error ("IDASolve failed"); - yout = NVecToCol (yy, m_num); ypout = NVecToCol (yyp, m_num); tout(j) = tsol; @@ -733,6 +855,7 @@ output.resize (cont + 1, m_num); // This may be not efficient tout.resize (cont + 1); tout(cont) = tsol; + yout = NVecToCol (yy, m_num); for (octave_idx_type i = 0; i < m_num; i++) @@ -752,12 +875,24 @@ if (status == 0) { // Interpolate in tend - N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + N_Vector dky; + if (m_havejacsparse){ + # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + # else + dky = N_VNew_Octave (m_num OCTAVE_SUNCONTEXT); + # endif + } + else + { + dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + } if (IDAGetDky (m_mem, tend, 0, dky) != 0) error ("IDA failed to interpolate y"); tout(cont) = tend; + yout = NVecToCol (dky, m_num); for (octave_idx_type i = 0; i < m_num; i++) @@ -777,7 +912,7 @@ yold, num_event_args); } - N_VDestroy_Serial (dky); + N_VDestroy(dky); } // Cleanup plotter @@ -917,9 +1052,21 @@ realtype h = 0, tcur = 0; bool status = 0; - N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); - - N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + N_Vector dky, dkyp; + if (m_havejacsparse){ +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); +# else + dky = N_VNew_Octave (m_num OCTAVE_SUNCONTEXT); + dkyp = N_VNew_Octave (m_num OCTAVE_SUNCONTEXT); +# endif + } + else + { + dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); + } ColumnVector yout (m_num); ColumnVector ypout (m_num); @@ -967,7 +1114,7 @@ tout(cont-1), yold, num_event_args); } - N_VDestroy_Serial (dky); + N_VDestroy(dky); return status; } @@ -1351,4 +1498,4 @@ %!assert (1) */ -OCTAVE_NAMESPACE_END +OCTAVE_NAMESPACE_END \ No newline at end of file diff --git a/libinterp/dldfcn/config-module.awk b/libinterp/dldfcn/config-module.awk --- a/libinterp/dldfcn/config-module.awk +++ b/libinterp/dldfcn/config-module.awk @@ -68,8 +68,16 @@ basename = files[i]; sub (/\.cc$/, "", basename); print ""; - printf ("%%canon_reldir%%_%s_la_SOURCES = %%reldir%%/%s\n", + if (files[i] == "__ode15__.cc") + { + printf ("%%canon_reldir%%_%s_la_SOURCES = %%reldir%%/__ode15__.cc %%reldir%%/oct-sundials.cc \n", + basename); + } + else + { + printf ("%%canon_reldir%%_%s_la_SOURCES = %%reldir%%/%s\n", basename, files[i]); + } if (cppflags[i]) { printf ("%%canon_reldir%%_%s_la_CPPFLAGS = $(libinterp_liboctinterp_la_CPPFLAGS) %s\n", diff --git a/libinterp/dldfcn/oct-sundials.cc b/libinterp/dldfcn/oct-sundials.cc new file mode 100644 --- /dev/null +++ b/libinterp/dldfcn/oct-sundials.cc @@ -0,0 +1,1134 @@ +# include "oct-sundials.h" + +#ifdef __cplusplus +extern "C" { +#endif + +using namespace octave; + +/* -------------------------------------------------------------------------- + * Utility functions + * -------------------------------------------------------------------------- */ +/* Private functions for special cases of vector operations */ +static void VCopy_Octave(N_Vector x, N_Vector z); + +/* -------------------------------------------------------------------------- + * Constructors + * -------------------------------------------------------------------------- */ + +/* Returns vector type ID. Used to identify vector implementation + from abstract N_Vector interface. */ +N_Vector_ID N_VGetVectorID_Octave([[maybe_unused]] N_Vector v) +{ + return SUNDIALS_NVEC_CUSTOM; +} + +/* Function to create a new empty octave vector */ + + N_Vector N_VNewEmpty_Octave(SP_ARG_SUNCONTEXT) +{ + N_Vector v; + + /* Create an empty vector object */ + v = NULL; + v = N_VNewEmpty(SP_OCTAVE_SUNCONTEXT); + if (v == NULL) + return (NULL); + + /* Attach operations */ + + /* constructors, destructors, and utility operations */ + v->ops->nvgetvectorid = N_VGetVectorID_Octave; + v->ops->nvclone = N_VClone_Octave; + v->ops->nvcloneempty = N_VCloneEmpty_Octave; + v->ops->nvdestroy = N_VDestroy_Octave; + v->ops->nvspace = N_VSpace_Octave; + v->ops->nvgetarraypointer = N_VGetArrayPointer_Octave; + v->ops->nvsetarraypointer = NULL; + v->ops->nvgetlength = N_VGetLength_Octave; + + /* standard vector operations */ + v->ops->nvlinearsum = N_VLinearSum_Octave; + v->ops->nvconst = N_VConst_Octave; + v->ops->nvprod = N_VProd_Octave; + v->ops->nvdiv = N_VDiv_Octave; + v->ops->nvscale = N_VScale_Octave; + v->ops->nvabs = N_VAbs_Octave; + v->ops->nvinv = N_VInv_Octave; + v->ops->nvaddconst = N_VAddConst_Octave; + v->ops->nvdotprod = N_VDotProd_Octave; + v->ops->nvmaxnorm = N_VMaxNorm_Octave; + v->ops->nvwrmsnormmask = N_VWrmsNormMask_Octave; + v->ops->nvwrmsnorm = N_VWrmsNorm_Octave; + v->ops->nvmin = N_VMin_Octave; + v->ops->nvwl2norm = NULL; + v->ops->nvl1norm = NULL; + v->ops->nvcompare = N_VCompare_Octave; + v->ops->nvinvtest = N_VInvTest_Octave; + v->ops->nvconstrmask = N_VConstrMask_Octave; + v->ops->nvminquotient = N_VMinQuotient_Octave; + + /* fused vector operations (optional) */ + v->ops->nvlinearcombination = NULL; + v->ops->nvscaleaddmulti = N_VScaleAddMulti_Octave; + v->ops->nvdotprodmulti = NULL; + + // /* vector array operations (optional) */ + v->ops->nvlinearsumvectorarray = N_VLinearSumVectorArray_Octave; + v->ops->nvscalevectorarray = N_VScaleVectorArray_Octave; + v->ops->nvconstvectorarray = NULL; + v->ops->nvwrmsnormvectorarray = NULL; + v->ops->nvwrmsnormmaskvectorarray = NULL; + v->ops->nvscaleaddmultivectorarray = NULL; + v->ops->nvlinearcombinationvectorarray = NULL; + + /* + * OPTIONAL operations with no default implementation. + */ + + /* local reduction operations */ + v->ops->nvdotprodlocal = NULL; + v->ops->nvmaxnormlocal = NULL; + v->ops->nvminlocal = NULL; + v->ops->nvl1normlocal = NULL; + v->ops->nvinvtestlocal = NULL; + v->ops->nvconstrmasklocal = NULL; + v->ops->nvminquotientlocal = NULL; + v->ops->nvwsqrsumlocal = N_VWSqrSumLocal_Octave; + v->ops->nvwsqrsummasklocal = N_VWSqrSumMaskLocal_Octave; + + /* debugging functions */ + v->ops->nvprint = N_VPrint_Octave; + v->ops->nvprintfile = NULL; + + return v; +} + +/* constructor to create a new nvector implementation + with a required length */ +N_Vector N_VNew_Octave(int length ARG_SUNCONTEXT) +{ + N_Vector v; + void *content; + # if defined (HAVE_SUNDIALS_SUNCONTEXT) + if (m_sunContext == NULL) + return (NULL); + #endif + + /* Create an empty vector object */ + v = NULL; + v = N_VNewEmpty_Octave(SP_OCTAVE_SUNCONTEXT); + if (v == NULL) + return (NULL); + + /* Create data */ + if (length > 0) { + + /* Allocate memory & Create content */ + content = NULL; + content = new ColumnVector (length); + if (content == NULL) + { + N_VDestroy(v); + return (NULL); + } + /* Attach content */ + v->content = content; + if(NV_DATA_C(v) == NULL) { N_VDestroy_Octave(v); return(NULL); } + +} +return (v); +} +/* Constructor to create a serial N_Vector + from existing ColumnVector component */ + +N_Vector N_VMake_Octave(const ColumnVector& cv ARG_SUNCONTEXT) +{ + N_Vector v; + void *content; + + /* Create an empty vector object */ + v = NULL; + v = N_VNewEmpty_Octave(SP_OCTAVE_SUNCONTEXT); + if (v == NULL) + return (NULL); + + /* Create content to point to allocated memory */ + content = NULL; + content = const_cast (&cv); + + if (content == NULL) + { + N_VDestroy(v); + return (NULL); + } + + /* Attach content */ + v->content = content; + return v; +} + +/* -------------------------------------------------------------------------- + * Vector get, set, and utility functions + * -------------------------------------------------------------------------- */ + +/* Function to return the global length of the vector.*/ +sunindextype N_VGetLength_Octave(N_Vector v) +{ + return NV_LENGTH_C(v); +} + +/* Function to print the vector vector to stdout */ +void N_VPrint_Octave(N_Vector x) +{ + ColumnVector *xv; + xv = static_cast NV_CONTENT_C(x); + + std::cout<<(*xv); + + return; +} + +/* -------------------------------------------------------------------------- + * Vector operations implementations + * -------------------------------------------------------------------------- */ + + +N_Vector N_VCloneEmpty_Octave(N_Vector w) +{ + N_Vector v; + + /* Check input */ + if (w == NULL) return(NULL); + + /* Create vector */ + v = NULL; + # if defined (HAVE_SUNDIALS_SUNCONTEXT) + v = N_VNewEmpty_Octave(w->sunctx); + #else + v = N_VNewEmpty_Octave(); + #endif + if (v == NULL) return(NULL); + + /* Attach operations */ + if (N_VCopyOps(w, v)) + { + N_VDestroy(v); + return(NULL); + } + + return(v); +} + +N_Vector N_VClone_Octave(N_Vector w) +{ + void* content; + + /* check inputs */ + if (w == NULL) return(NULL); + + /* Create an empty clone vector */ + N_Vector v = N_VCloneEmpty_Octave(w); + if (v == NULL) return(NULL); + /* Create content */ + content = NULL; + content = new ColumnVector (NV_LENGTH_C(w)); + + if (content == NULL) { N_VDestroy(v); return(NULL); } + + /* Attach content */ + v->content = content; + + return(v); +} + +void N_VDestroy_Octave(N_Vector v) +{ + if (v == NULL) return; + + /* free content */ + if (v->content != NULL) { + // our ColumnVector was created with 'new', so use 'delete' to destroy it. + delete NV_CONTENT_C(v); + v->content = NULL; + } + + /* free ops and vector */ + if (v->ops != NULL) { free(v->ops); v->ops = NULL; } + free(v); + v = NULL; + + return; +} + +void N_VSpace_Octave(N_Vector v, sunindextype *lrw, sunindextype *liw) +{ + *lrw = NV_LENGTH_C(v); + *liw = 1; + + return; +} + +realtype *N_VGetArrayPointer_Octave(N_Vector v) +{ + return((realtype *) NV_DATA_C(v)); +} + +/* + * Linear combination of two vectors: z = a*x + b*y + */ +void N_VLinearSum_Octave(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) +{ + + ColumnVector *xv, *yv, *zv; + xv = static_cast NV_CONTENT_C(x); + yv = static_cast NV_CONTENT_C(y); + zv = static_cast NV_CONTENT_C(z); + booleantype test; + + /* Case: a == b == 1.0 */ + + if ((a == ONE) && (b == ONE)) { + (*zv) = (*xv) + (*yv); + return; + } + + /* Cases: (1) a == 1.0, b = -1.0, (2) a == -1.0, b == 1.0 */ + + if ((test = ((a == ONE) && (b == -ONE))) || ((a == -ONE) && (b == ONE))) { + if(test == 0) + (*zv) = (*yv) - (*xv); + else + (*zv) = (*xv) - (*yv); + return; + } + + /* Cases: (1) a == 1.0, b == other or 0.0, (2) a == other or 0.0, b == 1.0 */ + /* if a or b is 0.0, then user should have called N_VScale */ + + if ((test = (a == ONE)) || (b == ONE)) { + if(test == 1) + (*zv) = (*xv) + b * (*yv); + else + (*zv) = a * (*xv) + (*yv); + return; + } + + /* Cases: (1) a == -1.0, b != 1.0, (2) a != 1.0, b == -1.0 */ + + if ((test = (a == -ONE)) || (b == -ONE)) { + if(test == 1) + (*zv) = b * (*yv) - (*xv); + else + (*zv) = a * (*xv) - (*yv); + return; + } + + /* Case: a == b */ + /* catches case both a and b are 0.0 - user should have called N_VConst */ + + if (a == b) { + (*zv) = a * ((*yv) + (*xv)); + return; + } + + /* Case: a == -b */ + + if (a == -b) { + (*zv) = b * ((*yv) - (*xv)); + return; + } + + /* Do all cases not handled above: + (1) a == other, b == 0.0 - user should have called N_VScale + (2) a == 0.0, b == other - user should have called N_VScale + (3) a,b == other, a !=b, a != -b */ + + (*zv) = a * (*xv) + b * (*yv); + return; +} + +/* + * Set all vector elements to a constant: z[i] = c + */ +void N_VConst_Octave(realtype c, N_Vector z) +{ + ColumnVector *zv; + zv = static_cast NV_CONTENT_C(z); + + zv->fill(c); + + return; +} + +/* + * Elementwise multiply vectors: z[i] = x[i]*y[i] + */ +void N_VProd_Octave(N_Vector x, N_Vector y, N_Vector z) +{ + + ColumnVector *xv, *zv, *yv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + yv = static_cast NV_CONTENT_C(y); + *zv = product((*xv),(*yv)); + + return; +} + +/* + * Elementwise divide vectors: z[i] = x[i]/y[i] + */ +void N_VDiv_Octave(N_Vector x, N_Vector y, N_Vector z) +{ + + ColumnVector *xv, *zv, *yv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + yv = static_cast NV_CONTENT_C(y); + *zv = quotient((*xv),(*yv)); + + return; +} + +/* + * Scale vector: z = c*x + */ +void N_VScale_Octave(realtype c, N_Vector x, N_Vector z) +{ + ColumnVector *xv = const_cast NV_CONTENT_C(x); + ColumnVector *zv = const_cast NV_CONTENT_C(z); + + // checks if both NVectors z and x + // are stored in the same place + if (zv == xv) { /* BLAS usage: scale x <- cx */ + (*xv) *= c; + // Not using this results in poor iterative + // algorthm performance + return; + } + if (c == ONE) { + VCopy_Octave(x, z); // substituting with (*zv) = (*xv) leads to wrong results. + } else if (c == -ONE) { + (*zv) -= (*xv); + } else { + (*zv) = c * (*xv); + } + + return; + +} + +/* + * Elementwise absolute value: z[i] = |x[i]| + */ +void N_VAbs_Octave(N_Vector x, N_Vector z) +{ + ColumnVector *xv,*zv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + + (*zv) = xv->abs(); + return; +} + +/* + * Elementwise inverse: z[i] = 1/x[i] + */ +void N_VInv_Octave(N_Vector x, N_Vector z) +{ + ColumnVector *xv,*zv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + + (*zv) = ONE / (*xv); + + return; +} + +/* + * Add constant: z = x + b + */ +void N_VAddConst_Octave(N_Vector x, realtype b, N_Vector z) +{ + ColumnVector *xv,*zv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + + (*zv) = (*xv) + b; + return; +} + +/* + * Scalar product of vectors x and y + */ +realtype N_VDotProd_Octave(N_Vector x, N_Vector y) +{ + ColumnVector *xv,*yv; + realtype sum; + xv = static_cast NV_CONTENT_C(x); + yv = static_cast NV_CONTENT_C(y); + sum = (*xv).transpose() * (*yv); + return(sum); +} + +/* + * Max norm (L infinity) of vector x + */ +realtype N_VMaxNorm_Octave(N_Vector x) +{ + ColumnVector *xv, abret; + realtype ret; + xv = static_cast NV_CONTENT_C(x); + abret = static_cast (xv->abs()); + ColumnVector *abp = &(abret); + ret = abp->max(); + + return(ret); +} + +/* + * Weighted RMS norm + */ +realtype N_VWrmsNorm_Octave(N_Vector x, N_Vector w) +{ + octave_value_list ov = ovl((N_VWSqrSumLocal_Octave(x, w)/(NV_LENGTH_C(x)))); + return((octave::Fsqrt(ov,1)(0)).double_value()); +} + +/* + * local weighted squared sum + */ +realtype N_VWSqrSumLocal_Octave(N_Vector x, N_Vector w) +{ + N_Vector prod = N_VClone_Octave(x); + ColumnVector *xv,*yv, *pv; + realtype sum; + xv = static_cast NV_CONTENT_C(x); + yv = static_cast NV_CONTENT_C(w); + pv = static_cast NV_CONTENT_C(prod); + + (*pv) = product((*xv),(*yv)); + octave_value_list ov = ovl((*pv)); + sum = (octave::Fsumsq(ov,1)(0)).double_value(); + return(sum); +} + +/* + * Masked weighted RMS norm + */ +realtype N_VWrmsNormMask_Octave(N_Vector x, N_Vector w, N_Vector id) +{ + octave_value_list ov = ovl((N_VWSqrSumMaskLocal_Octave(x, w, id)/(NV_LENGTH_C(x)))); + return((octave::Fsqrt(ov,1)(0)).double_value()); +} + +/* + * local weighted masked squared sum + */ +realtype N_VWSqrSumMaskLocal_Octave(N_Vector x, N_Vector w, N_Vector id) +{ + N_Vector prod = N_VClone_Octave(x); + N_Vector prod2 = N_VClone_Octave(x); + ColumnVector *xv,*yv, *mv, *pv, *pv2; + realtype sum; + xv = static_cast NV_CONTENT_C(x); + yv = static_cast NV_CONTENT_C(w); + mv = static_cast NV_CONTENT_C(id); //mask vector + pv = static_cast NV_CONTENT_C(prod); + pv2 = static_cast NV_CONTENT_C(prod2); + + const octave_value_list ov = octave_value_list({(*mv),(*xv),(*mv)}); + const octave_value_list ov2 = octave_value_list({(*mv),(*yv),(*mv)}); + (*pv) = (octave::Fmerge(ov2,1)(0)).column_vector_value(); + (*pv2) = (octave::Fmerge(ov2,1)(0)).column_vector_value(); + sum = N_VWSqrSumLocal_Octave(prod, prod2); + + return(sum); +} + +/* Minimum value in Vector */ +realtype N_VMin_Octave(N_Vector x) +{ + ColumnVector *xv; + realtype min; + xv = static_cast NV_CONTENT_C(x); + + min = xv->min(); + return(min); +} + +/* + * Elementwise z[i] = |x[i]| >= c ? 1 : 0 + */ +void N_VCompare_Octave(realtype c, N_Vector x, N_Vector z) +{ + ColumnVector *xv,*zv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + /* + * FIXME : The comparisons here shouldn't be low-level, + * but octave::Fge which should be used gives an interpreter error + */ + for (int i = 0; i < (xv)->numel(); i++) { + (*zv)(i) = ((xv->abs())(i) >= c) ? ONE : ZERO; + } + + return; +} + +/* + * Elementwise inverse with zero checking: z[i] = 1/x[i], x[i] != 0 + */ +booleantype N_VInvTest_Octave(N_Vector x, N_Vector z) +{ + N_Vector y = N_VClone_Octave(x); + + ColumnVector *xv, *zv, *yv; + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + yv = static_cast NV_CONTENT_C(y); + booleantype no_zero_found = SUNTRUE ; + (*zv) = ONE /(*xv); + if(xv->nnz() != xv->numel()) + no_zero_found = SUNFALSE; + if(no_zero_found==SUNTRUE) + (*zv) = ONE /(*xv); + else { + (*yv) = (*zv); + const octave_value_list ov = ovl((*zv)); + (*yv) = (octave::Fisfinite(ov,1)(0)).column_vector_value(); + const octave_value_list ov2 = octave_value_list({(*yv),(*zv),(*yv)}); + (*zv) = (octave::Fmerge(ov2,1)(0)).column_vector_value(); + } + + return no_zero_found; +} + +/* FIXME : Same problem as Compare */ +booleantype N_VConstrMask_Octave(N_Vector c, N_Vector x, N_Vector m) +{ + realtype temp; + booleantype test; + + ColumnVector *cv, *xv, *mv; + xv = static_cast NV_CONTENT_C(x); + cv = static_cast NV_CONTENT_C(c); + mv = static_cast NV_CONTENT_C(m); + + for (sunindextype i = 0; i < NV_LENGTH_C(x); i++) { + (*mv)(i) = ZERO; + + /* Continue if no constraints were set for the variable */ + if ((*cv)(i) == ZERO) + continue; + + /* Check if a set constraint has been violated */ + test = ((cv->abs())(i) > ONEPT5 && (*xv)(i)*(*cv)(i) <= ZERO) || + ((cv->abs())(i) > HALF && (*xv)(i)*(*cv)(i) < ZERO); + if (test) { + temp = (*mv)(i) = ONE; + } + } + + /* Return false if any constraint was violated */ + return (temp == ONE) ? SUNFALSE : SUNTRUE; +} + +/* + * Find minimum quotient: minq = min ( num[i]/denom[i]), denom[i] != 0. + */ +realtype N_VMinQuotient_Octave(N_Vector num, N_Vector denom) +{ + N_Vector z = N_VClone_Octave(num); + ColumnVector *xv, *zv, *yv; + realtype min; + min = BIG_REAL; + + xv = static_cast NV_CONTENT_C(num); + yv = static_cast NV_CONTENT_C(denom); + zv = static_cast NV_CONTENT_C(z); + *zv = quotient((*xv),(*yv)); + min = (zv)->min(); + + return(min); +} + +/* + * ----------------------------------------------------------------- + * fused vector operations + * ----------------------------------------------------------------- + */ + +int N_VLinearCombination_Octave(int nvec, realtype* c, N_Vector* X, N_Vector z) +{ + /* invalid number of vectors */ + if (nvec < 1) return(-1); + + /* should have called N_VScale */ + if (nvec == 1) { + N_VScale_Octave(c[0], X[0], z); + return(0); + } + + /* should have called N_VLinearSum */ + if (nvec == 2) { + N_VLinearSum_Octave(c[0], X[0], c[1], X[1], z); + return(0); + } + + ColumnVector *xv, *zv; + + xv = static_cast NV_CONTENT_C(X[0]); + zv = static_cast NV_CONTENT_C(z); + + /* + * z = sum{ c[i] * X[i] }, i = 0,...,nvec-1 + */ + + (*zv) = (*c)*(*xv); + return(0); +} + + +int N_VScaleAddMulti_Octave(int nvec, realtype* a, N_Vector x, N_Vector* Y, N_Vector* Z) +{ + /* invalid number of vectors */ + if (nvec < 1) return(-1); + + /* should have called N_VLinearSum */ + if (nvec == 1) { + N_VLinearSum_Octave(a[0], x, ONE, Y[0], Z[0]); + return(0); + } + + ColumnVector *xv, *zv, *yv; + + yv = static_cast NV_CONTENT_C(Y[0]); + zv = static_cast NV_CONTENT_C(Z[0]); + xv = static_cast NV_CONTENT_C(x); + + /* + * Z[i][j] = Y[i][j] + a[i] * x[j] + */ + (*zv) = (*yv) + (*a)*(*xv); + return(0); +} + +/* + * ----------------------------------------------------------------- + * vector array operations + * ----------------------------------------------------------------- + */ + +int N_VLinearSumVectorArray_Octave(int nvec, + realtype a, N_Vector* X, + realtype b, N_Vector* Y, + N_Vector* Z) +{ + + ColumnVector *xv, *zv, *yv; + int i; + + /* invalid number of vectors */ + if (nvec < 1) return(-1); + + /* should have called N_VLinearSum */ + if (nvec == 1) { + N_VLinearSum_Octave(a, X[0], b, Y[0], Z[0]); + return(0); + } + + for (i=0; i NV_CONTENT_C(Y[i]); + zv = static_cast NV_CONTENT_C(Z[i]); + xv = static_cast NV_CONTENT_C(X[i]); + (*zv) = a * (*xv) + b * (*yv); + } + + return(0); +} + + +int N_VScaleVectorArray_Octave(int nvec, realtype* c, N_Vector* X, N_Vector* Z) +{ + sunindextype i,j; + realtype *xd=NULL; + ColumnVector *xv, *zv; + + xv = static_cast NV_CONTENT_C(X[0]); + zv = static_cast NV_CONTENT_C(Z[0]); + + /* invalid number of vectors */ + if (nvec < 1) return(-1); + + /* should have called N_VScale */ + if (nvec == 1) { + N_VScale_Octave(c[0], X[0], Z[0]); + return(0); + } + + /* + * X[i] *= c[i] + */ + /* only checking if NVectors X and Z + are stored in the same address */ + if (xv == zv) { + // (*xv) *= (*c); //this should work but doesn't, in both N_VScale + // // instead results in convergence failures and wrong answer + for (i=0; inumel(); j++) { + xd[j] *= c[i]; + } + } + return(0); + } + + /* + * Z[i] = c[i] * X[i] + */ + (*zv) = (*c) * (*xv); + return(0); +} + +/* + * ----------------------------------------------------------------- + * private functions for special cases of vector operations + * ----------------------------------------------------------------- + */ + + +static void VCopy_Octave(N_Vector x, N_Vector z) +// used in N_VScale and if just do (*zv)=(*xv) results in wrong answer. +{ + ColumnVector *xv, *zv; + + xv = static_cast NV_CONTENT_C(x); + zv = static_cast NV_CONTENT_C(z); + + for (sunindextype i = 0; i < xv->numel(); i++) + (*zv)(i) = (*xv)(i); + + return; +} + +/* + * ----------------------------------------------------------------- + * Enable / Disable fused and vector array operations + * ----------------------------------------------------------------- + */ + +int N_VEnableFusedOps_Octave(N_Vector v, booleantype tf) +{ + /* check that vector is non-NULL */ + if (v == NULL) return(-1); + + /* check that ops structure is non-NULL */ + if (v->ops == NULL) return(-1); + + if (tf) { + /* enable all fused vector operations */ + v->ops->nvlinearcombination = NULL; + v->ops->nvscaleaddmulti = N_VScaleAddMulti_Octave; + v->ops->nvdotprodmulti = NULL; + /* enable all vector array operations */ + v->ops->nvlinearsumvectorarray = N_VLinearSumVectorArray_Octave; + v->ops->nvscalevectorarray = N_VScaleVectorArray_Octave; + v->ops->nvconstvectorarray = NULL; + v->ops->nvwrmsnormvectorarray = NULL; + v->ops->nvwrmsnormmaskvectorarray = NULL; + v->ops->nvscaleaddmultivectorarray = NULL; + v->ops->nvlinearcombinationvectorarray = NULL; + } else { + /* disable all fused vector operations */ + v->ops->nvlinearcombination = NULL; + v->ops->nvscaleaddmulti = NULL; + v->ops->nvdotprodmulti = NULL; + /* disable all vector array operations */ + v->ops->nvlinearsumvectorarray = NULL; + v->ops->nvscalevectorarray = NULL; + v->ops->nvconstvectorarray = NULL; + v->ops->nvwrmsnormvectorarray = NULL; + v->ops->nvwrmsnormmaskvectorarray = NULL; + v->ops->nvscaleaddmultivectorarray = NULL; + v->ops->nvlinearcombinationvectorarray = NULL; + } + + /* return success */ + return(0); +} + + +/* +* ----------------------------------------------------------------- +* exported functions for octave implementation of SUNMATRIX_SPARSE +* ----------------------------------------------------------------- +*/ + +/* Constructor function to create a new sparse matrix + based on dimensions */ + +SUNMatrix OCTSparseMatrix(sunindextype M, sunindextype N, + sunindextype NNZ ARG_SUNCONTEXT) +{ + SUNMatrix A; + void *content; + + /* return with NULL matrix on illegal input */ + if ((M <= 0) || (N <= 0) || (NNZ < 0)) + return (NULL); + + /* Create an empty matrix object */ + A = NULL; + A = SUNMatNewEmpty(SP_OCTAVE_SUNCONTEXT); + if (A == NULL) + return (NULL); + + // /* Attach required operations */ + A->ops->getid = OCTMatGetID_Sparse; + A->ops->clone = NULL; + A->ops->destroy = OCTMatDestroy_Sparse; + A->ops->zero = OCTMatZero_Sparse; + A->ops->copy = NULL; + A->ops->scaleadd = NULL; + A->ops->scaleaddi = NULL; + A->ops->matvec = NULL; + A->ops->space = NULL; + + /* Create content */ + content = NULL; + // creating a new SparseMatrix Octave data structure with rows and columns as specified + content = new SparseMatrix(M, N, NNZ); + if (content == NULL) + { + SUNMatDestroy(A); + return (NULL); + } + + /* Attach content */ + A->content = content; + /* Fill content */ + return (A); +} + +/* Functions to access the contents + of the sparse matrix structure */ + +sunindextype OCTSparseMatrix_Rows(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_ROWS_O(A); + else + return SUNMAT_ILL_INPUT; +} + +sunindextype OCTSparseMatrix_Columns(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_COLUMNS_O(A); + else + return SUNMAT_ILL_INPUT; +} + +sunindextype OCTSparseMatrix_NNZ(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_NNZ_O(A); + else + return SUNMAT_ILL_INPUT; +} + +sunindextype OCTSparseMatrix_NP(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_NP_O(A); + else + return SUNMAT_ILL_INPUT; +} + +int OCTSparseMatrix_SparseType(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_SPARSETYPE_O(A); + else + return SUNMAT_ILL_INPUT; +} + +realtype *OCTSparseMatrix_Data(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_DATA_O(A); + else + return NULL; +} + +sunindextype *OCTSparseMatrix_IndexValues(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_INDEXVALS_O(A); + else + return NULL; +} + +sunindextype *OCTSparseMatrix_IndexPointers(SUNMatrix A) +{ + if (SUNMatGetID(A) == SUNMATRIX_CUSTOM) + return SM_INDEXPTRS_O(A); + else + return NULL; +} + +/* + * ----------------------------------------------------------------- + * implementation of matrix operations + * ----------------------------------------------------------------- + */ + +SUNMatrix_ID OCTMatGetID_Sparse([[maybe_unused]] SUNMatrix A) +{ + return SUNMATRIX_CUSTOM; +} + +void OCTMatDestroy_Sparse(SUNMatrix A) +{ + if (A == NULL) + return; + + // /* free content */ + if (A->content != NULL) + { + /* deleting object created using new */ + delete SM_CONTENT_O(A); + } + A->content = NULL; + + /* free ops and matrix */ + if (A->ops) + { + free(A->ops); + A->ops = NULL; + } + free(A); + A = NULL; + + return; +} + +int OCTMatZero_Sparse(SUNMatrix A) +{ + SparseMatrix *ptr = new SparseMatrix(SM_ROWS_O(A),SM_COLS_O(A)); + SparseMatrix *am = static_cast SM_CONTENT_O(A); + *am = *ptr; + + return SUNMAT_SUCCESS; +} + +/* Function to print the sparse matrix */ + +void OCTSparseMatrix_Print(SUNMatrix A) +{ + SparseMatrix *am = static_cast SM_CONTENT_O(A); + std::cout<<(*am); + return; +} + + +/* ---------------------------------------------------------------------------- + * Function to create a new General Octave linear solver + */ + +SUNLinearSolver OCTLinSol_Gen(N_Vector y, SUNMatrix A ARG_SUNCONTEXT) +{ + SUNLinearSolver S; + OCTLinearSolverContent_GEN content; + + if (OCTSparseMatrix_Rows(A) != OCTSparseMatrix_Columns(A)) return(NULL); + + if ( N_VGetVectorID(y) != SUNDIALS_NVEC_CUSTOM ) + return(NULL); + + if (OCTSparseMatrix_Rows(A) != N_VGetLength(y)) return(NULL); + + + /* Create an empty linear solver */ + S = NULL; + S = SUNLinSolNewEmpty(SP_OCTAVE_SUNCONTEXT); + if (S == NULL) return(NULL); + + /* Attach operations */ + S->ops->gettype = OCTLinSolGetType_Gen; + S->ops->getid = OCTLinSolGetID_Gen; + S->ops->initialize = NULL; + S->ops->setup = NULL; + S->ops->solve = OCTLinSolSolve_Gen; + S->ops->lastflag = NULL; + S->ops->space = NULL; + S->ops->free = OCTLinSolFree_Gen; + + /* Create content */ + content = NULL; + content = (OCTLinearSolverContent_GEN) malloc(sizeof *content); + if (content == NULL) { SUNLinSolFree(S); return(NULL); } + + /* Attach content */ + S->content = content; + + /* Fill content */ + content->last_flag = 0; + content->first_factorize = 1; + + return(S); +} + +/* + * ----------------------------------------------------------------- + * implementation of linear solver operations + * ----------------------------------------------------------------- + */ + +SUNLinearSolver_Type OCTLinSolGetType_Gen([[maybe_unused]] SUNLinearSolver S) +{ + return(SUNLINEARSOLVER_DIRECT); +} + + +SUNLinearSolver_ID OCTLinSolGetID_Gen([[maybe_unused]] SUNLinearSolver S) +{ + return(SUNLINEARSOLVER_CUSTOM); +} + + +int OCTLinSolSolve_Gen(SUNLinearSolver S, SUNMatrix A, N_Vector x, + N_Vector b, [[maybe_unused]] realtype tol) +{ + /* check for valid inputs */ + if ( (A == NULL) || (S == NULL) || (x == NULL) || (b == NULL) ) + return(SUNLS_MEM_NULL); + ColumnVector *xv = static_cast NV_CONTENT_C(x); + ColumnVector *zv = static_cast NV_CONTENT_C(b); + + SparseMatrix *am = static_cast SM_CONTENT_O(A); + + (*xv) = am->solve((*zv)); + + return 1; +} + +int OCTLinSolFree_Gen(SUNLinearSolver S) +{ + /* return with success if already freed */ + if (S == NULL) return(SUNLS_SUCCESS); + + /* delete generic structures */ + if (S->ops) { + free(S->ops); + S->ops = NULL; + } + free(S); S = NULL; + return(SUNLS_SUCCESS); +} + +#ifdef __cplusplus +} /* extern "C" */ +#endif \ No newline at end of file diff --git a/libinterp/dldfcn/oct-sundials.h b/libinterp/dldfcn/oct-sundials.h new file mode 100644 --- /dev/null +++ b/libinterp/dldfcn/oct-sundials.h @@ -0,0 +1,244 @@ +// Copyright (C) 2016-2022 The Octave Project Developers +// +// See the file COPYRIGHT.md in the top-level directory of this +// distribution or . +// +// This file is part of Octave. +// +// Octave is free software: you can redistribute it and/or modify it +// under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// Octave is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with Octave; see the file COPYING. If not, see +// . +// +//////////////////////////////////////////////////////////////////////// + +/* ----------------------------------------------------------------- + * generating custom SUNDIALS implementations + * to be able to allow Octave's native classes and + * solvers to be utilized by the SUNDIALS IDA Class + * ----------------------------------------------------------------- */ + +#ifndef __CONFIG_H_INCLUDE_GUARD +#define __CONFIG_H_INCLUDE_GUARD +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif +#endif + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "oct-locbuf.h" +#include "quit.h" + +/* FIXME : This flag currently checks for a + lot of Sundials features, such as nvecserial, + etc, which might not be needed for just these files. + Maybe a simpler configure flag should be made for + this. */ + +#if defined (HAVE_SUNDIALS) + +# include +# include +# include +# include +# include + +#endif + +/* ----------------------------------------------------------------- + * Useful mathematical macros for computation + * ----------------------------------------------------------------- */ +#define ZERO RCONST(0.0) +#define HALF RCONST(0.5) +#define ONE RCONST(1.0) +#define ONEPT5 RCONST(1.5) +#define TWO RCONST(2.0) +#define TWOTHIRDS RCONST(0.666666666666666666666666666666667) + +/* ----------------------------------------------------------------- + * SUNContext : Since SUNDIALS 6.0.0, all SUNDIALS objects + * (vectors, linear and nonlinear solvers, matrices, etc.) + * hold a reference to a common simulation context object + * defined by the SUNContext class. The macro ensures these + * declarations maintain backward compatibility with versions + * before 6.0.0 as well. + * ----------------------------------------------------------------- */ +# if (HAVE_SUNDIALS_SUNCONTEXT) +# define SP_ARG_SUNCONTEXT SUNContext m_sunContext +# define ARG_SUNCONTEXT , SUNContext m_sunContext +# define OCTAVE_SUNCONTEXT , m_sunContext +# define SP_OCTAVE_SUNCONTEXT m_sunContext +# else +# define ARG_SUNCONTEXT +# define OCTAVE_SUNCONTEXT +# define SP_ARG_SUNCONTEXT +# define SP_OCTAVE_SUNCONTEXT +# endif + +#ifdef __cplusplus /* wrapper to enable C++ usage */ +extern "C" { +#endif + +/* + * ----------------------------------------------------------------- + * Macros for access to octave implementation of N_Vector + * ----------------------------------------------------------------- + */ + +#define NV_CONTENT_C(v) ( (ColumnVector *)(v->content) ) +#define NV_LENGTH_C(v) ( NV_CONTENT_C(v)->numel() ) +#define NV_DATA_C(v) ( NV_CONTENT_C(v)->fortran_vec() ) +#define NV_Ith_C(v,i) ( NV_DATA_C(v)[i]) + +/* ----------------------------------------------------------------- + * NVECTOR API functions exported for Octave's custom implementation + * ----------------------------------------------------------------- */ + +SUNDIALS_EXPORT N_Vector N_VNew_Octave(int vec_length ARG_SUNCONTEXT ); +SUNDIALS_EXPORT N_Vector N_VNewEmpty_Octave ( SP_ARG_SUNCONTEXT ); +SUNDIALS_EXPORT N_Vector N_VMake_Octave(const ColumnVector& v ARG_SUNCONTEXT ); + +SUNDIALS_EXPORT sunindextype N_VGetLength_Octave(N_Vector v); +SUNDIALS_EXPORT void N_VPrint_Octave(N_Vector v); +SUNDIALS_EXPORT N_Vector_ID N_VGetVectorID_Octave(N_Vector v); +SUNDIALS_EXPORT N_Vector N_VCloneEmpty_Octave(N_Vector w); +SUNDIALS_EXPORT N_Vector N_VClone_Octave(N_Vector w); +SUNDIALS_EXPORT void N_VDestroy_Octave(N_Vector v); +SUNDIALS_EXPORT void N_VSpace_Octave(N_Vector v, sunindextype *lrw, sunindextype *liw); +SUNDIALS_EXPORT double *N_VGetArrayPointer_Octave(N_Vector v); +SUNDIALS_EXPORT void N_VSetArrayPointer_Octave(realtype *v_data, N_Vector v); + +/* standard vector operations */ +SUNDIALS_EXPORT void N_VLinearSum_Octave(double a, N_Vector x, double b, N_Vector y, N_Vector z); +SUNDIALS_EXPORT void N_VConst_Octave(double c, N_Vector z); +SUNDIALS_EXPORT void N_VProd_Octave(N_Vector x, N_Vector y, N_Vector z); +SUNDIALS_EXPORT void N_VDiv_Octave(N_Vector x, N_Vector y, N_Vector z); +SUNDIALS_EXPORT void N_VScale_Octave(double c, N_Vector x, N_Vector z); +SUNDIALS_EXPORT void N_VAbs_Octave(N_Vector x, N_Vector z); +SUNDIALS_EXPORT void N_VInv_Octave(N_Vector x, N_Vector z); +SUNDIALS_EXPORT void N_VAddConst_Octave(N_Vector x, double b, N_Vector z); +SUNDIALS_EXPORT double N_VDotProd_Octave(N_Vector x, N_Vector y); +SUNDIALS_EXPORT double N_VMaxNorm_Octave(N_Vector x); +SUNDIALS_EXPORT double N_VWrmsNorm_Octave(N_Vector x, N_Vector w); +SUNDIALS_EXPORT double N_VWrmsNormMask_Octave(N_Vector x, N_Vector w, N_Vector id); +SUNDIALS_EXPORT double N_VMin_Octave(N_Vector x); +SUNDIALS_EXPORT void N_VCompare_Octave(double c, N_Vector x, N_Vector z); +SUNDIALS_EXPORT booleantype N_VInvTest_Octave(N_Vector x, N_Vector z); +SUNDIALS_EXPORT booleantype N_VConstrMask_Octave(N_Vector c, N_Vector x, N_Vector m); +SUNDIALS_EXPORT double N_VMinQuotient_Octave(N_Vector num, N_Vector denom); + +// /* fused vector operations */ +SUNDIALS_EXPORT int N_VScaleAddMulti_Octave(int nvec, double* a, N_Vector x, + N_Vector* Y, N_Vector* Z); + +/* vector array operations */ +SUNDIALS_EXPORT int N_VLinearSumVectorArray_Octave(int nvec, + double a, N_Vector* X, + double b, N_Vector* Y, + N_Vector* Z); +SUNDIALS_EXPORT int N_VScaleVectorArray_Octave(int nvec, double* c, + N_Vector* X, N_Vector* Z); + +/* OPTIONAL local reduction kernels (no parallel communication) */ +SUNDIALS_EXPORT double N_VWSqrSumLocal_Octave(N_Vector x, N_Vector w); +SUNDIALS_EXPORT double N_VWSqrSumMaskLocal_Octave(N_Vector x, N_Vector w, N_Vector id); + +/* + * ----------------------------------------------------------------- + * Enable / disable fused vector operations + * ----------------------------------------------------------------- + */ +SUNDIALS_EXPORT int N_VEnableFusedOps_Octave(N_Vector v, booleantype tf); + + +/* ------------------------------------------------------ + * Macros for access to Octave implementation of SUNMATRIX + * ----------------------------------------------------- */ +#define CSC_MAT 0 /* Matrix Type Definition */ +#define SM_CONTENT_O(A) ( (SparseMatrix *)(A->content) ) +#define SM_ROWS_O(A) ( SM_CONTENT_O(A)->rows() ) +#define SM_COLS_O(A) ( SM_CONTENT_O(A)->cols() ) +#define SM_COLUMNS_O(A) ( SM_CONTENT_O(A)->columns() ) +#define SM_NNZ_O(A) ( SM_CONTENT_O(A)->nnz() ) +/* directly setting =N as Octave only supports CSC Representation */ +#define SM_NP_O(A) ( SM_CONTENT_O(A)->columns() ) +#define SM_SPARSETYPE_O(A) ( CSC_MAT ) +#define SM_DATA_O(A) ( SM_CONTENT_O(A)->data() ) +#define SM_INDEXVALS_O(A) ( SM_CONTENT_O(A)->ridx() ) +#define SM_INDEXPTRS_O(A) ( SM_CONTENT_O(A)->cidx() ) + +/* ---------------------------------------- + * SUNMATRIX API functions exported for Octave's + * custom implementation + * ---------------------------------------- */ + +SUNDIALS_EXPORT SUNMatrix OCTSparseMatrix(sunindextype M, sunindextype N, + sunindextype NNZ ARG_SUNCONTEXT ); +SUNDIALS_EXPORT void OCTSparseMatrix_Print(SUNMatrix A); + +SUNDIALS_EXPORT sunindextype OCTSparseMatrix_Rows(SUNMatrix A); +SUNDIALS_EXPORT sunindextype OCTSparseMatrix_Columns(SUNMatrix A); +SUNDIALS_EXPORT sunindextype OCTSparseMatrix_NNZ(SUNMatrix A); +SUNDIALS_EXPORT sunindextype OCTSparseMatrix_NP(SUNMatrix A); +SUNDIALS_EXPORT int OCTSparseMatrix_SparseType(SUNMatrix A); +SUNDIALS_EXPORT realtype* OCTSparseMatrix_Data(SUNMatrix A); +SUNDIALS_EXPORT sunindextype* OCTSparseMatrix_IndexValues(SUNMatrix A); +SUNDIALS_EXPORT sunindextype* OCTSparseMatrix_IndexPointers(SUNMatrix A); + +SUNDIALS_EXPORT SUNMatrix_ID OCTMatGetID_Sparse(SUNMatrix A); +SUNDIALS_EXPORT void OCTMatDestroy_Sparse (SUNMatrix A); +SUNDIALS_EXPORT int OCTMatZero_Sparse (SUNMatrix A); + +/* + * ---------------------------------------------------------------------------- + * Octave implementation of SUNLinearSolver + * ---------------------------------------------------------------------------- + */ +struct _OCTLinearSolverContent_GEN { + int last_flag; + int first_factorize; +}; + +/* -------------------------------------------------- + * Declarations of Octave's custom implementation of + * a Matrix Linear Solver + * ------------------------------------------------- */ +typedef struct _OCTLinearSolverContent_GEN *OCTLinearSolverContent_GEN; + +/* ------------------------------------- + * SUNLINEARSOLVER API functions exported for Octave's + * custom implementation + * -----------------------------------*/ + +SUNDIALS_EXPORT SUNLinearSolver OCTLinSol_Gen(N_Vector y, SUNMatrix A ARG_SUNCONTEXT ); +SUNDIALS_EXPORT SUNLinearSolver_Type OCTLinSolGetType_Gen(SUNLinearSolver S); +SUNDIALS_EXPORT SUNLinearSolver_ID OCTLinSolGetID_Gen(SUNLinearSolver S); +SUNDIALS_EXPORT int OCTLinSolSolve_Gen(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol); +SUNDIALS_EXPORT int OCTLinSolFree_Gen(SUNLinearSolver S); + + +#ifdef __cplusplus +} +#endif \ No newline at end of file diff --git a/scripts/ode/ode15i.m b/scripts/ode/ode15i.m --- a/scripts/ode/ode15i.m +++ b/scripts/ode/ode15i.m @@ -455,7 +455,7 @@ %! assert ([t(end), y(end,:)], fref, 1e-3); ## Jacobian fcn sparse -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! opt = odeset ("Jacobian", @jacfunsparse, "AbsTol", 1e-7, "RelTol", 1e-7); %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt); %! assert ([t(end), y(end,:)], fref, 1e-3); @@ -550,7 +550,7 @@ %! '"Jacobian" matrices must be real square matrices'); ## Jacobian cell sparse wrong dimension -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! DFDY = sparse ([-0.04, 1; %! 0.04, 1]); %! DFDYP = sparse ([-1, 0, 0; diff --git a/scripts/ode/ode15s.m b/scripts/ode/ode15s.m --- a/scripts/ode/ode15s.m +++ b/scripts/ode/ode15s.m @@ -554,21 +554,21 @@ %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! opt = odeset ("MStateDependence", "none", %! "Mass", [1, 0, 0; 0, 1, 0; 0, 0, 0], %! "Jacobian", @jacfunsparse); %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! opt = odeset ("MStateDependence", "none", %! "Mass", sparse ([1, 0, 0; 0, 1, 0; 0, 0, 0]), %! "Jacobian", @jacfunsparse); %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! warning ("off", "ode15s:mass_state_dependent_provided", "local"); %! opt = odeset ("MStateDependence", "none", %! "Mass", @massdensefunstate, @@ -584,14 +584,14 @@ %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! opt = odeset ("MStateDependence", "none", %! "Mass", @massdensefuntime, %! "Jacobian", @jacfunsparse); %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_SUNLINSOL_KLU +%!testif HAVE_SUNDIALS %! opt = odeset ("MStateDependence", "none", %! "Mass", @masssparsefuntime, %! "Jacobian", @jacfunsparse);