GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
__ode15__.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2016-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include "dColVector.h"
31#include "dMatrix.h"
32#include "dSparse.h"
33#include "f77-fcn.h"
34#include "lo-utils.h"
35
36#include "Cell.h"
37#include "defun-dld.h"
38#include "error.h"
39#include "errwarn.h"
40#include "oct-map.h"
41#include "ov.h"
42#include "ovl.h"
43#include "pager.h"
44#include "parse.h"
45
46#if defined (HAVE_SUNDIALS)
47
48# if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H)
49# include <nvector/nvector_serial.h>
50# endif
51
52# if defined (HAVE_IDA_IDA_H)
53# include <ida/ida.h>
54# elif defined (HAVE_IDA_H)
55# include <ida.h>
56# endif
57# if defined (HAVE_IDA_IDA_DIRECT_H)
58# include <ida/ida_direct.h>
59# elif defined (HAVE_IDA_DIRECT_H)
60# include <ida_direct.h>
61# endif
62
63# if defined (HAVE_SUNLINSOL_SUNLINSOL_DENSE_H)
64# include <sunlinsol/sunlinsol_dense.h>
65# endif
66
67# if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
68# if defined (HAVE_KLU_H)
69# include <klu.h>
70# endif
71# if defined (HAVE_KLU_KLU_H)
72# include <klu/klu.h>
73# endif
74# if defined (HAVE_SUITESPARSE_KLU_H)
75# include <suitesparse/klu.h>
76# endif
77# if defined (HAVE_UFPARSE_KLU_H)
78# include <ufsparse/klu.h>
79# endif
80# include <sunlinsol/sunlinsol_klu.h>
81# endif
82
83#endif
84
85OCTAVE_NAMESPACE_BEGIN
86
87#if defined (HAVE_SUNDIALS)
88
89# if ! defined (HAVE_IDASETJACFN) && defined (HAVE_IDADLSSETJACFN)
90 static inline int
91 IDASetJacFn (void *ida_mem, IDADlsJacFn jac)
92 {
93 return IDADlsSetJacFn (ida_mem, jac);
94 }
95# endif
96
97# if ! defined (HAVE_IDASETLINEARSOLVER) && defined (HAVE_IDADLSSETLINEARSOLVER)
98 static inline int
99 IDASetLinearSolver (void *ida_mem, SUNLinearSolver LS, SUNMatrix A)
100 {
101 return IDADlsSetLinearSolver (ida_mem, LS, A);
102 }
103# endif
104
105# if ! defined (HAVE_SUNLINSOL_DENSE) && defined (HAVE_SUNDENSELINEARSOLVER)
106 static inline SUNLinearSolver
107 SUNLinSol_Dense (N_Vector y, SUNMatrix A)
108 {
109 return SUNDenseLinearSolver (y, A);
110 }
111# endif
112
113# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
114# if ! defined (HAVE_SUNLINSOL_KLU) && defined (HAVE_SUNKLU)
115 static inline SUNLinearSolver
116 SUNLinSol_KLU (N_Vector y, SUNMatrix A)
117 {
118 return SUNKLU (y, A);
119 }
120# endif
121# endif
122
123 static inline realtype *
124 nv_data_s (N_Vector& v)
125 {
126# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
127 // Disable warning from GCC about old-style casts in Sundials
128 // macro expansions. Do this in a function so that this
129 // diagnostic may still be enabled for the rest of the file.
130# pragma GCC diagnostic push
131# pragma GCC diagnostic ignored "-Wold-style-cast"
132# endif
133
134 return NV_DATA_S (v);
135
136# if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
137 // Restore prevailing warning state for remainder of the file.
138# pragma GCC diagnostic pop
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_havejacfun (false),
173 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fun (),
174 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
175 m_spdfdyp (nullptr), m_fun (nullptr), m_jacfun (nullptr),
176 m_jacspfun (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_havejacfun (false),
184 m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fun (ida_fcn),
185 m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
186 m_spdfdyp (nullptr), m_fun (daefun), m_jacfun (nullptr),
187 m_jacspfun (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
200 }
201
202 IDA&
203 set_jacobian (const octave_value& jac, DAEJacFuncDense j)
204 {
205 m_jacfun = j;
206 m_ida_jac = jac;
207 m_havejac = true;
208 m_havejacfun = 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_jacspfun = j;
218 m_ida_jac = jac;
219 m_havejac = true;
220 m_havejacfun = 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_havejacfun = 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_havejacfun = 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
262 static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
263# endif
264
265 void
266 set_up (const ColumnVector& y);
267
268 void
269 set_tolerance (ColumnVector& abstol, realtype reltol);
270
271 void
272 set_tolerance (realtype abstol, realtype reltol);
273
274 static int
275 resfun (realtype t, N_Vector yy, N_Vector yyp,
276 N_Vector rr, void *user_data);
277
278 void
279 resfun_impl (realtype t, N_Vector& yy,
280 N_Vector& yyp, N_Vector& rr);
281 static int
282 jacdense (realtype t, realtype cj, N_Vector yy, N_Vector yyp,
283 N_Vector, SUNMatrix JJ, void *user_data, 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_havejacfun;
362 bool m_havejacsparse;
363 void *m_mem;
364 octave_f77_int_type m_num;
365 octave_value m_ida_fun;
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_fun;
372 DAEJacFuncDense m_jacfun;
373 DAEJacFuncSparse m_jacspfun;
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_fun) (y, yp, t, m_ida_fun);
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
408# if defined (HAVE_SUNDIALS_SUNCONTEXT)
409# define OCTAVE_SUNCONTEXT , m_sunContext
410# else
411# define OCTAVE_SUNCONTEXT
412# endif
413
414 void
415 IDA::set_up (const ColumnVector& y)
416 {
417 N_Vector yy = ColToNVec (y, m_num);
418
419 if (m_havejacsparse)
420 {
421# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
422# if defined (HAVE_SUNSPARSEMATRIX_REALLOCATE)
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
425 // entries.
426 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT
427 OCTAVE_SUNCONTEXT);
428# else
429 octave_f77_int_type max_elems;
430 if (math::int_multiply_overflow (m_num, m_num, &max_elems))
431 error ("Unable to allocate memory for sparse Jacobian");
432
433 m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT
434 OCTAVE_SUNCONTEXT);
435# endif
436 if (! m_sunJacMatrix)
437 error ("Unable to create sparse Jacobian for Sundials");
438
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
464 OCTAVE_SUNCONTEXT);
465 if (! m_sunLinearSolver)
466 error ("Unable to create dense linear solver");
467
468 if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
469 error ("Unable to set dense linear solver");
470
471 if (m_havejac && IDASetJacFn (m_mem, IDA::jacdense) != 0)
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_havejacfun)
491 jac = (*m_jacfun) (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_havejacfun)
514 jac = (*m_jacspfun) (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);
590# else
591 m_mem = IDACreate ();
592# endif
593
594 N_Vector yy = ColToNVec (m_y0, m_num);
595
596 N_Vector yyp = ColToNVec (m_yp0, m_num);
597
598 IDA::set_userdata ();
599
600 if (IDAInit (m_mem, IDA::resfun, m_t0, yy, yyp) != 0)
601 error ("IDA not initialized");
602 }
603
604 void
605 IDA::set_tolerance (ColumnVector& abstol, realtype reltol)
606 {
607 N_Vector abs_tol = ColToNVec (abstol, m_num);
608
609 if (IDASVtolerances (m_mem, reltol, abs_tol) != 0)
610 error ("IDA: Tolerance not set");
611
612 N_VDestroy_Serial (abs_tol);
613 }
614
615 void
616 IDA::set_tolerance (realtype abstol, realtype reltol)
617 {
618 if (IDASStolerances (m_mem, reltol, abstol) != 0)
619 error ("IDA: Tolerance not set");
620 }
621
622 octave_value_list
623 IDA::integrate (const octave_idx_type numt, const ColumnVector& tspan,
624 const ColumnVector& y, const ColumnVector& yp,
625 const int refine, bool haverefine, bool haveoutputfcn,
626 const octave_value& output_fcn, bool haveoutputsel,
627 ColumnVector& outputsel, bool haveeventfunction,
628 const octave_value& event_fcn,
629 const octave_idx_type num_event_args)
630 {
631 // Set up output
632 ColumnVector tout, yout (m_num), ypout (m_num), ysel (outputsel.numel ());
633 ColumnVector ie, te, oldval, oldisterminal, olddir;
634 Matrix output, ye;
635 octave_idx_type cont = 0, temp = 0;
636 bool status = 0;
637 std::string string = "";
638 ColumnVector yold = y;
639
640 realtype tsol = tspan(0);
641 realtype tend = tspan(numt-1);
642
643 N_Vector yyp = ColToNVec (yp, m_num);
644
645 N_Vector yy = ColToNVec (y, m_num);
646
647 // Initialize OutputFcn
648 if (haveoutputfcn)
649 status = IDA::outputfun (output_fcn, haveoutputsel, y,
650 tsol, tend, outputsel, "init");
651
652 // Initialize Events
653 if (haveeventfunction)
654 status = IDA::event (event_fcn, te, ye, ie, tsol, y,
655 "init", yp, oldval, oldisterminal,
656 olddir, cont, temp, tsol, yold, num_event_args);
657
658 if (numt > 2)
659 {
660 // First output value
661 tout.resize (numt);
662 tout(0) = tsol;
663 output.resize (numt, m_num);
664
665 for (octave_idx_type i = 0; i < m_num; i++)
666 output.elem (0, i) = y.elem (i);
667
668 //Main loop
669 for (octave_idx_type j = 1; j < numt && status == 0; j++)
670 {
671 // IDANORMAL already interpolates tspan(j)
672
673 if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
674 error ("IDASolve failed");
675
676 yout = NVecToCol (yy, m_num);
677 ypout = NVecToCol (yyp, m_num);
678 tout(j) = tsol;
679
680 for (octave_idx_type i = 0; i < m_num; i++)
681 output.elem (j, i) = yout.elem (i);
682
683 if (haveoutputfcn)
684 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol,
685 tend, outputsel, string);
686
687 if (haveeventfunction)
688 status = IDA::event (event_fcn, te, ye, ie, tout(j), yout,
689 string, ypout, oldval, oldisterminal,
690 olddir, j, temp, tout(j-1), yold,
691 num_event_args);
692
693 // If integration is stopped, return only the reached steps
694 if (status == 1)
695 {
696 output.resize (j + 1, m_num);
697 tout.resize (j + 1);
698 }
699
700 }
701 }
702 else // numel (tspan) == 2
703 {
704 // First output value
705 tout.resize (1);
706 tout(0) = tsol;
707 output.resize (1, m_num);
708
709 for (octave_idx_type i = 0; i < m_num; i++)
710 output.elem (0, i) = y.elem (i);
711
712 bool posdirection = (tend > tsol);
713
714 //main loop
715 while (((posdirection == 1 && tsol < tend)
716 || (posdirection == 0 && tsol > tend))
717 && status == 0)
718 {
719 if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)
720 error ("IDASolve failed");
721
722 if (haverefine)
723 status = IDA::interpolate (cont, output, tout, refine, tend,
724 haveoutputfcn, haveoutputsel,
725 output_fcn, outputsel,
726 haveeventfunction, event_fcn, te,
727 ye, ie, oldval, oldisterminal,
728 olddir, temp, yold,
729 num_event_args);
730
731 ypout = NVecToCol (yyp, m_num);
732 cont += 1;
733 output.resize (cont + 1, m_num); // This may be not efficient
734 tout.resize (cont + 1);
735 tout(cont) = tsol;
736 yout = NVecToCol (yy, m_num);
737
738 for (octave_idx_type i = 0; i < m_num; i++)
739 output.elem (cont, i) = yout.elem (i);
740
741 if (haveoutputfcn && ! haverefine && tout(cont) < tend)
742 status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol,
743 tend, outputsel, string);
744
745 if (haveeventfunction && ! haverefine && tout(cont) < tend)
746 status = IDA::event (event_fcn, te, ye, ie, tout(cont), yout,
747 string, ypout, oldval, oldisterminal,
748 olddir, cont, temp, tout(cont-1), yold,
749 num_event_args);
750 }
751
752 if (status == 0)
753 {
754 // Interpolate in tend
755 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
756
757 if (IDAGetDky (m_mem, tend, 0, dky) != 0)
758 error ("IDA failed to interpolate y");
759
760 tout(cont) = tend;
761 yout = NVecToCol (dky, m_num);
762
763 for (octave_idx_type i = 0; i < m_num; i++)
764 output.elem (cont, i) = yout.elem (i);
765
766 // Plot final value
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 = 0;
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 if stop integration or not
867 {
868 te.resize (temp + index.numel ());
869 ye.resize (temp + index.numel (), m_num);
870 ie.resize (temp + index.numel ());
871
872 for (octave_idx_type i = 0; i < index.numel (); i++)
873 {
874
875 if (isterminal (index(i)) == 1)
876 status = 1; // Stop integration
877
878 // Linear interpolation
879 ie(temp+i) = index(i);
880 te(temp+i) = tsol - val(index(i)) * (tsol - told)
881 / (val(index(i)) - oldval(index(i)));
882
883 ColumnVector ytemp
884 = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told);
885
886 for (octave_idx_type j = 0; j < m_num; j++)
887 ye.elem (temp + i, j) = ytemp.elem (j);
888
889 }
890
891 temp += index.numel ();
892 }
893
894 // Update variables
895 yold = y;
896 told = tsol;
897 olddir = dir;
898 oldval = val;
899 oldisterminal = isterminal;
900 }
901
902 return status;
903 }
904
905 bool
906 IDA::interpolate (octave_idx_type& cont, Matrix& output, ColumnVector& tout,
907 int refine, realtype tend, bool haveoutputfcn,
908 bool haveoutputsel, const octave_value& output_fcn,
909 ColumnVector& outputsel, bool haveeventfunction,
910 const octave_value& event_fcn, ColumnVector& te,
911 Matrix& ye, ColumnVector& ie, ColumnVector& oldval,
912 ColumnVector& oldisterminal, ColumnVector& olddir,
913 octave_idx_type& temp, ColumnVector& yold,
914 const octave_idx_type num_event_args)
915 {
916 realtype h = 0, tcur = 0;
917 bool status = 0;
918
919 N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
920
921 N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
922
923 ColumnVector yout (m_num);
924 ColumnVector ypout (m_num);
925 std::string string = "";
926
927 if (IDAGetLastStep (m_mem, &h) != 0)
928 error ("IDA failed to return last step");
929
930 if (IDAGetCurrentTime (m_mem, &tcur) != 0)
931 error ("IDA failed to return the current time");
932
933 realtype tin = tcur - h;
934
935 realtype step = h / refine;
936
937 for (octave_idx_type i = 1;
938 i < refine && tin + step * i < tend && status == 0;
939 i++)
940 {
941 if (IDAGetDky (m_mem, tin + step*i, 0, dky) != 0)
942 error ("IDA failed to interpolate y");
943
944 if (IDAGetDky (m_mem, tin + step*i, 1, dkyp) != 0)
945 error ("IDA failed to interpolate yp");
946
947 cont += 1;
948 output.resize (cont + 1, m_num);
949 tout.resize (cont + 1);
950
951 tout(cont) = tin + step * i;
952 yout = NVecToCol (dky, m_num);
953 ypout = NVecToCol (dkyp, m_num);
954
955 for (octave_idx_type j = 0; j < m_num; j++)
956 output.elem (cont, j) = yout.elem (j);
957
958 if (haveoutputfcn)
959 status = IDA::outputfun (output_fcn, haveoutputsel, yout,
960 tout(cont), tend, outputsel, "");
961
962 if (haveeventfunction)
963 status = IDA::event (event_fcn, te, ye, ie, tout(cont),
964 yout, string, ypout, oldval,
965 oldisterminal, olddir, cont, temp,
966 tout(cont-1), yold, num_event_args);
967 }
968
969 N_VDestroy_Serial (dky);
970
971 return status;
972 }
973
974 bool
975 IDA::outputfun (const octave_value& output_fcn, bool haveoutputsel,
976 const ColumnVector& yout, realtype tsol,
977 realtype tend, ColumnVector& outputsel,
978 const std::string& flag)
979 {
980 bool status = 0;
981
982 octave_value_list output;
983 output(2) = flag;
984
985 ColumnVector ysel (outputsel.numel ());
986 if (haveoutputsel)
987 {
988 for (octave_idx_type i = 0; i < outputsel.numel (); i++)
989 ysel(i) = yout(outputsel(i));
990
991 output(1) = ysel;
992 }
993 else
994 output(1) = yout;
995
996 if (flag == "init")
997 {
998 ColumnVector toutput (2);
999 toutput(0) = tsol;
1000 toutput(1) = tend;
1001 output(0) = toutput;
1002
1003 feval (output_fcn, output, 0);
1004 }
1005 else if (flag == "")
1006 {
1007 output(0) = tsol;
1008 octave_value_list val = feval (output_fcn, output, 1);
1009 status = val(0).bool_value ();
1010 }
1011 else
1012 {
1013 // Cleanup plotter
1014 output(0) = tend;
1015 feval (output_fcn, output, 0);
1016 }
1017
1018 return status;
1019 }
1020
1021 void
1022 IDA::set_maxstep (realtype maxstep)
1023 {
1024 if (IDASetMaxStep (m_mem, maxstep) != 0)
1025 error ("IDA: Max Step not set");
1026 }
1027
1028 void
1029 IDA::set_initialstep (realtype initialstep)
1030 {
1031 if (IDASetInitStep (m_mem, initialstep) != 0)
1032 error ("IDA: Initial Step not set");
1033 }
1034
1035 void
1036 IDA::set_maxorder (int maxorder)
1037 {
1038 if (IDASetMaxOrd (m_mem, maxorder) != 0)
1039 error ("IDA: Max Order not set");
1040 }
1041
1042 void
1043 IDA::print_stat (void)
1044 {
1045 long int nsteps = 0, netfails = 0, nrevals = 0;
1046
1047 if (IDAGetNumSteps (m_mem, &nsteps) != 0)
1048 error ("IDA failed to return the number of internal steps");
1049
1050 if (IDAGetNumErrTestFails (m_mem, &netfails) != 0)
1051 error ("IDA failed to return the number of internal errors");
1052
1053 if (IDAGetNumResEvals (m_mem, &nrevals) != 0)
1054 error ("IDA failed to return the number of residual evaluations");
1055
1056 octave_stdout << nsteps << " successful steps\n";
1057 octave_stdout << netfails << " failed attempts\n";
1058 octave_stdout << nrevals << " function evaluations\n";
1059 // octave_stdout << " partial derivatives\n";
1060 // octave_stdout << " LU decompositions\n";
1061 // octave_stdout << " solutions of linear systems\n";
1062 }
1063
1064 static ColumnVector
1065 ida_user_function (const ColumnVector& x, const ColumnVector& xdot,
1066 double t, const octave_value& ida_fc)
1067 {
1068 octave_value_list tmp;
1069
1070 try
1071 {
1072 tmp = feval (ida_fc, ovl (t, x, xdot), 1);
1073 }
1074 catch (execution_exception& ee)
1075 {
1076 err_user_supplied_eval (ee, "__ode15__");
1077 }
1078
1079 return tmp(0).vector_value ();
1080 }
1081
1082 static Matrix
1083 ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot,
1084 double t, double cj, const octave_value& ida_jc)
1085 {
1086 octave_value_list tmp;
1087
1088 try
1089 {
1090 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
1091 }
1092 catch (execution_exception& ee)
1093 {
1094 err_user_supplied_eval (ee, "__ode15__");
1095 }
1096
1097 return tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
1098 }
1099
1100 static SparseMatrix
1101 ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot,
1102 double t, double cj, const octave_value& ida_jc)
1103 {
1104 octave_value_list tmp;
1105
1106 try
1107 {
1108 tmp = feval (ida_jc, ovl (t, x, xdot), 2);
1109 }
1110 catch (execution_exception& ee)
1111 {
1112 err_user_supplied_eval (ee, "__ode15__");
1113 }
1114
1115 return tmp(0).sparse_matrix_value () + cj * tmp(1).sparse_matrix_value ();
1116 }
1117
1118 static Matrix
1119 ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj)
1120 {
1121 return (*dfdy) + cj * (*dfdyp);
1122 }
1123
1124 static SparseMatrix
1125 ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp,
1126 double cj)
1127 {
1128 return (*spdfdy) + cj * (*spdfdyp);
1129 }
1130
1131 static octave_value_list
1132 do_ode15 (const octave_value& ida_fcn,
1133 const ColumnVector& tspan,
1134 const octave_idx_type numt,
1135 const realtype t0,
1136 const ColumnVector& y0,
1137 const ColumnVector& yp0,
1138 const octave_scalar_map& options,
1139 const octave_idx_type num_event_args)
1140 {
1141 octave_value_list retval;
1142
1143 // Create object
1144 IDA dae (t0, y0, yp0, ida_fcn, ida_user_function);
1145
1146 // Set Jacobian
1147 bool havejac = options.getfield ("havejac").bool_value ();
1148
1149 bool havejacsparse = options.getfield ("havejacsparse").bool_value ();
1150
1151 bool havejacfun = options.getfield ("havejacfun").bool_value ();
1152
1153 Matrix ida_dfdy, ida_dfdyp;
1154 SparseMatrix ida_spdfdy, ida_spdfdyp;
1155
1156 if (havejac)
1157 {
1158 if (havejacfun)
1159 {
1160 octave_value ida_jac = options.getfield ("Jacobian");
1161
1162 if (havejacsparse)
1163 dae.set_jacobian (ida_jac, ida_sparse_jac);
1164 else
1165 dae.set_jacobian (ida_jac, ida_dense_jac);
1166 }
1167 else
1168 {
1169 Cell jaccell = options.getfield ("Jacobian").cell_value ();
1170
1171 if (havejacsparse)
1172 {
1173 ida_spdfdy = jaccell(0).sparse_matrix_value ();
1174 ida_spdfdyp = jaccell(1).sparse_matrix_value ();
1175
1176 dae.set_jacobian (&ida_spdfdy, &ida_spdfdyp,
1177 ida_sparse_cell_jac);
1178 }
1179 else
1180 {
1181 ida_dfdy = jaccell(0).matrix_value ();
1182 ida_dfdyp = jaccell(1).matrix_value ();
1183
1184 dae.set_jacobian (&ida_dfdy, &ida_dfdyp, ida_dense_cell_jac);
1185 }
1186 }
1187 }
1188
1189 // Initialize IDA
1190 dae.initialize ();
1191
1192 // Set tolerances
1193 realtype rel_tol = options.getfield ("RelTol").double_value ();
1194
1195 bool haveabstolvec = options.getfield ("haveabstolvec").bool_value ();
1196
1197 if (haveabstolvec)
1198 {
1199 ColumnVector abs_tol = options.getfield ("AbsTol").vector_value ();
1200
1201 dae.set_tolerance (abs_tol, rel_tol);
1202 }
1203 else
1204 {
1205 realtype abs_tol = options.getfield ("AbsTol").double_value ();
1206
1207 dae.set_tolerance (abs_tol, rel_tol);
1208 }
1209
1210 //Set max step
1211 realtype maxstep = options.getfield ("MaxStep").double_value ();
1212
1213 dae.set_maxstep (maxstep);
1214
1215 //Set initial step
1216 if (! options.getfield ("InitialStep").isempty ())
1217 {
1218 realtype initialstep = options.getfield ("InitialStep").double_value ();
1219
1220 dae.set_initialstep (initialstep);
1221 }
1222
1223 //Set max order FIXME: it doesn't work
1224 int maxorder = options.getfield ("MaxOrder").int_value ();
1225
1226 dae.set_maxorder (maxorder);
1227
1228 //Set Refine
1229 const int refine = options.getfield ("Refine").int_value ();
1230
1231 bool haverefine = (refine > 1);
1232
1233 octave_value output_fcn;
1234 ColumnVector outputsel;
1235
1236 // OutputFcn
1237 bool haveoutputfunction
1238 = options.getfield ("haveoutputfunction").bool_value ();
1239
1240 if (haveoutputfunction)
1241 output_fcn = options.getfield ("OutputFcn");
1242
1243 // OutputSel
1244 bool haveoutputsel = options.getfield ("haveoutputselection").bool_value ();
1245
1246 if (haveoutputsel)
1247 outputsel = options.getfield ("OutputSel").vector_value ();
1248
1249 octave_value event_fcn;
1250
1251 // Events
1252 bool haveeventfunction
1253 = options.getfield ("haveeventfunction").bool_value ();
1254
1255 if (haveeventfunction)
1256 event_fcn = options.getfield ("Events");
1257
1258 // Set up linear solver
1259 dae.set_up (y0);
1260
1261 // Integrate
1262 retval = dae.integrate (numt, tspan, y0, yp0, refine,
1263 haverefine, haveoutputfunction,
1264 output_fcn, haveoutputsel, outputsel,
1265 haveeventfunction, event_fcn, num_event_args);
1266
1267 // Statistics
1268 bool havestats = options.getfield ("havestats").bool_value ();
1269
1270 if (havestats)
1271 dae.print_stat ();
1272
1273 return retval;
1274 }
1275
1276#endif
1277
1278DEFUN_DLD (__ode15__, args, ,
1279 doc: /* -*- texinfo -*-
1280@deftypefn {} {@var{t}, @var{y} =} __ode15__ (@var{fun}, @var{tspan}, @var{y0}, @var{yp0}, @var{options}, @var{num_event_args})
1281Undocumented internal function.
1282@end deftypefn */)
1283{
1284
1285#if defined (HAVE_SUNDIALS)
1286
1287 // Check number of parameters
1288 if (args.length () != 6)
1289 print_usage ();
1290
1291 // Check odefun
1292 octave_value ida_fcn = args(0);
1293
1294 if (! ida_fcn.is_function_handle ())
1295 error ("__ode15__: odefun must be a function handle");
1296
1297 // Check input tspan
1298 ColumnVector tspan
1299 = args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers");
1300
1301 octave_idx_type numt = tspan.numel ();
1302
1303 realtype t0 = tspan (0);
1304
1305 if (numt < 2)
1306 error ("__ode15__: TRANGE must contain at least 2 elements");
1307 else if (tspan.issorted () == UNSORTED || tspan(0) == tspan(numt - 1))
1308 error ("__ode15__: TRANGE must be strictly monotonic");
1309
1310 // input y0 and yp0
1311 ColumnVector y0
1312 = args(2).xvector_value ("__ode15__: initial state y0 must be a vector");
1313
1314 ColumnVector yp0
1315 = args(3).xvector_value ("__ode15__: initial state yp0 must be a vector");
1316
1317
1318 if (y0.numel () != yp0.numel ())
1319 error ("__ode15__: initial state y0 and yp0 must have the same length");
1320 else if (y0.numel () < 1)
1321 error ("__ode15__: initial state yp0 must be a vector or a scalar");
1322
1323
1324 if (! args(4).isstruct ())
1325 error ("__ode15__: OPTS argument must be a structure");
1326
1327 octave_scalar_map options
1328 = args(4).xscalar_map_value ("__ode15__: OPTS argument must be a scalar structure");
1329
1330 // Provided number of arguments in the ode callback function
1331 octave_idx_type num_event_args
1332 = args(5).xidx_type_value ("__ode15__: NUM_EVENT_ARGS must be an integer");
1333
1334 if (num_event_args != 2 && num_event_args != 3)
1335 error ("__ode15__: number of input arguments in event callback must be 2 or 3");
1336
1337 return do_ode15 (ida_fcn, tspan, numt, t0, y0, yp0, options, num_event_args);
1338
1339#else
1340
1341 octave_unused_parameter (args);
1342
1343 err_disabled_feature ("__ode15__", "sundials_ida, sundials_nvecserial");
1344
1345#endif
1346}
1347
1348/*
1349## No test needed for internal helper function.
1350%!assert (1)
1351*/
1352
1353OCTAVE_NAMESPACE_END
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:411
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:534
Array::issorted
OCTARRAY_API sortmode issorted(sortmode mode=UNSORTED) const
Ordering is auto-detected or can be specified.
Definition: Array.cc:2027
Cell
Definition: Cell.h:43
ColumnVector::resize
void resize(octave_idx_type n, const double &rfv=0)
Definition: dColVector.h:112
Matrix
Definition: dMatrix.h:42
Matrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
execution_exception
octave_idx_type
octave_scalar_map
Definition: oct-map.h:162
octave_scalar_map::getfield
octave_value getfield(const std::string &key) const
Definition: oct-map.cc:183
octave_value_list
Definition: ovl.h:44
octave_value_list::resize
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
Definition: ovl.h:117
octave_value_list::length
octave_idx_type length(void) const
Definition: ovl.h:113
octave_value
Definition: ov.h:73
octave_value::vector_value
OCTINTERP_API Array< double > vector_value(bool frc_str_conv=false, bool frc_vec_conv=false) const
octave_value::bool_value
bool bool_value(bool warn=false) const
Definition: ov.h:930
octave_value::int_value
int int_value(bool req_int=false, bool frc_str_conv=false) const
Definition: ov.h:857
octave_value::cell_value
Cell cell_value(void) const
octave_value::is_function_handle
bool is_function_handle(void) const
Definition: ov.h:813
octave_value::isempty
bool isempty(void) const
Definition: ov.h:646
octave_value::double_value
double double_value(bool frc_str_conv=false) const
Definition: ov.h:886
DEFUN_DLD
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Macro to define an at run time dynamically loadable builtin function.
Definition: defun-dld.h:61
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
error
void error(const char *fmt,...)
Definition: error.cc:980
err_disabled_feature
void err_disabled_feature(const std::string &fcn, const std::string &feature, const std::string &pkg)
Definition: errwarn.cc:53
err_user_supplied_eval
void err_user_supplied_eval(const char *name)
Definition: errwarn.cc:152
A
F77_RET_T const F77_INT F77_CMPLX * A
Definition: lo-lapack-proto.h:858
d
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
Definition: lo-slatec-proto.h:39
x
F77_RET_T const F77_DBLE * x
Definition: lo-slatec-proto.h:38
Matrix
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
SparseMatrix
class OCTAVE_API SparseMatrix
Definition: mx-fwd.h:55
ColumnVector
class OCTAVE_API ColumnVector
Definition: mx-fwd.h:45
octave::math::int_multiply_overflow
bool int_multiply_overflow(int a, int b, int *r)
Definition: lo-utils.cc:480
feval
octave_value_list feval(const char *name, const octave_value_list &args, int nargout)
Evaluate an Octave function (built-in or interpreted) and return the list of result values.
Definition: oct-parse.cc:10300
UNSORTED
@ UNSORTED
Definition: oct-sort.h:97
initialize
static void initialize(void)
Definition: octave-config.in.cc:102
ovl
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition: ovl.h:211
octave_stdout
#define octave_stdout
Definition: pager.h:314