23 #if defined (HAVE_CONFIG_H) 36 const double&,
double*,
F77_INT&,
double*,
40 double*,
const double&,
double*,
F77_INT*);
43 const double*,
const double*,
44 const double*,
const double&,
45 const double*,
double*,
F77_INT*,
46 double*,
const double&,
F77_INT&,
68 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
76 tmp_deriv.
elem (
i) = deriv[
i];
82 tmp_delta =
user_fun (tmp_state, tmp_deriv, time, tmp_ires);
84 ires = octave::to_f77_int (tmp_ires);
93 delta[
i] = tmp_delta.
elem (
i);
97 END_INTERRUPT_WITH_EXCEPTIONS;
107 const double *,
const double *,
const double&,
108 const double *,
double *,
F77_INT *,
double *,
111 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
113 (*current_liboctave_error_handler) (
"daspk: PSOL is not implemented");
115 END_INTERRUPT_WITH_EXCEPTIONS;
122 double *pd,
const double& cj,
double *,
F77_INT *)
124 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
133 tmp_deriv.
elem (
i) = deriv[
i];
141 pd[
nn * j +
i] = tmp_pd.
elem (
i, j);
143 END_INTERRUPT_WITH_EXCEPTIONS;
194 (*current_liboctave_error_handler)
195 (
"daspk: inconsistent sizes for state and residual vectors");
204 (*current_liboctave_error_handler)
205 (
"daspk: no user supplied RHS subroutine!");
217 octave_idx_type eavfet = exclude_algebraic_variables_from_error_test ();
220 if (eiq == 1 || eiq == 3)
222 if (ccic == 1 || eavfet == 1)
225 lrw = 50 + 9*n + n*n;
234 abs_tol = absolute_tolerance ();
235 rel_tol = relative_tolerance ();
240 if (abs_tol_len == 1 && rel_tol_len == 1)
244 else if (abs_tol_len == n && rel_tol_len == n)
251 (*current_liboctave_error_handler)
252 (
"daspk: inconsistent sizes for tolerance arrays");
258 double hmax = maximum_step_size ();
267 double h0 = initial_step_size ();
279 if (maxord > 0 && maxord < 6)
282 iwork(2) = octave::to_f77_int (maxord);
287 (*current_liboctave_error_handler)
288 (
"daspk: invalid value for maximum order");
308 if (val < -2 || val > 2)
311 (*current_liboctave_error_handler)
312 (
"daspk: invalid value for inequality constraint type");
322 (*current_liboctave_error_handler)
323 (
"daspk: inequality constraint types size mismatch");
332 info(9) = octave::to_f77_int (eiq);
337 (*current_liboctave_error_handler)
338 (
"daspk: invalid value for enforce inequality constraints option");
356 if (eiq == 0 || eiq == 2)
358 else if (eiq == 1 || eiq == 3)
362 (
"daspk: invalid value for eiq: %d", eiq);
365 iwork(lid+
i) = (av(
i) ? -1 : 1);
370 (*current_liboctave_error_handler)
371 (
"daspk: algebraic variables size mismatch");
379 (*current_liboctave_error_handler)
380 (
"daspk: invalid value for compute consistent initial condition option");
385 info(10) = octave::to_f77_int (ccic);
401 if (eiq == 0 || eiq == 2)
403 else if (eiq == 1 || eiq == 3)
407 (
"daspk: invalid value for eiq: %d", eiq);
410 iwork(lid+
i) = (av(
i) ? -1 : 1);
414 if (use_initial_condition_heuristics ())
418 if (ich.
numel () == 6)
431 (*current_liboctave_error_handler)
432 (
"daspk: invalid initial condition heuristics option");
446 info(17) = octave::to_f77_int (pici);
451 (*current_liboctave_error_handler)
452 (
"daspk: invalid value for print initial condition info option");
458 DASPK_options::reset =
false;
474 double *dummy =
nullptr;
480 prel_tol, pabs_tol, tmp_istate, prwork,
lrw,
535 (*current_liboctave_error_handler)
536 (
"unrecognized value of istate (= %d) returned from ddaspk",
istate);
558 if (n_out > 0 && n > 0)
561 xdot_out.
resize (n_out, n);
603 if (n_out > 0 && n > 0)
606 xdot_out.
resize (n_out, n);
620 double next_crit = tcrit.
elem (0);
622 while (i_out < n_out)
624 bool do_restart =
false;
626 next_out = tout.
elem (i_out);
628 next_crit = tcrit.
elem (i_crit);
633 if (next_crit == next_out)
642 else if (next_crit < next_out)
703 std::ostringstream buf;
710 retval =
"a step was successfully taken in intermediate-output mode.";
714 retval =
"integration completed by stepping exactly to TOUT";
718 retval =
"integration to tout completed by stepping past TOUT";
722 retval =
"initial condition calculation completed successfully";
726 retval =
"a large amount of work has been expended (t =" + t_curr +
')';
730 retval =
"the error tolerances are too stringent";
734 retval =
"error weight became zero during problem. (t = " + t_curr +
735 "; solution component i vanished, and atol or atol(i) == 0)";
739 retval =
"repeated error test failures on the last attempted step (t = " 744 retval =
"the corrector could not converge (t = " + t_curr +
')';
748 retval =
"the matrix of partial derivatives is singular (t = " + t_curr +
753 retval =
"the corrector could not converge (t = " + t_curr +
754 "; repeated test failures)";
758 retval =
"corrector could not converge because IRES was -1 (t = " 763 retval =
"return requested in user-supplied function (t = " + t_curr +
768 retval =
"failed to compute consistent initial conditions";
772 retval =
"unrecoverable error encountered inside user's PSOL function (t = " 777 retval =
"the Krylov linear system solver failed to converge (t = " 782 retval =
"unrecoverable error (see printed message)";
786 retval =
"unknown error state";
static DAEFunc::DAERHSFunc user_fun
F77_INT(* daspk_jac_ptr)(const double &, const double *, const double *, double *, const double &, double *, F77_INT *)
ColumnVector do_integrate(double t)
Matrix integrate(const ColumnVector &tout, Matrix &xdot_out)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
identity matrix If supplied two scalar respectively For allows like xample val
std::string error_message(void) const
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
static DAEFunc::DAEJacFunc user_jac
const T * fortran_vec(void) const
T & elem(octave_idx_type n)
F77_INT(* daspk_psol_ptr)(const F77_INT &, const double &, const double *, const double *, const double *, const double &, const double *, double *, F77_INT *, double *, const double &, F77_INT &, double *, F77_INT *)
#define F77_XFCN(f, F, args)
subroutine DDASPK(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
F77_RET_T(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_DBLE *, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
static F77_INT ddaspk_f(const double &time, const double *state, const double *deriv, const double &, double *delta, F77_INT &ires, double *, F77_INT *)
octave_value resize(const dim_vector &dv, bool fill=false) const
octave_idx_type size(void) const
Matrix(* DAEJacFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, double cj)
virtual void force_restart(void)
void set_stop_time(double tt)
DAEJacFunc jacobian_function(void) const
F77_INT(* daspk_fcn_ptr)(const double &, const double *, const double *, const double &, double *, F77_INT &, double *, F77_INT *)
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
void clear_stop_time(void)
static uint32_t state[624]
static F77_INT ddaspk_j(const double &time, const double *state, const double *deriv, double *pd, const double &cj, double *, F77_INT *)
Array< octave_f77_int_type > iwork
Array< octave_f77_int_type > info
F77_RET_T F77_FUNC(ddaspk, DDASPK)(daspk_fcn_ptr
octave_f77_int_type F77_INT
octave_idx_type nint_big(double x)
static F77_INT ddaspk_psol(const F77_INT &, const double &, const double *, const double *, const double *, const double &, const double *, double *, F77_INT *, double *, const double &, F77_INT &, double *, F77_INT *)
DAERHSFunc function(void) const
ColumnVector(* DAERHSFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, octave_idx_type &ires)
octave_idx_type numel(void) const
Number of elements in the array.
Vector representing the dimensions (size) of an Array.
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string