Mercurial > octave

comparison libinterp/dldfcn/__ode15__.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
85 OCTAVE_BEGIN_NAMESPACE(octave) 85 OCTAVE_BEGIN_NAMESPACE(octave)
86 86
87 #if defined (HAVE_SUNDIALS) 87 #if defined (HAVE_SUNDIALS)
88 88
89 # if ! defined (HAVE_IDASETJACFN) && defined (HAVE_IDADLSSETJACFN) 89 # if ! defined (HAVE_IDASETJACFN) && defined (HAVE_IDADLSSETJACFN)
90 static inline int 90 static inline int
91 IDASetJacFn (void *ida_mem, IDADlsJacFn jac) 91 IDASetJacFn (void *ida_mem, IDADlsJacFn jac)
92 { 92 {
93 return IDADlsSetJacFn (ida_mem, jac); 93 return IDADlsSetJacFn (ida_mem, jac);
94 } 94 }
95 # endif 95 # endif
96 96
97 # if ! defined (HAVE_IDASETLINEARSOLVER) && defined (HAVE_IDADLSSETLINEARSOLVER) 97 # if ! defined (HAVE_IDASETLINEARSOLVER) && defined (HAVE_IDADLSSETLINEARSOLVER)
98 static inline int 98 static inline int
99 IDASetLinearSolver (void *ida_mem, SUNLinearSolver LS, SUNMatrix A) 99 IDASetLinearSolver (void *ida_mem, SUNLinearSolver LS, SUNMatrix A)
100 { 100 {
101 return IDADlsSetLinearSolver (ida_mem, LS, A); 101 return IDADlsSetLinearSolver (ida_mem, LS, A);
102 } 102 }
103 # endif 103 # endif
104 104
105 # if ! defined (HAVE_SUNLINSOL_DENSE) && defined (HAVE_SUNDENSELINEARSOLVER) 105 # if ! defined (HAVE_SUNLINSOL_DENSE) && defined (HAVE_SUNDENSELINEARSOLVER)
106 static inline SUNLinearSolver 106 static inline SUNLinearSolver
107 SUNLinSol_Dense (N_Vector y, SUNMatrix A) 107 SUNLinSol_Dense (N_Vector y, SUNMatrix A)
108 { 108 {
109 return SUNDenseLinearSolver (y, A); 109 return SUNDenseLinearSolver (y, A);
110 } 110 }
111 # endif 111 # endif
112 112
113 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) 113 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
114 # if ! defined (HAVE_SUNLINSOL_KLU) && defined (HAVE_SUNKLU) 114 # if ! defined (HAVE_SUNLINSOL_KLU) && defined (HAVE_SUNKLU)
115 static inline SUNLinearSolver 115 static inline SUNLinearSolver
116 SUNLinSol_KLU (N_Vector y, SUNMatrix A) 116 SUNLinSol_KLU (N_Vector y, SUNMatrix A)
117 { 117 {
118 return SUNKLU (y, A); 118 return SUNKLU (y, A);
119 } 119 }
120 # endif 120 # endif
121 # endif 121 # endif
122 122
123 static inline realtype * 123 static inline realtype *
124 nv_data_s (N_Vector& v) 124 nv_data_s (N_Vector& v)
125 { 125 {
126 # if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) 126 # if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
127 // Disable warning from GCC about old-style casts in Sundials 127 // Disable warning from GCC about old-style casts in Sundials
128 // macro expansions. Do this in a function so that this 128 // macro expansions. Do this in a function so that this
129 // diagnostic may still be enabled for the rest of the file. 129 // diagnostic may still be enabled for the rest of the file.
130 # pragma GCC diagnostic push 130 # pragma GCC diagnostic push
131 # pragma GCC diagnostic ignored "-Wold-style-cast" 131 # pragma GCC diagnostic ignored "-Wold-style-cast"
132 # endif 132 # endif
133 133
134 return NV_DATA_S (v); 134 return NV_DATA_S (v);
135 135
136 # if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC) 136 # if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
137 // Restore prevailing warning state for remainder of the file. 137 // Restore prevailing warning state for remainder of the file.
138 # pragma GCC diagnostic pop 138 # pragma GCC diagnostic pop
139 # endif 139 # endif
140 }
141
142 class IDA
143 {
144 public:
145
146 typedef
147 ColumnVector (*DAERHSFuncIDA) (const ColumnVector& x,
148 const ColumnVector& xdot,
149 realtype t, const octave_value& idaf);
150
151 typedef
152 Matrix (*DAEJacFuncDense) (const ColumnVector& x,
153 const ColumnVector& xdot, realtype t,
154 realtype cj, const octave_value& idaj);
155
156 typedef
157 SparseMatrix (*DAEJacFuncSparse) (const ColumnVector& x,
158 const ColumnVector& xdot,
159 realtype t, realtype cj,
160 const octave_value& idaj);
161
162 typedef
163 Matrix (*DAEJacCellDense) (Matrix *dfdy, Matrix *dfdyp,
164 realtype cj);
165
166 typedef
167 SparseMatrix (*DAEJacCellSparse) (SparseMatrix *dfdy,
168 SparseMatrix *dfdyp, realtype cj);
169
170 //Default
171 IDA (void)
172 : m_t0 (0.0), m_y0 (), m_yp0 (), m_havejac (false), m_havejacfcn (false),
173 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fcn (),
174 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
175 m_spdfdyp (nullptr), m_fcn (nullptr), m_jacfcn (nullptr),
176 m_jacspfcn (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr),
177 m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
178 { }
179
180
181 IDA (realtype t, ColumnVector y, ColumnVector yp,
182 const octave_value& ida_fcn, DAERHSFuncIDA daefun)
183 : m_t0 (t), m_y0 (y), m_yp0 (yp), m_havejac (false), m_havejacfcn (false),
184 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fcn (ida_fcn),
185 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
186 m_spdfdyp (nullptr), m_fcn (daefun), m_jacfcn (nullptr),
187 m_jacspfcn (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr),
188 m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
189 { }
190
191
192 ~IDA (void)
193 {
194 IDAFree (&m_mem);
195 SUNLinSolFree (m_sunLinearSolver);
196 SUNMatDestroy (m_sunJacMatrix);
197 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
198 SUNContext_Free (&m_sunContext);
199 # endif
140 } 200 }
141 201
142 class IDA 202 IDA&
203 set_jacobian (const octave_value& jac, DAEJacFuncDense j)
143 { 204 {
144 public: 205 m_jacfcn = j;
145 206 m_ida_jac = jac;
146 typedef 207 m_havejac = true;
147 ColumnVector (*DAERHSFuncIDA) (const ColumnVector& x, 208 m_havejacfcn = true;
148 const ColumnVector& xdot, 209 m_havejacsparse = false;
149 realtype t, const octave_value& idaf); 210
150 211 return *this;
151 typedef 212 }
152 Matrix (*DAEJacFuncDense) (const ColumnVector& x, 213
153 const ColumnVector& xdot, realtype t, 214 IDA&
154 realtype cj, const octave_value& idaj); 215 set_jacobian (const octave_value& jac, DAEJacFuncSparse j)
155 216 {
156 typedef 217 m_jacspfcn = j;
157 SparseMatrix (*DAEJacFuncSparse) (const ColumnVector& x, 218 m_ida_jac = jac;
158 const ColumnVector& xdot, 219 m_havejac = true;
159 realtype t, realtype cj, 220 m_havejacfcn = true;
160 const octave_value& idaj); 221 m_havejacsparse = true;
161 222
162 typedef 223 return *this;
163 Matrix (*DAEJacCellDense) (Matrix *dfdy, Matrix *dfdyp, 224 }
164 realtype cj); 225
165 226 IDA&
166 typedef 227 set_jacobian (Matrix *dy, Matrix *dyp, DAEJacCellDense j)
167 SparseMatrix (*DAEJacCellSparse) (SparseMatrix *dfdy, 228 {
168 SparseMatrix *dfdyp, realtype cj); 229 m_jacdcell = j;
169 230 m_dfdy = dy;
170 //Default 231 m_dfdyp = dyp;
171 IDA (void) 232 m_havejac = true;
172 : m_t0 (0.0), m_y0 (), m_yp0 (), m_havejac (false), m_havejacfcn (false), 233 m_havejacfcn = false;
173 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fcn (), 234 m_havejacsparse = false;
174 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr), 235
175 m_spdfdyp (nullptr), m_fcn (nullptr), m_jacfcn (nullptr), 236 return *this;
176 m_jacspfcn (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr), 237 }
177 m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr) 238
178 { } 239 IDA&
179 240 set_jacobian (SparseMatrix *dy, SparseMatrix *dyp,
180 241 DAEJacCellSparse j)
181 IDA (realtype t, ColumnVector y, ColumnVector yp, 242 {
182 const octave_value& ida_fcn, DAERHSFuncIDA daefun) 243 m_jacspcell = j;
183 : m_t0 (t), m_y0 (y), m_yp0 (yp), m_havejac (false), m_havejacfcn (false), 244 m_spdfdy = dy;
184 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fcn (ida_fcn), 245 m_spdfdyp = dyp;
185 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr), 246 m_havejac = true;
186 m_spdfdyp (nullptr), m_fcn (daefun), m_jacfcn (nullptr), 247 m_havejacfcn = false;
187 m_jacspfcn (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr), 248 m_havejacsparse = true;
188 m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr) 249
189 { } 250 return *this;
190 251 }
191 252
192 ~IDA (void) 253 void set_userdata (void);
193 { 254
194 IDAFree (&m_mem); 255 void initialize (void);
195 SUNLinSolFree (m_sunLinearSolver); 256
196 SUNMatDestroy (m_sunJacMatrix); 257 static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n);
258
197 # if defined (HAVE_SUNDIALS_SUNCONTEXT) 259 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
198 SUNContext_Free (&m_sunContext); 260 N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
199 # endif
200 }
201
202 IDA&
203 set_jacobian (const octave_value& jac, DAEJacFuncDense j)
204 {
205 m_jacfcn = j;
206 m_ida_jac = jac;
207 m_havejac = true;
208 m_havejacfcn = true;
209 m_havejacsparse = false;
210
211 return *this;
212 }
213
214 IDA&
215 set_jacobian (const octave_value& jac, DAEJacFuncSparse j)
216 {
217 m_jacspfcn = j;
218 m_ida_jac = jac;
219 m_havejac = true;
220 m_havejacfcn = true;
221 m_havejacsparse = true;
222
223 return *this;
224 }
225
226 IDA&
227 set_jacobian (Matrix *dy, Matrix *dyp, DAEJacCellDense j)
228 {
229 m_jacdcell = j;
230 m_dfdy = dy;
231 m_dfdyp = dyp;
232 m_havejac = true;
233 m_havejacfcn = false;
234 m_havejacsparse = false;
235
236 return *this;
237 }
238
239 IDA&
240 set_jacobian (SparseMatrix *dy, SparseMatrix *dyp,
241 DAEJacCellSparse j)
242 {
243 m_jacspcell = j;
244 m_spdfdy = dy;
245 m_spdfdyp = dyp;
246 m_havejac = true;
247 m_havejacfcn = false;
248 m_havejacsparse = true;
249
250 return *this;
251 }
252
253 void set_userdata (void);
254
255 void initialize (void);
256
257 static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n);
258
259 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
260 N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
261 # else 261 # else
262 static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n); 262 static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
263 # endif 263 # endif
264 264
265 void 265 void
266 set_up (const ColumnVector& y); 266 set_up (const ColumnVector& y);
267 267
268 void 268 void
269 set_tolerance (ColumnVector& abstol, realtype reltol); 269 set_tolerance (ColumnVector& abstol, realtype reltol);
270 270
271 void 271 void
272 set_tolerance (realtype abstol, realtype reltol); 272 set_tolerance (realtype abstol, realtype reltol);
273 273
274 static int 274 static int
275 resfun (realtype t, N_Vector yy, N_Vector yyp, 275 resfun (realtype t, N_Vector yy, N_Vector yyp,
276 N_Vector rr, void *user_data); 276 N_Vector rr, void *user_data);
277 277
278 void 278 void
279 resfun_impl (realtype t, N_Vector& yy, 279 resfun_impl (realtype t, N_Vector& yy,
280 N_Vector& yyp, N_Vector& rr); 280 N_Vector& yyp, N_Vector& rr);
281 static int 281 static int
282 jacdense (realtype t, realtype cj, N_Vector yy, N_Vector yyp, 282 jacdense (realtype t, realtype cj, N_Vector yy, N_Vector yyp,
283 N_Vector, SUNMatrix JJ, void *user_data, N_Vector, 283 N_Vector, SUNMatrix JJ, void *user_data, N_Vector,
284 N_Vector, N_Vector) 284 N_Vector, N_Vector)
285 {
286 IDA *self = static_cast <IDA *> (user_data);
287 self->jacdense_impl (t, cj, yy, yyp, JJ);
288 return 0;
289 }
290
291 void
292 jacdense_impl (realtype t, realtype cj,
293 N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ);
294
295 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
296 static int
297 jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp,
298 N_Vector, SUNMatrix Jac, void *user_data, N_Vector,
299 N_Vector, N_Vector)
300 {
301 IDA *self = static_cast <IDA *> (user_data);
302 self->jacsparse_impl (t, cj, yy, yyp, Jac);
303 return 0;
304 }
305
306 void
307 jacsparse_impl (realtype t, realtype cj,
308 N_Vector& yy, N_Vector& yyp, SUNMatrix& Jac);
309 # endif
310
311 void set_maxstep (realtype maxstep);
312
313 void set_initialstep (realtype initialstep);
314
315 bool
316 interpolate (octave_idx_type& cont, Matrix& output, ColumnVector& tout,
317 int refine, realtype tend, bool haveoutputfcn,
318 bool haveoutputsel, const octave_value& output_fcn,
319 ColumnVector& outputsel, bool haveeventfunction,
320 const octave_value& event_fcn, ColumnVector& te,
321 Matrix& ye, ColumnVector& ie, ColumnVector& oldval,
322 ColumnVector& oldisterminal, ColumnVector& olddir,
323 octave_idx_type& temp, ColumnVector& yold,
324 const octave_idx_type num_event_args);
325
326 bool
327 outputfun (const octave_value& output_fcn, bool haveoutputsel,
328 const ColumnVector& output, realtype tout, realtype tend,
329 ColumnVector& outputsel, const std::string& flag);
330
331
332 bool
333 event (const octave_value& event_fcn,
334 ColumnVector& te, Matrix& ye, ColumnVector& ie,
335 realtype tsol, const ColumnVector& y, const std::string& flag,
336 const ColumnVector& yp, ColumnVector& oldval,
337 ColumnVector& oldisterminal, ColumnVector& olddir,
338 octave_idx_type cont, octave_idx_type& temp, realtype told,
339 ColumnVector& yold,
340 const octave_idx_type num_event_args);
341
342 void set_maxorder (int maxorder);
343
344 octave_value_list
345 integrate (const octave_idx_type numt, const ColumnVector& tt,
346 const ColumnVector& y0, const ColumnVector& yp0,
347 const int refine, bool haverefine, bool haveoutputfcn,
348 const octave_value& output_fcn, bool haveoutputsel,
349 ColumnVector& outputsel, bool haveeventfunction,
350 const octave_value& event_fcn,
351 const octave_idx_type num_event_args);
352
353 void print_stat (void);
354
355 private:
356
357 realtype m_t0;
358 ColumnVector m_y0;
359 ColumnVector m_yp0;
360 bool m_havejac;
361 bool m_havejacfcn;
362 bool m_havejacsparse;
363 void *m_mem;
364 octave_f77_int_type m_num;
365 octave_value m_ida_fcn;
366 octave_value m_ida_jac;
367 Matrix *m_dfdy;
368 Matrix *m_dfdyp;
369 SparseMatrix *m_spdfdy;
370 SparseMatrix *m_spdfdyp;
371 DAERHSFuncIDA m_fcn;
372 DAEJacFuncDense m_jacfcn;
373 DAEJacFuncSparse m_jacspfcn;
374 DAEJacCellDense m_jacdcell;
375 DAEJacCellSparse m_jacspcell;
376 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
377 SUNContext m_sunContext;
378 # endif
379 SUNMatrix m_sunJacMatrix;
380 SUNLinearSolver m_sunLinearSolver;
381 };
382
383 int
384 IDA::resfun (realtype t, N_Vector yy, N_Vector yyp, N_Vector rr,
385 void *user_data)
386 { 285 {
387 IDA *self = static_cast <IDA *> (user_data); 286 IDA *self = static_cast <IDA *> (user_data);
388 self->resfun_impl (t, yy, yyp, rr); 287 self->jacdense_impl (t, cj, yy, yyp, JJ);
389 return 0; 288 return 0;
390 } 289 }
391 290
392 void 291 void
393 IDA::resfun_impl (realtype t, N_Vector& yy, 292 jacdense_impl (realtype t, realtype cj,
394 N_Vector& yyp, N_Vector& rr) 293 N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ);
294
295 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
296 static int
297 jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp,
298 N_Vector, SUNMatrix Jac, void *user_data, N_Vector,
299 N_Vector, N_Vector)
395 { 300 {
396 ColumnVector y = IDA::NVecToCol (yy, m_num); 301 IDA *self = static_cast <IDA *> (user_data);
397 302 self->jacsparse_impl (t, cj, yy, yyp, Jac);
398 ColumnVector yp = IDA::NVecToCol (yyp, m_num); 303 return 0;
399
400 ColumnVector res = (*m_fcn) (y, yp, t, m_ida_fcn);
401
402 realtype *puntrr = nv_data_s (rr);
403
404 for (octave_idx_type i = 0; i < m_num; i++)
405 puntrr[i] = res(i);
406 } 304 }
305
306 void
307 jacsparse_impl (realtype t, realtype cj,
308 N_Vector& yy, N_Vector& yyp, SUNMatrix& Jac);
309 # endif
310
311 void set_maxstep (realtype maxstep);
312
313 void set_initialstep (realtype initialstep);
314
315 bool
316 interpolate (octave_idx_type& cont, Matrix& output, ColumnVector& tout,
317 int refine, realtype tend, bool haveoutputfcn,
318 bool haveoutputsel, const octave_value& output_fcn,
319 ColumnVector& outputsel, bool haveeventfunction,
320 const octave_value& event_fcn, ColumnVector& te,
321 Matrix& ye, ColumnVector& ie, ColumnVector& oldval,
322 ColumnVector& oldisterminal, ColumnVector& olddir,
323 octave_idx_type& temp, ColumnVector& yold,
324 const octave_idx_type num_event_args);
325
326 bool
327 outputfun (const octave_value& output_fcn, bool haveoutputsel,
328 const ColumnVector& output, realtype tout, realtype tend,
329 ColumnVector& outputsel, const std::string& flag);
330
331
332 bool
333 event (const octave_value& event_fcn,
334 ColumnVector& te, Matrix& ye, ColumnVector& ie,
335 realtype tsol, const ColumnVector& y, const std::string& flag,
336 const ColumnVector& yp, ColumnVector& oldval,
337 ColumnVector& oldisterminal, ColumnVector& olddir,
338 octave_idx_type cont, octave_idx_type& temp, realtype told,
339 ColumnVector& yold,
340 const octave_idx_type num_event_args);
341
342 void set_maxorder (int maxorder);
343
344 octave_value_list
345 integrate (const octave_idx_type numt, const ColumnVector& tt,
346 const ColumnVector& y0, const ColumnVector& yp0,
347 const int refine, bool haverefine, bool haveoutputfcn,
348 const octave_value& output_fcn, bool haveoutputsel,
349 ColumnVector& outputsel, bool haveeventfunction,
350 const octave_value& event_fcn,
351 const octave_idx_type num_event_args);
352
353 void print_stat (void);
354
355 private:
356
357 realtype m_t0;
358 ColumnVector m_y0;
359 ColumnVector m_yp0;
360 bool m_havejac;
361 bool m_havejacfcn;
362 bool m_havejacsparse;
363 void *m_mem;
364 octave_f77_int_type m_num;
365 octave_value m_ida_fcn;
366 octave_value m_ida_jac;
367 Matrix *m_dfdy;
368 Matrix *m_dfdyp;
369 SparseMatrix *m_spdfdy;
370 SparseMatrix *m_spdfdyp;
371 DAERHSFuncIDA m_fcn;
372 DAEJacFuncDense m_jacfcn;
373 DAEJacFuncSparse m_jacspfcn;
374 DAEJacCellDense m_jacdcell;
375 DAEJacCellSparse m_jacspcell;
376 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
377 SUNContext m_sunContext;
378 # endif
379 SUNMatrix m_sunJacMatrix;
380 SUNLinearSolver m_sunLinearSolver;
381 };
382
383 int
384 IDA::resfun (realtype t, N_Vector yy, N_Vector yyp, N_Vector rr,
385 void *user_data)
386 {
387 IDA *self = static_cast <IDA *> (user_data);
388 self->resfun_impl (t, yy, yyp, rr);
389 return 0;
390 }
391
392 void
393 IDA::resfun_impl (realtype t, N_Vector& yy,
394 N_Vector& yyp, N_Vector& rr)
395 {
396 ColumnVector y = IDA::NVecToCol (yy, m_num);
397
398 ColumnVector yp = IDA::NVecToCol (yyp, m_num);
399
400 ColumnVector res = (*m_fcn) (y, yp, t, m_ida_fcn);
401
402 realtype *puntrr = nv_data_s (rr);
403
404 for (octave_idx_type i = 0; i < m_num; i++)
405 puntrr[i] = res(i);
406 }
407 407
408 # if defined (HAVE_SUNDIALS_SUNCONTEXT) 408 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
409 # define OCTAVE_SUNCONTEXT , m_sunContext 409 # define OCTAVE_SUNCONTEXT , m_sunContext
410 # else 410 # else
411 # define OCTAVE_SUNCONTEXT 411 # define OCTAVE_SUNCONTEXT
412 # endif 412 # endif
413 413
414 void 414 void
415 IDA::set_up (const ColumnVector& y) 415 IDA::set_up (const ColumnVector& y)
416 { 416 {
417 N_Vector yy = ColToNVec (y, m_num); 417 N_Vector yy = ColToNVec (y, m_num);
418 418
419 if (m_havejacsparse) 419 if (m_havejacsparse)
420 { 420 {
421 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) 421 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
422 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE) 422 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE)
423 // Initially allocate memory for 0 entries. We will reallocate when we 423 // Initially allocate memory for 0 entries. We will reallocate when we
424 // get the Jacobian matrix from the user and know the actual number of 424 // get the Jacobian matrix from the user and know the actual number of
425 // entries. 425 // entries.
426 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT 426 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT
427 OCTAVE_SUNCONTEXT); 427 OCTAVE_SUNCONTEXT);
428 # else 428 # else
429 octave_f77_int_type max_elems; 429 octave_f77_int_type max_elems;
430 if (math::int_multiply_overflow (m_num, m_num, &max_elems)) 430 if (math::int_multiply_overflow (m_num, m_num, &max_elems))
431 error ("Unable to allocate memory for sparse Jacobian"); 431 error ("Unable to allocate memory for sparse Jacobian");
432 432
433 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT 433 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT
434 OCTAVE_SUNCONTEXT); 434 OCTAVE_SUNCONTEXT);
435 # endif 435 # endif
436 if (! m_sunJacMatrix) 436 if (! m_sunJacMatrix)
437 error ("Unable to create sparse Jacobian for Sundials"); 437 error ("Unable to create sparse Jacobian for Sundials");
438 438
439 m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix 439 m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix
440 OCTAVE_SUNCONTEXT);
441 if (! m_sunLinearSolver)
442 error ("Unable to create KLU sparse solver");
443
444 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
445 error ("Unable to set sparse linear solver");
446
447 IDASetJacFn (m_mem, IDA::jacsparse);
448
449 # else
450 error ("SUNDIALS SUNLINSOL KLU was unavailable or disabled when "
451 "Octave was built");
452
453 # endif
454
455 }
456 else
457 {
458
459 m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT);
460 if (! m_sunJacMatrix)
461 error ("Unable to create dense Jacobian for Sundials");
462
463 m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix
440 OCTAVE_SUNCONTEXT); 464 OCTAVE_SUNCONTEXT);
441 if (! m_sunLinearSolver) 465 if (! m_sunLinearSolver)
442 error ("Unable to create KLU sparse solver"); 466 error ("Unable to create dense linear solver");
443 467
444 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix)) 468 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
445 error ("Unable to set sparse linear solver"); 469 error ("Unable to set dense linear solver");
446 470
447 IDASetJacFn (m_mem, IDA::jacsparse); 471 if (m_havejac && IDASetJacFn (m_mem, IDA::jacdense) != 0)
448 472 error ("Unable to set dense Jacobian function");
473
474 }
475 }
476
477 void
478 IDA::jacdense_impl (realtype t, realtype cj,
479 N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ)
480
481 {
482 octave_f77_int_type Neq = NV_LENGTH_S (yy);
483
484 ColumnVector y = NVecToCol (yy, Neq);
485
486 ColumnVector yp = NVecToCol (yyp, Neq);
487
488 Matrix jac;
489
490 if (m_havejacfcn)
491 jac = (*m_jacfcn) (y, yp, t, cj, m_ida_jac);
492 else
493 jac = (*m_jacdcell) (m_dfdy, m_dfdyp, cj);
494
495 octave_f77_int_type num_jac = to_f77_int (jac.numel ());
496 std::copy (jac.fortran_vec (),
497 jac.fortran_vec () + num_jac,
498 SUNDenseMatrix_Data (JJ));
499 }
500
501 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
502 void
503 IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp,
504 SUNMatrix& Jac)
505
506 {
507 ColumnVector y = NVecToCol (yy, m_num);
508
509 ColumnVector yp = NVecToCol (yyp, m_num);
510
511 SparseMatrix jac;
512
513 if (m_havejacfcn)
514 jac = (*m_jacspfcn) (y, yp, t, cj, m_ida_jac);
515 else
516 jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj);
517
518 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE)
519 octave_f77_int_type nnz = to_f77_int (jac.nnz ());
520 if (nnz > SUNSparseMatrix_NNZ (Jac))
521 {
522 // Allocate memory for sparse Jacobian defined in user function.
523 // This will always be required at least once since we set the number
524 // of non-zero elements to zero initially.
525 if (SUNSparseMatrix_Reallocate (Jac, nnz))
526 error ("Unable to allocate sufficient memory for IDA sparse matrix");
527 }
528 # endif
529
530 SUNMatZero_Sparse (Jac);
531 // We have to use "sunindextype *" here but still need to check that
532 // conversion of each element to "octave_f77_int_type" is save.
533 sunindextype *colptrs = SUNSparseMatrix_IndexPointers (Jac);
534 sunindextype *rowvals = SUNSparseMatrix_IndexValues (Jac);
535
536 for (octave_f77_int_type i = 0; i < m_num + 1; i++)
537 colptrs[i] = to_f77_int (jac.cidx (i));
538
539 double *d = SUNSparseMatrix_Data (Jac);
540 for (octave_f77_int_type i = 0; i < to_f77_int (jac.nnz ()); i++)
541 {
542 rowvals[i] = to_f77_int (jac.ridx (i));
543 d[i] = jac.data (i);
544 }
545 }
546 # endif
547
548 ColumnVector
549 IDA::NVecToCol (N_Vector& v, octave_f77_int_type n)
550 {
551 ColumnVector data (n);
552 realtype *punt = nv_data_s (v);
553
554 for (octave_f77_int_type i = 0; i < n; i++)
555 data(i) = punt[i];
556
557 return data;
558 }
559
560 N_Vector
561 IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n)
562 {
563 N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT);
564
565 realtype *punt = nv_data_s (v);
566
567 for (octave_f77_int_type i = 0; i < n; i++)
568 punt[i] = data(i);
569
570 return v;
571 }
572
573 void
574 IDA::set_userdata (void)
575 {
576 void *userdata = this;
577
578 if (IDASetUserData (m_mem, userdata) != 0)
579 error ("User data not set");
580 }
581
582 void
583 IDA::initialize (void)
584 {
585 m_num = to_f77_int (m_y0.numel ());
586 # if defined (HAVE_SUNDIALS_SUNCONTEXT)
587 if (SUNContext_Create (nullptr, &m_sunContext) < 0)
588 error ("__ode15__: unable to create context for SUNDIALS");
589 m_mem = IDACreate (m_sunContext);
449 # else 590 # else
450 error ("SUNDIALS SUNLINSOL KLU was unavailable or disabled when " 591 m_mem = IDACreate ();
451 "Octave was built"); 592 # endif
452 593
453 # endif 594 N_Vector yy = ColToNVec (m_y0, m_num);
454 595
455 } 596 N_Vector yyp = ColToNVec (m_yp0, m_num);
456 else 597
457 { 598 IDA::set_userdata ();
458 599
459 m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT); 600 if (IDAInit (m_mem, IDA::resfun, m_t0, yy, yyp) != 0)
460 if (! m_sunJacMatrix) 601 error ("IDA not initialized");
461 error ("Unable to create dense Jacobian for Sundials"); 602 }
462 603
463 m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix 604 void
464 OCTAVE_SUNCONTEXT); 605 IDA::set_tolerance (ColumnVector& abstol, realtype reltol)
465 if (! m_sunLinearSolver) 606 {
466 error ("Unable to create dense linear solver"); 607 N_Vector abs_tol = ColToNVec (abstol, m_num);
467 608
468 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix)) 609 if (IDASVtolerances (m_mem, reltol, abs_tol) != 0)
469 error ("Unable to set dense linear solver"); 610 error ("IDA: Tolerance not set");
470 611
471 if (m_havejac && IDASetJacFn (m_mem, IDA::jacdense) != 0) 612 N_VDestroy_Serial (abs_tol);
472 error ("Unable to set dense Jacobian function"); 613 }
473 614
474 } 615 void
475 } 616 IDA::set_tolerance (realtype abstol, realtype reltol)
476 617 {
477 void 618 if (IDASStolerances (m_mem, reltol, abstol) != 0)
478 IDA::jacdense_impl (realtype t, realtype cj, 619 error ("IDA: Tolerance not set");
479 N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ) 620 }
480 621
481 { 622 octave_value_list
482 octave_f77_int_type Neq = NV_LENGTH_S (yy); 623 IDA::integrate (const octave_idx_type numt, const ColumnVector& tspan,
483 624 const ColumnVector& y, const ColumnVector& yp,
484 ColumnVector y = NVecToCol (yy, Neq); 625 const int refine, bool haverefine, bool haveoutputfcn,
485 626 const octave_value& output_fcn, bool haveoutputsel,
486 ColumnVector yp = NVecToCol (yyp, Neq); 627 ColumnVector& outputsel, bool haveeventfunction,
487 628 const octave_value& event_fcn,
488 Matrix jac; 629 const octave_idx_type num_event_args)
489 630 {
490 if (m_havejacfcn) 631 // Set up output
491 jac = (*m_jacfcn) (y, yp, t, cj, m_ida_jac); 632 ColumnVector tout, yout (m_num), ypout (m_num), ysel (outputsel.numel ());
492 else 633 ColumnVector ie, te, oldval, oldisterminal, olddir;
493 jac = (*m_jacdcell) (m_dfdy, m_dfdyp, cj); 634 Matrix output, ye;
494 635 octave_idx_type cont = 0, temp = 0;
495 octave_f77_int_type num_jac = to_f77_int (jac.numel ()); 636 bool status = false;
496 std::copy (jac.fortran_vec (), 637 std::string string = "";
497 jac.fortran_vec () + num_jac, 638 ColumnVector yold = y;
498 SUNDenseMatrix_Data (JJ)); 639
499 } 640 realtype tsol = tspan(0);
500 641 realtype tend = tspan(numt-1);
501 # if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) 642
502 void 643 N_Vector yyp = ColToNVec (yp, m_num);
503 IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, 644
504 SUNMatrix& Jac) 645 N_Vector yy = ColToNVec (y, m_num);
505 646
506 { 647 // Initialize OutputFcn
507 ColumnVector y = NVecToCol (yy, m_num); 648 if (haveoutputfcn)
508 649 status = IDA::outputfun (output_fcn, haveoutputsel, y,
509 ColumnVector yp = NVecToCol (yyp, m_num); 650 tsol, tend, outputsel, "init");
510 651
511 SparseMatrix jac; 652 // Initialize Events
512 653 if (haveeventfunction)
513 if (m_havejacfcn) 654 status = IDA::event (event_fcn, te, ye, ie, tsol, y,
514 jac = (*m_jacspfcn) (y, yp, t, cj, m_ida_jac); 655 "init", yp, oldval, oldisterminal,
515 else 656 olddir, cont, temp, tsol, yold, num_event_args);
516 jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj); 657
517 658 if (numt > 2)
518 # if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE) 659 {
519 octave_f77_int_type nnz = to_f77_int (jac.nnz ()); 660 // First output value
520 if (nnz > SUNSparseMatrix_NNZ (Jac)) 661 tout.resize (numt);
521 { 662 tout(0) = tsol;
522 // Allocate memory for sparse Jacobian defined in user function. 663 output.resize (numt, m_num);
523 // This will always be required at least once since we set the number 664
524 // of non-zero elements to zero initially. 665 for (octave_idx_type i = 0; i < m_num; i++)
525 if (SUNSparseMatrix_Reallocate (Jac, nnz)) 666 output.elem (0, i) = y.elem (i);
526 error ("Unable to allocate sufficient memory for IDA sparse matrix"); 667
527 } 668 //Main loop
528 # endif 669 for (octave_idx_type j = 1; j < numt && status == 0; j++)
529 670 {
530 SUNMatZero_Sparse (Jac); 671 // IDANORMAL already interpolates tspan(j)
531 // We have to use "sunindextype *" here but still need to check that 672
532 // conversion of each element to "octave_f77_int_type" is save. 673 if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
533 sunindextype *colptrs = SUNSparseMatrix_IndexPointers (Jac); 674 error ("IDASolve failed");
534 sunindextype *rowvals = SUNSparseMatrix_IndexValues (Jac); 675
535 676 yout = NVecToCol (yy, m_num);
536 for (octave_f77_int_type i = 0; i < m_num + 1; i++) 677 ypout = NVecToCol (yyp, m_num);
537 colptrs[i] = to_f77_int (jac.cidx (i)); 678 tout(j) = tsol;
538 679
539 double *d = SUNSparseMatrix_Data (Jac); 680 for (octave_idx_type i = 0; i < m_num; i++)
540 for (octave_f77_int_type i = 0; i < to_f77_int (jac.nnz ()); i++) 681 output.elem (j, i) = yout.elem (i);
541 { 682
542 rowvals[i] = to_f77_int (jac.ridx (i)); 683 if (haveoutputfcn)
543 d[i] = jac.data (i); 684 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol,
544 } 685 tend, outputsel, string);
545 } 686
546 # endif 687 if (haveeventfunction)
547 688 status = IDA::event (event_fcn, te, ye, ie, tout(j), yout,
548 ColumnVector 689 string, ypout, oldval, oldisterminal,
549 IDA::NVecToCol (N_Vector& v, octave_f77_int_type n) 690 olddir, j, temp, tout(j-1), yold,
550 { 691 num_event_args);
551 ColumnVector data (n); 692
552 realtype *punt = nv_data_s (v); 693 // If integration is stopped, return only the reached steps
553 694 if (status == 1)
554 for (octave_f77_int_type i = 0; i < n; i++) 695 {
555 data(i) = punt[i]; 696 output.resize (j + 1, m_num);
556 697 tout.resize (j + 1);
557 return data; 698 }
558 } 699
559 700 }
560 N_Vector 701 }
561 IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n) 702 else // numel (tspan) == 2
562 { 703 {
563 N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT); 704 // First output value
564 705 tout.resize (1);
565 realtype *punt = nv_data_s (v); 706 tout(0) = tsol;
566 707 output.resize (1, m_num);
567 for (octave_f77_int_type i = 0; i < n; i++) 708
568 punt[i] = data(i); 709 for (octave_idx_type i = 0; i < m_num; i++)
569 710 output.elem (0, i) = y.elem (i);
570 return v; 711
571 } 712 bool posdirection = (tend > tsol);
572 713
573 void 714 //main loop
574 IDA::set_userdata (void) 715 while (((posdirection == 1 && tsol < tend)
575 { 716 || (posdirection == 0 && tsol > tend))
576 void *userdata = this; 717 && status == 0)
577 718 {
578 if (IDASetUserData (m_mem, userdata) != 0) 719 if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)
579 error ("User data not set"); 720 error ("IDASolve failed");
580 } 721
581 722 if (haverefine)
582 void 723 status = IDA::interpolate (cont, output, tout, refine, tend,
583 IDA::initialize (void) 724 haveoutputfcn, haveoutputsel,
584 { 725 output_fcn, outputsel,
585 m_num = to_f77_int (m_y0.numel ()); 726 haveeventfunction, event_fcn, te,
586 # if defined (HAVE_SUNDIALS_SUNCONTEXT) 727 ye, ie, oldval, oldisterminal,
587 if (SUNContext_Create (nullptr, &m_sunContext) < 0) 728 olddir, temp, yold,
588 error ("__ode15__: unable to create context for SUNDIALS"); 729 num_event_args);
589 m_mem = IDACreate (m_sunContext); 730
590 # else 731 ypout = NVecToCol (yyp, m_num);
591 m_mem = IDACreate (); 732 cont += 1;
592 # endif 733 output.resize (cont + 1, m_num); // This may be not efficient
593 734 tout.resize (cont + 1);
594 N_Vector yy = ColToNVec (m_y0, m_num); 735 tout(cont) = tsol;
595 736 yout = NVecToCol (yy, m_num);
596 N_Vector yyp = ColToNVec (m_yp0, m_num); 737
597 738 for (octave_idx_type i = 0; i < m_num; i++)
598 IDA::set_userdata (); 739 output.elem (cont, i) = yout.elem (i);
599 740
600 if (IDAInit (m_mem, IDA::resfun, m_t0, yy, yyp) != 0) 741 if (haveoutputfcn && ! haverefine && tout(cont) < tend)
601 error ("IDA not initialized"); 742 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol,
602 } 743 tend, outputsel, string);
603 744
604 void 745 if (haveeventfunction && ! haverefine && tout(cont) < tend)
605 IDA::set_tolerance (ColumnVector& abstol, realtype reltol) 746 status = IDA::event (event_fcn, te, ye, ie, tout(cont), yout,
606 { 747 string, ypout, oldval, oldisterminal,
607 N_Vector abs_tol = ColToNVec (abstol, m_num); 748 olddir, cont, temp, tout(cont-1), yold,
608 749 num_event_args);
609 if (IDASVtolerances (m_mem, reltol, abs_tol) != 0) 750 }
610 error ("IDA: Tolerance not set"); 751
611 752 if (status == 0)
612 N_VDestroy_Serial (abs_tol); 753 {
613 } 754 // Interpolate in tend
614 755 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
615 void 756
616 IDA::set_tolerance (realtype abstol, realtype reltol) 757 if (IDAGetDky (m_mem, tend, 0, dky) != 0)
617 { 758 error ("IDA failed to interpolate y");
618 if (IDASStolerances (m_mem, reltol, abstol) != 0) 759
619 error ("IDA: Tolerance not set"); 760 tout(cont) = tend;
620 } 761 yout = NVecToCol (dky, m_num);
621 762
622 octave_value_list 763 for (octave_idx_type i = 0; i < m_num; i++)
623 IDA::integrate (const octave_idx_type numt, const ColumnVector& tspan, 764 output.elem (cont, i) = yout.elem (i);
624 const ColumnVector& y, const ColumnVector& yp, 765
625 const int refine, bool haverefine, bool haveoutputfcn, 766 // Plot final value
626 const octave_value& output_fcn, bool haveoutputsel, 767 if (haveoutputfcn)
768 {
769 status = IDA::outputfun (output_fcn, haveoutputsel, yout,
770 tend, tend, outputsel, string);
771
772 // Events during last step
773 if (haveeventfunction)
774 status = IDA::event (event_fcn, te, ye, ie, tend, yout,
775 string, ypout, oldval, oldisterminal,
776 olddir, cont, temp, tout(cont-1),
777 yold, num_event_args);
778 }
779
780 N_VDestroy_Serial (dky);
781 }
782
783 // Cleanup plotter
784 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tend, tend,
785 outputsel, "done");
786
787 }
788
789 // Index of Events (ie) variable must use 1-based indexing
790 return ovl (tout, output, te, ye, ie + 1.0);
791 }
792
793 bool
794 IDA::event (const octave_value& event_fcn,
795 ColumnVector& te, Matrix& ye, ColumnVector& ie, realtype tsol,
796 const ColumnVector& y, const std::string& flag,
797 const ColumnVector& yp, ColumnVector& oldval,
798 ColumnVector& oldisterminal, ColumnVector& olddir,
799 octave_idx_type cont, octave_idx_type& temp, realtype told,
800 ColumnVector& yold,
801 const octave_idx_type num_event_args)
802 {
803 bool status = false;
804
805 octave_value_list args;
806 if (num_event_args == 2)
807 args = ovl (tsol, y);
808 else
809 args = ovl (tsol, y, yp);
810
811 // cont is the number of steps reached by the solver
812 // temp is the number of events registered
813
814 if (flag == "init")
815 {
816 octave_value_list output = feval (event_fcn, args, 3);
817 oldval = output(0).vector_value ();
818 oldisterminal = output(1).vector_value ();
819 olddir = output(2).vector_value ();
820 }
821 else if (flag == "")
822 {
823 ColumnVector index (0);
824 octave_value_list output = feval (event_fcn, args, 3);
825 ColumnVector val = output(0).vector_value ();
826 ColumnVector isterminal = output(1).vector_value ();
827 ColumnVector dir = output(2).vector_value ();
828
829 // Get the index of the changed values
830 for (octave_idx_type i = 0; i < val.numel (); i++)
831 {
832 // Check for sign change and whether a rising / falling edge
833 // either passes through zero or detaches from zero (bug #59063)
834 if ((dir(i) != -1
835 && ((val(i) >= 0 && oldval(i) < 0)
836 || (val(i) > 0 && oldval(i) <= 0))) // increasing
837 || (dir(i) != 1
838 && ((val(i) <= 0 && oldval(i) > 0)
839 || (val(i) < 0 && oldval(i) >= 0)))) // decreasing
840 {
841 index.resize (index.numel () + 1);
842 index (index.numel () - 1) = i;
843 }
844 }
845
846 if (cont == 1 && index.numel () > 0) // Events in first step
847 {
848 temp = 1; // register only the first event
849 te.resize (1);
850 ye.resize (1, m_num);
851 ie.resize (1);
852
853 // Linear interpolation
854 ie(0) = index(0);
855 te(0) = tsol - val (index(0)) * (tsol - told)
856 / (val (index(0)) - oldval (index(0)));
857
858 ColumnVector ytemp
859 = y - ((tsol - te(0)) * (y - yold) / (tsol - told));
860
861 for (octave_idx_type i = 0; i < m_num; i++)
862 ye.elem (0, i) = ytemp.elem (i);
863
864 }
865 else if (index.numel () > 0)
866 // Not first step: register all events and test
867 // if stop integration or not
868 {
869 te.resize (temp + index.numel ());
870 ye.resize (temp + index.numel (), m_num);
871 ie.resize (temp + index.numel ());
872
873 for (octave_idx_type i = 0; i < index.numel (); i++)
874 {
875
876 if (isterminal (index(i)) == 1)
877 status = 1; // Stop integration
878
879 // Linear interpolation
880 ie(temp+i) = index(i);
881 te(temp+i) = tsol - val(index(i)) * (tsol - told)
882 / (val(index(i)) - oldval(index(i)));
883
884 ColumnVector ytemp
885 = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told);
886
887 for (octave_idx_type j = 0; j < m_num; j++)
888 ye.elem (temp + i, j) = ytemp.elem (j);
889
890 }
891
892 temp += index.numel ();
893 }
894
895 // Update variables
896 yold = y;
897 told = tsol;
898 olddir = dir;
899 oldval = val;
900 oldisterminal = isterminal;
901 }
902
903 return status;
904 }
905
906 bool
907 IDA::interpolate (octave_idx_type& cont, Matrix& output, ColumnVector& tout,
908 int refine, realtype tend, bool haveoutputfcn,
909 bool haveoutputsel, const octave_value& output_fcn,
627 ColumnVector& outputsel, bool haveeventfunction, 910 ColumnVector& outputsel, bool haveeventfunction,
628 const octave_value& event_fcn, 911 const octave_value& event_fcn, ColumnVector& te,
912 Matrix& ye, ColumnVector& ie, ColumnVector& oldval,
913 ColumnVector& oldisterminal, ColumnVector& olddir,
914 octave_idx_type& temp, ColumnVector& yold,
629 const octave_idx_type num_event_args) 915 const octave_idx_type num_event_args)
630 { 916 {
631 // Set up output 917 realtype h = 0, tcur = 0;
632 ColumnVector tout, yout (m_num), ypout (m_num), ysel (outputsel.numel ()); 918 bool status = false;
633 ColumnVector ie, te, oldval, oldisterminal, olddir; 919
634 Matrix output, ye; 920 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
635 octave_idx_type cont = 0, temp = 0; 921
636 bool status = false; 922 N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
637 std::string string = ""; 923
638 ColumnVector yold = y; 924 ColumnVector yout (m_num);
639 925 ColumnVector ypout (m_num);
640 realtype tsol = tspan(0); 926 std::string string = "";
641 realtype tend = tspan(numt-1); 927
642 928 if (IDAGetLastStep (m_mem, &h) != 0)
643 N_Vector yyp = ColToNVec (yp, m_num); 929 error ("IDA failed to return last step");
644 930
645 N_Vector yy = ColToNVec (y, m_num); 931 if (IDAGetCurrentTime (m_mem, &tcur) != 0)
646 932 error ("IDA failed to return the current time");
647 // Initialize OutputFcn 933
648 if (haveoutputfcn) 934 realtype tin = tcur - h;
649 status = IDA::outputfun (output_fcn, haveoutputsel, y, 935
650 tsol, tend, outputsel, "init"); 936 realtype step = h / refine;
651 937
652 // Initialize Events 938 for (octave_idx_type i = 1;
653 if (haveeventfunction) 939 i < refine && tin + step * i < tend && status == 0;
654 status = IDA::event (event_fcn, te, ye, ie, tsol, y, 940 i++)
655 "init", yp, oldval, oldisterminal, 941 {
656 olddir, cont, temp, tsol, yold, num_event_args); 942 if (IDAGetDky (m_mem, tin + step*i, 0, dky) != 0)
657 943 error ("IDA failed to interpolate y");
658 if (numt > 2) 944
659 { 945 if (IDAGetDky (m_mem, tin + step*i, 1, dkyp) != 0)
660 // First output value 946 error ("IDA failed to interpolate yp");
661 tout.resize (numt); 947
662 tout(0) = tsol; 948 cont += 1;
663 output.resize (numt, m_num); 949 output.resize (cont + 1, m_num);
664 950 tout.resize (cont + 1);
665 for (octave_idx_type i = 0; i < m_num; i++) 951
666 output.elem (0, i) = y.elem (i); 952 tout(cont) = tin + step * i;
667 953 yout = NVecToCol (dky, m_num);
668 //Main loop 954 ypout = NVecToCol (dkyp, m_num);
669 for (octave_idx_type j = 1; j < numt && status == 0; j++) 955
670 { 956 for (octave_idx_type j = 0; j < m_num; j++)
671 // IDANORMAL already interpolates tspan(j) 957 output.elem (cont, j) = yout.elem (j);
672 958
673 if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0) 959 if (haveoutputfcn)
674 error ("IDASolve failed"); 960 status = IDA::outputfun (output_fcn, haveoutputsel, yout,
675 961 tout(cont), tend, outputsel, "");
676 yout = NVecToCol (yy, m_num); 962
677 ypout = NVecToCol (yyp, m_num); 963 if (haveeventfunction)
678 tout(j) = tsol; 964 status = IDA::event (event_fcn, te, ye, ie, tout(cont),
679 965 yout, string, ypout, oldval,
680 for (octave_idx_type i = 0; i < m_num; i++) 966 oldisterminal, olddir, cont, temp,
681 output.elem (j, i) = yout.elem (i); 967 tout(cont-1), yold, num_event_args);
682 968 }
683 if (haveoutputfcn) 969
684 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, 970 N_VDestroy_Serial (dky);
685 tend, outputsel, string); 971
686 972 return status;
687 if (haveeventfunction) 973 }
688 status = IDA::event (event_fcn, te, ye, ie, tout(j), yout, 974
689 string, ypout, oldval, oldisterminal, 975 bool
690 olddir, j, temp, tout(j-1), yold, 976 IDA::outputfun (const octave_value& output_fcn, bool haveoutputsel,
691 num_event_args); 977 const ColumnVector& yout, realtype tsol,
692 978 realtype tend, ColumnVector& outputsel,
693 // If integration is stopped, return only the reached steps 979 const std::string& flag)
694 if (status == 1) 980 {
695 { 981 bool status = false;
696 output.resize (j + 1, m_num); 982
697 tout.resize (j + 1); 983 octave_value_list output;
698 } 984 output(2) = flag;
699 985
700 } 986 ColumnVector ysel (outputsel.numel ());
701 } 987 if (haveoutputsel)
702 else // numel (tspan) == 2 988 {
703 { 989 for (octave_idx_type i = 0; i < outputsel.numel (); i++)
704 // First output value 990 ysel(i) = yout(outputsel(i));
705 tout.resize (1); 991
706 tout(0) = tsol; 992 output(1) = ysel;
707 output.resize (1, m_num); 993 }
708 994 else
709 for (octave_idx_type i = 0; i < m_num; i++) 995 output(1) = yout;
710 output.elem (0, i) = y.elem (i); 996
711 997 if (flag == "init")
712 bool posdirection = (tend > tsol); 998 {
713 999 ColumnVector toutput (2);
714 //main loop 1000 toutput(0) = tsol;
715 while (((posdirection == 1 && tsol < tend) 1001 toutput(1) = tend;
716 || (posdirection == 0 && tsol > tend)) 1002 output(0) = toutput;
717 && status == 0) 1003
718 { 1004 feval (output_fcn, output, 0);
719 if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0) 1005 }
720 error ("IDASolve failed"); 1006 else if (flag == "")
721 1007 {
722 if (haverefine) 1008 output(0) = tsol;
723 status = IDA::interpolate (cont, output, tout, refine, tend, 1009 octave_value_list val = feval (output_fcn, output, 1);
724 haveoutputfcn, haveoutputsel, 1010 status = val(0).bool_value ();
725 output_fcn, outputsel, 1011 }
726 haveeventfunction, event_fcn, te, 1012 else
727 ye, ie, oldval, oldisterminal, 1013 {
728 olddir, temp, yold, 1014 // Cleanup plotter
729 num_event_args); 1015 output(0) = tend;
730 1016 feval (output_fcn, output, 0);
731 ypout = NVecToCol (yyp, m_num); 1017 }
732 cont += 1; 1018
733 output.resize (cont + 1, m_num); // This may be not efficient 1019 return status;
734 tout.resize (cont + 1); 1020 }
735 tout(cont) = tsol; 1021
736 yout = NVecToCol (yy, m_num); 1022 void
737 1023 IDA::set_maxstep (realtype maxstep)
738 for (octave_idx_type i = 0; i < m_num; i++) 1024 {
739 output.elem (cont, i) = yout.elem (i); 1025 if (IDASetMaxStep (m_mem, maxstep) != 0)
740 1026 error ("IDA: Max Step not set");
741 if (haveoutputfcn && ! haverefine && tout(cont) < tend) 1027 }
742 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, 1028
743 tend, outputsel, string); 1029 void
744 1030 IDA::set_initialstep (realtype initialstep)
745 if (haveeventfunction && ! haverefine && tout(cont) < tend) 1031 {
746 status = IDA::event (event_fcn, te, ye, ie, tout(cont), yout, 1032 if (IDASetInitStep (m_mem, initialstep) != 0)
747 string, ypout, oldval, oldisterminal, 1033 error ("IDA: Initial Step not set");
748 olddir, cont, temp, tout(cont-1), yold, 1034 }
749 num_event_args); 1035
750 } 1036 void
751 1037 IDA::set_maxorder (int maxorder)
752 if (status == 0) 1038 {
753 { 1039 if (IDASetMaxOrd (m_mem, maxorder) != 0)
754 // Interpolate in tend 1040 error ("IDA: Max Order not set");
755 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); 1041 }
756 1042
757 if (IDAGetDky (m_mem, tend, 0, dky) != 0) 1043 void
758 error ("IDA failed to interpolate y"); 1044 IDA::print_stat (void)
759 1045 {
760 tout(cont) = tend; 1046 long int nsteps = 0, netfails = 0, nrevals = 0;
761 yout = NVecToCol (dky, m_num); 1047
762 1048 if (IDAGetNumSteps (m_mem, &nsteps) != 0)
763 for (octave_idx_type i = 0; i < m_num; i++) 1049 error ("IDA failed to return the number of internal steps");
764 output.elem (cont, i) = yout.elem (i); 1050
765 1051 if (IDAGetNumErrTestFails (m_mem, &netfails) != 0)
766 // Plot final value 1052 error ("IDA failed to return the number of internal errors");
767 if (haveoutputfcn) 1053
768 { 1054 if (IDAGetNumResEvals (m_mem, &nrevals) != 0)
769 status = IDA::outputfun (output_fcn, haveoutputsel, yout, 1055 error ("IDA failed to return the number of residual evaluations");
770 tend, tend, outputsel, string); 1056
771 1057 octave_stdout << nsteps << " successful steps\n";
772 // Events during last step 1058 octave_stdout << netfails << " failed attempts\n";
773 if (haveeventfunction) 1059 octave_stdout << nrevals << " function evaluations\n";
774 status = IDA::event (event_fcn, te, ye, ie, tend, yout, 1060 // octave_stdout << " partial derivatives\n";
775 string, ypout, oldval, oldisterminal, 1061 // octave_stdout << " LU decompositions\n";
776 olddir, cont, temp, tout(cont-1), 1062 // octave_stdout << " solutions of linear systems\n";
777 yold, num_event_args); 1063 }
778 } 1064
779 1065 static ColumnVector
780 N_VDestroy_Serial (dky); 1066 ida_user_function (const ColumnVector& x, const ColumnVector& xdot,
781 } 1067 double t, const octave_value& ida_fc)
782 1068 {
783 // Cleanup plotter 1069 octave_value_list tmp;
784 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tend, tend, 1070
785 outputsel, "done"); 1071 try
786 1072 {
787 } 1073 tmp = feval (ida_fc, ovl (t, x, xdot), 1);
788 1074 }
789 // Index of Events (ie) variable must use 1-based indexing 1075 catch (execution_exception& ee)
790 return ovl (tout, output, te, ye, ie + 1.0); 1076 {
791 } 1077 err_user_supplied_eval (ee, "__ode15__");
792 1078 }
793 bool 1079
794 IDA::event (const octave_value& event_fcn, 1080 return tmp(0).vector_value ();
795 ColumnVector& te, Matrix& ye, ColumnVector& ie, realtype tsol, 1081 }
796 const ColumnVector& y, const std::string& flag, 1082
797 const ColumnVector& yp, ColumnVector& oldval, 1083 static Matrix
798 ColumnVector& oldisterminal, ColumnVector& olddir, 1084 ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot,
799 octave_idx_type cont, octave_idx_type& temp, realtype told, 1085 double t, double cj, const octave_value& ida_jc)
800 ColumnVector& yold, 1086 {
801 const octave_idx_type num_event_args) 1087 octave_value_list tmp;
802 { 1088
803 bool status = false; 1089 try
804 1090 {
805 octave_value_list args; 1091 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
806 if (num_event_args == 2) 1092 }
807 args = ovl (tsol, y); 1093 catch (execution_exception& ee)
808 else 1094 {
809 args = ovl (tsol, y, yp); 1095 err_user_supplied_eval (ee, "__ode15__");
810 1096 }
811 // cont is the number of steps reached by the solver 1097
812 // temp is the number of events registered 1098 return tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
813 1099 }
814 if (flag == "init") 1100
815 { 1101 static SparseMatrix
816 octave_value_list output = feval (event_fcn, args, 3); 1102 ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot,
817 oldval = output(0).vector_value (); 1103 double t, double cj, const octave_value& ida_jc)
818 oldisterminal = output(1).vector_value (); 1104 {
819 olddir = output(2).vector_value (); 1105 octave_value_list tmp;
820 } 1106
821 else if (flag == "") 1107 try
822 { 1108 {
823 ColumnVector index (0); 1109 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
824 octave_value_list output = feval (event_fcn, args, 3); 1110 }
825 ColumnVector val = output(0).vector_value (); 1111 catch (execution_exception& ee)
826 ColumnVector isterminal = output(1).vector_value (); 1112 {
827 ColumnVector dir = output(2).vector_value (); 1113 err_user_supplied_eval (ee, "__ode15__");
828 1114 }
829 // Get the index of the changed values 1115
830 for (octave_idx_type i = 0; i < val.numel (); i++) 1116 return tmp(0).sparse_matrix_value () + cj * tmp(1).sparse_matrix_value ();
831 { 1117 }
832 // Check for sign change and whether a rising / falling edge 1118
833 // either passes through zero or detaches from zero (bug #59063) 1119 static Matrix
834 if ((dir(i) != -1 1120 ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj)
835 && ((val(i) >= 0 && oldval(i) < 0) 1121 {
836 || (val(i) > 0 && oldval(i) <= 0))) // increasing 1122 return (*dfdy) + cj * (*dfdyp);
837 || (dir(i) != 1 1123 }
838 && ((val(i) <= 0 && oldval(i) > 0) 1124
839 || (val(i) < 0 && oldval(i) >= 0)))) // decreasing 1125 static SparseMatrix
840 { 1126 ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp,
841 index.resize (index.numel () + 1); 1127 double cj)
842 index (index.numel () - 1) = i; 1128 {
843 } 1129 return (*spdfdy) + cj * (*spdfdyp);
844 } 1130 }
845 1131
846 if (cont == 1 && index.numel () > 0) // Events in first step 1132 static octave_value_list
847 { 1133 do_ode15 (const octave_value& ida_fcn,
848 temp = 1; // register only the first event 1134 const ColumnVector& tspan,
849 te.resize (1); 1135 const octave_idx_type numt,
850 ye.resize (1, m_num); 1136 const realtype t0,
851 ie.resize (1); 1137 const ColumnVector& y0,
852 1138 const ColumnVector& yp0,
853 // Linear interpolation 1139 const octave_scalar_map& options,
854 ie(0) = index(0); 1140 const octave_idx_type num_event_args)
855 te(0) = tsol - val (index(0)) * (tsol - told) 1141 {
856 / (val (index(0)) - oldval (index(0))); 1142 octave_value_list retval;
857 1143
858 ColumnVector ytemp 1144 // Create object
859 = y - ((tsol - te(0)) * (y - yold) / (tsol - told)); 1145 IDA dae (t0, y0, yp0, ida_fcn, ida_user_function);
860 1146
861 for (octave_idx_type i = 0; i < m_num; i++) 1147 // Set Jacobian
862 ye.elem (0, i) = ytemp.elem (i); 1148 bool havejac = options.getfield ("havejac").bool_value ();
863 1149
864 } 1150 bool havejacsparse = options.getfield ("havejacsparse").bool_value ();
865 else if (index.numel () > 0) 1151
866 // Not first step: register all events and test 1152 bool havejacfcn = options.getfield ("havejacfcn").bool_value ();
867 // if stop integration or not 1153
868 { 1154 Matrix ida_dfdy, ida_dfdyp;
869 te.resize (temp + index.numel ()); 1155 SparseMatrix ida_spdfdy, ida_spdfdyp;
870 ye.resize (temp + index.numel (), m_num); 1156
871 ie.resize (temp + index.numel ()); 1157 if (havejac)
872 1158 {
873 for (octave_idx_type i = 0; i < index.numel (); i++) 1159 if (havejacfcn)
874 { 1160 {
875 1161 octave_value ida_jac = options.getfield ("Jacobian");
876 if (isterminal (index(i)) == 1) 1162
877 status = 1; // Stop integration 1163 if (havejacsparse)
878 1164 dae.set_jacobian (ida_jac, ida_sparse_jac);
879 // Linear interpolation 1165 else
880 ie(temp+i) = index(i); 1166 dae.set_jacobian (ida_jac, ida_dense_jac);
881 te(temp+i) = tsol - val(index(i)) * (tsol - told) 1167 }
882 / (val(index(i)) - oldval(index(i))); 1168 else
883 1169 {
884 ColumnVector ytemp 1170 Cell jaccell = options.getfield ("Jacobian").cell_value ();
885 = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told); 1171
886 1172 if (havejacsparse)
887 for (octave_idx_type j = 0; j < m_num; j++) 1173 {
888 ye.elem (temp + i, j) = ytemp.elem (j); 1174 ida_spdfdy = jaccell(0).sparse_matrix_value ();
889 1175 ida_spdfdyp = jaccell(1).sparse_matrix_value ();
890 } 1176
891 1177 dae.set_jacobian (&ida_spdfdy, &ida_spdfdyp,
892 temp += index.numel (); 1178 ida_sparse_cell_jac);
893 } 1179 }
894 1180 else
895 // Update variables 1181 {
896 yold = y; 1182 ida_dfdy = jaccell(0).matrix_value ();
897 told = tsol; 1183 ida_dfdyp = jaccell(1).matrix_value ();
898 olddir = dir; 1184
899 oldval = val; 1185 dae.set_jacobian (&ida_dfdy, &ida_dfdyp, ida_dense_cell_jac);
900 oldisterminal = isterminal; 1186 }
901 } 1187 }
902 1188 }
903 return status; 1189
904 } 1190 // Initialize IDA
905 1191 dae.initialize ();
906 bool 1192
907 IDA::interpolate (octave_idx_type& cont, Matrix& output, ColumnVector& tout, 1193 // Set tolerances
908 int refine, realtype tend, bool haveoutputfcn, 1194 realtype rel_tol = options.getfield ("RelTol").double_value ();
909 bool haveoutputsel, const octave_value& output_fcn, 1195
910 ColumnVector& outputsel, bool haveeventfunction, 1196 bool haveabstolvec = options.getfield ("haveabstolvec").bool_value ();
911 const octave_value& event_fcn, ColumnVector& te, 1197
912 Matrix& ye, ColumnVector& ie, ColumnVector& oldval, 1198 if (haveabstolvec)
913 ColumnVector& oldisterminal, ColumnVector& olddir, 1199 {
914 octave_idx_type& temp, ColumnVector& yold, 1200 ColumnVector abs_tol = options.getfield ("AbsTol").vector_value ();
915 const octave_idx_type num_event_args) 1201
916 { 1202 dae.set_tolerance (abs_tol, rel_tol);
917 realtype h = 0, tcur = 0; 1203 }
918 bool status = false; 1204 else
919 1205 {
920 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); 1206 realtype abs_tol = options.getfield ("AbsTol").double_value ();
921 1207
922 N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT); 1208 dae.set_tolerance (abs_tol, rel_tol);
923 1209 }
924 ColumnVector yout (m_num); 1210
925 ColumnVector ypout (m_num); 1211 //Set max step
926 std::string string = ""; 1212 realtype maxstep = options.getfield ("MaxStep").double_value ();
927 1213
928 if (IDAGetLastStep (m_mem, &h) != 0) 1214 dae.set_maxstep (maxstep);
929 error ("IDA failed to return last step"); 1215
930 1216 //Set initial step
931 if (IDAGetCurrentTime (m_mem, &tcur) != 0) 1217 if (! options.getfield ("InitialStep").isempty ())
932 error ("IDA failed to return the current time"); 1218 {
933 1219 realtype initialstep = options.getfield ("InitialStep").double_value ();
934 realtype tin = tcur - h; 1220
935 1221 dae.set_initialstep (initialstep);
936 realtype step = h / refine; 1222 }
937 1223
938 for (octave_idx_type i = 1; 1224 //Set max order FIXME: it doesn't work
939 i < refine && tin + step * i < tend && status == 0; 1225 int maxorder = options.getfield ("MaxOrder").int_value ();
940 i++) 1226
941 { 1227 dae.set_maxorder (maxorder);
942 if (IDAGetDky (m_mem, tin + step*i, 0, dky) != 0) 1228
943 error ("IDA failed to interpolate y"); 1229 //Set Refine
944 1230 const int refine = options.getfield ("Refine").int_value ();
945 if (IDAGetDky (m_mem, tin + step*i, 1, dkyp) != 0) 1231
946 error ("IDA failed to interpolate yp"); 1232 bool haverefine = (refine > 1);
947 1233
948 cont += 1; 1234 octave_value output_fcn;
949 output.resize (cont + 1, m_num); 1235 ColumnVector outputsel;
950 tout.resize (cont + 1); 1236
951 1237 // OutputFcn
952 tout(cont) = tin + step * i; 1238 bool haveoutputfunction
953 yout = NVecToCol (dky, m_num); 1239 = options.getfield ("haveoutputfunction").bool_value ();
954 ypout = NVecToCol (dkyp, m_num); 1240
955 1241 if (haveoutputfunction)
956 for (octave_idx_type j = 0; j < m_num; j++) 1242 output_fcn = options.getfield ("OutputFcn");
957 output.elem (cont, j) = yout.elem (j); 1243
958 1244 // OutputSel
959 if (haveoutputfcn) 1245 bool haveoutputsel = options.getfield ("haveoutputselection").bool_value ();
960 status = IDA::outputfun (output_fcn, haveoutputsel, yout, 1246
961 tout(cont), tend, outputsel, ""); 1247 if (haveoutputsel)
962 1248 outputsel = options.getfield ("OutputSel").vector_value ();
963 if (haveeventfunction) 1249
964 status = IDA::event (event_fcn, te, ye, ie, tout(cont), 1250 octave_value event_fcn;
965 yout, string, ypout, oldval, 1251
966 oldisterminal, olddir, cont, temp, 1252 // Events
967 tout(cont-1), yold, num_event_args); 1253 bool haveeventfunction
968 } 1254 = options.getfield ("haveeventfunction").bool_value ();
969 1255
970 N_VDestroy_Serial (dky); 1256 if (haveeventfunction)
971 1257 event_fcn = options.getfield ("Events");
972 return status; 1258
973 } 1259 // Set up linear solver
974 1260 dae.set_up (y0);
975 bool 1261
976 IDA::outputfun (const octave_value& output_fcn, bool haveoutputsel, 1262 // Integrate
977 const ColumnVector& yout, realtype tsol, 1263 retval = dae.integrate (numt, tspan, y0, yp0, refine,
978 realtype tend, ColumnVector& outputsel, 1264 haverefine, haveoutputfunction,
979 const std::string& flag) 1265 output_fcn, haveoutputsel, outputsel,
980 { 1266 haveeventfunction, event_fcn, num_event_args);
981 bool status = false; 1267
982 1268 // Statistics
983 octave_value_list output; 1269 bool havestats = options.getfield ("havestats").bool_value ();
984 output(2) = flag; 1270
985 1271 if (havestats)
986 ColumnVector ysel (outputsel.numel ()); 1272 dae.print_stat ();
987 if (haveoutputsel) 1273
988 { 1274 return retval;
989 for (octave_idx_type i = 0; i < outputsel.numel (); i++) 1275 }
990 ysel(i) = yout(outputsel(i));
991
992 output(1) = ysel;
993 }
994 else
995 output(1) = yout;
996
997 if (flag == "init")
998 {
999 ColumnVector toutput (2);
1000 toutput(0) = tsol;
1001 toutput(1) = tend;
1002 output(0) = toutput;
1003
1004 feval (output_fcn, output, 0);
1005 }
1006 else if (flag == "")
1007 {
1008 output(0) = tsol;
1009 octave_value_list val = feval (output_fcn, output, 1);
1010 status = val(0).bool_value ();
1011 }
1012 else
1013 {
1014 // Cleanup plotter
1015 output(0) = tend;
1016 feval (output_fcn, output, 0);
1017 }
1018
1019 return status;
1020 }
1021
1022 void
1023 IDA::set_maxstep (realtype maxstep)
1024 {
1025 if (IDASetMaxStep (m_mem, maxstep) != 0)
1026 error ("IDA: Max Step not set");
1027 }
1028
1029 void
1030 IDA::set_initialstep (realtype initialstep)
1031 {
1032 if (IDASetInitStep (m_mem, initialstep) != 0)
1033 error ("IDA: Initial Step not set");
1034 }
1035
1036 void
1037 IDA::set_maxorder (int maxorder)
1038 {
1039 if (IDASetMaxOrd (m_mem, maxorder) != 0)
1040 error ("IDA: Max Order not set");
1041 }
1042
1043 void
1044 IDA::print_stat (void)
1045 {
1046 long int nsteps = 0, netfails = 0, nrevals = 0;
1047
1048 if (IDAGetNumSteps (m_mem, &nsteps) != 0)
1049 error ("IDA failed to return the number of internal steps");
1050
1051 if (IDAGetNumErrTestFails (m_mem, &netfails) != 0)
1052 error ("IDA failed to return the number of internal errors");
1053
1054 if (IDAGetNumResEvals (m_mem, &nrevals) != 0)
1055 error ("IDA failed to return the number of residual evaluations");
1056
1057 octave_stdout << nsteps << " successful steps\n";
1058 octave_stdout << netfails << " failed attempts\n";
1059 octave_stdout << nrevals << " function evaluations\n";
1060 // octave_stdout << " partial derivatives\n";
1061 // octave_stdout << " LU decompositions\n";
1062 // octave_stdout << " solutions of linear systems\n";
1063 }
1064
1065 static ColumnVector
1066 ida_user_function (const ColumnVector& x, const ColumnVector& xdot,
1067 double t, const octave_value& ida_fc)
1068 {
1069 octave_value_list tmp;
1070
1071 try
1072 {
1073 tmp = feval (ida_fc, ovl (t, x, xdot), 1);
1074 }
1075 catch (execution_exception& ee)
1076 {
1077 err_user_supplied_eval (ee, "__ode15__");
1078 }
1079
1080 return tmp(0).vector_value ();
1081 }
1082
1083 static Matrix
1084 ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot,
1085 double t, double cj, const octave_value& ida_jc)
1086 {
1087 octave_value_list tmp;
1088
1089 try
1090 {
1091 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
1092 }
1093 catch (execution_exception& ee)
1094 {
1095 err_user_supplied_eval (ee, "__ode15__");
1096 }
1097
1098 return tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
1099 }
1100
1101 static SparseMatrix
1102 ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot,
1103 double t, double cj, const octave_value& ida_jc)
1104 {
1105 octave_value_list tmp;
1106
1107 try
1108 {
1109 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
1110 }
1111 catch (execution_exception& ee)
1112 {
1113 err_user_supplied_eval (ee, "__ode15__");
1114 }
1115
1116 return tmp(0).sparse_matrix_value () + cj * tmp(1).sparse_matrix_value ();
1117 }
1118
1119 static Matrix
1120 ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj)
1121 {
1122 return (*dfdy) + cj * (*dfdyp);
1123 }
1124
1125 static SparseMatrix
1126 ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp,
1127 double cj)
1128 {
1129 return (*spdfdy) + cj * (*spdfdyp);
1130 }
1131
1132 static octave_value_list
1133 do_ode15 (const octave_value& ida_fcn,
1134 const ColumnVector& tspan,
1135 const octave_idx_type numt,
1136 const realtype t0,
1137 const ColumnVector& y0,
1138 const ColumnVector& yp0,
1139 const octave_scalar_map& options,
1140 const octave_idx_type num_event_args)
1141 {
1142 octave_value_list retval;
1143
1144 // Create object
1145 IDA dae (t0, y0, yp0, ida_fcn, ida_user_function);
1146
1147 // Set Jacobian
1148 bool havejac = options.getfield ("havejac").bool_value ();
1149
1150 bool havejacsparse = options.getfield ("havejacsparse").bool_value ();
1151
1152 bool havejacfcn = options.getfield ("havejacfcn").bool_value ();
1153
1154 Matrix ida_dfdy, ida_dfdyp;
1155 SparseMatrix ida_spdfdy, ida_spdfdyp;
1156
1157 if (havejac)
1158 {
1159 if (havejacfcn)
1160 {
1161 octave_value ida_jac = options.getfield ("Jacobian");
1162
1163 if (havejacsparse)
1164 dae.set_jacobian (ida_jac, ida_sparse_jac);
1165 else
1166 dae.set_jacobian (ida_jac, ida_dense_jac);
1167 }
1168 else
1169 {
1170 Cell jaccell = options.getfield ("Jacobian").cell_value ();
1171
1172 if (havejacsparse)
1173 {
1174 ida_spdfdy = jaccell(0).sparse_matrix_value ();
1175 ida_spdfdyp = jaccell(1).sparse_matrix_value ();
1176
1177 dae.set_jacobian (&ida_spdfdy, &ida_spdfdyp,
1178 ida_sparse_cell_jac);
1179 }
1180 else
1181 {
1182 ida_dfdy = jaccell(0).matrix_value ();
1183 ida_dfdyp = jaccell(1).matrix_value ();
1184
1185 dae.set_jacobian (&ida_dfdy, &ida_dfdyp, ida_dense_cell_jac);
1186 }
1187 }
1188 }
1189
1190 // Initialize IDA
1191 dae.initialize ();
1192
1193 // Set tolerances
1194 realtype rel_tol = options.getfield ("RelTol").double_value ();
1195
1196 bool haveabstolvec = options.getfield ("haveabstolvec").bool_value ();
1197
1198 if (haveabstolvec)
1199 {
1200 ColumnVector abs_tol = options.getfield ("AbsTol").vector_value ();
1201
1202 dae.set_tolerance (abs_tol, rel_tol);
1203 }
1204 else
1205 {
1206 realtype abs_tol = options.getfield ("AbsTol").double_value ();
1207
1208 dae.set_tolerance (abs_tol, rel_tol);
1209 }
1210
1211 //Set max step
1212 realtype maxstep = options.getfield ("MaxStep").double_value ();
1213
1214 dae.set_maxstep (maxstep);
1215
1216 //Set initial step
1217 if (! options.getfield ("InitialStep").isempty ())
1218 {
1219 realtype initialstep = options.getfield ("InitialStep").double_value ();
1220
1221 dae.set_initialstep (initialstep);
1222 }
1223
1224 //Set max order FIXME: it doesn't work
1225 int maxorder = options.getfield ("MaxOrder").int_value ();
1226
1227 dae.set_maxorder (maxorder);
1228
1229 //Set Refine
1230 const int refine = options.getfield ("Refine").int_value ();
1231
1232 bool haverefine = (refine > 1);
1233
1234 octave_value output_fcn;
1235 ColumnVector outputsel;
1236
1237 // OutputFcn
1238 bool haveoutputfunction
1239 = options.getfield ("haveoutputfunction").bool_value ();
1240
1241 if (haveoutputfunction)
1242 output_fcn = options.getfield ("OutputFcn");
1243
1244 // OutputSel
1245 bool haveoutputsel = options.getfield ("haveoutputselection").bool_value ();
1246
1247 if (haveoutputsel)
1248 outputsel = options.getfield ("OutputSel").vector_value ();
1249
1250 octave_value event_fcn;
1251
1252 // Events
1253 bool haveeventfunction
1254 = options.getfield ("haveeventfunction").bool_value ();
1255
1256 if (haveeventfunction)
1257 event_fcn = options.getfield ("Events");
1258
1259 // Set up linear solver
1260 dae.set_up (y0);
1261
1262 // Integrate
1263 retval = dae.integrate (numt, tspan, y0, yp0, refine,
1264 haverefine, haveoutputfunction,
1265 output_fcn, haveoutputsel, outputsel,
1266 haveeventfunction, event_fcn, num_event_args);
1267
1268 // Statistics
1269 bool havestats = options.getfield ("havestats").bool_value ();
1270
1271 if (havestats)
1272 dae.print_stat ();
1273
1274 return retval;
1275 }
1276 1276
1277 #endif 1277 #endif
1278 1278
1279 DEFUN_DLD (__ode15__, args, , 1279 DEFUN_DLD (__ode15__, args, ,
1280 doc: /* -*- texinfo -*- 1280 doc: /* -*- texinfo -*-