23 #if defined (HAVE_CONFIG_H)
38 const double*,
const double&,
43 const double*,
double*,
44 const double&,
double*,
48 const double&,
const double*,
49 const double*,
const double*,
50 const double&,
const double*,
52 double*,
const double&,
60 F77_DBLE&, F77_DBLE*, F77_DBLE*, F77_DBLE&,
61 const F77_INT*,
const F77_DBLE*,
62 const F77_DBLE*, F77_INT&,
63 F77_DBLE*,
const F77_INT&,
64 F77_INT*,
const F77_INT&,
65 const F77_DBLE*,
const F77_INT*,
78 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
86 tmp_deriv.
elem (
i) = deriv[
i];
87 tmp_state.
elem (
i) = state[
i];
90 tmp_delta =
user_fun (tmp_state, tmp_deriv, time, ires);
99 delta[
i] = tmp_delta.
elem (
i);
103 END_INTERRUPT_WITH_EXCEPTIONS;
113 const double *,
const double *,
const double&,
117 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
121 END_INTERRUPT_WITH_EXCEPTIONS;
130 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
139 tmp_deriv.
elem (
i) = deriv[
i];
140 tmp_state.
elem (
i) = state[
i];
147 pd[nn * j +
i] = tmp_pd.
elem (
i, j);
149 END_INTERRUPT_WITH_EXCEPTIONS;
198 if (res.
numel () != x.numel ())
201 (*current_liboctave_error_handler)
202 (
"daspk: inconsistent sizes for state and residual vectors");
211 (*current_liboctave_error_handler)
212 (
"daspk: no user supplied RHS subroutine!");
227 if (eiq == 1 || eiq == 3)
229 if (ccic == 1 || eavfet == 1)
232 lrw = 50 + 9*n + n*n;
247 if (abs_tol_len == 1 && rel_tol_len == 1)
251 else if (abs_tol_len == n && rel_tol_len == n)
259 (*current_liboctave_error_handler)
260 (
"daspk: inconsistent sizes for tolerance arrays");
287 if (maxord > 0 && maxord < 6)
295 (*current_liboctave_error_handler)
296 (
"daspk: invalid value for maximum order");
309 if (ict.
numel () == n)
314 if (val < -2 || val > 2)
317 (*current_liboctave_error_handler)
318 (
"daspk: invalid value for inequality constraint type");
328 (*current_liboctave_error_handler)
329 (
"daspk: inequality constraint types size mismatch");
343 (*current_liboctave_error_handler)
344 (
"daspk: invalid value for enforce inequality constraints option");
357 if (av.
numel () == n)
360 if (eiq == 0 || eiq == 2)
362 else if (eiq == 1 || eiq == 3)
373 (*current_liboctave_error_handler)
374 (
"daspk: algebraic variables size mismatch");
382 (*current_liboctave_error_handler)
383 (
"daspk: invalid value for compute consistent initial condition option");
399 if (av.
numel () == n)
402 if (eiq == 0 || eiq == 2)
404 else if (eiq == 1 || eiq == 3)
418 if (ich.
numel () == 6)
431 (*current_liboctave_error_handler)
432 (
"daspk: invalid initial condition heuristics option");
451 (*current_liboctave_error_handler)
452 (
"daspk: invalid value for print initial condition info option");
531 (*current_liboctave_error_handler)
532 (
"unrecognized value of istate (= %d) returned from ddaspk",
istate);
554 if (n_out > 0 && n > 0)
557 xdot_out.
resize (n_out, n);
599 if (n_out > 0 && n > 0)
602 xdot_out.
resize (n_out, n);
616 double next_crit = tcrit.
elem (0);
618 while (i_out < n_out)
620 bool do_restart =
false;
622 next_out = tout.
elem (i_out);
624 next_crit = tcrit.
elem (i_crit);
629 if (next_crit == next_out)
638 else if (next_crit < next_out)
699 std::ostringstream buf;
706 retval =
"a step was successfully taken in intermediate-output mode.";
710 retval =
"integration completed by stepping exactly to TOUT";
714 retval =
"integration to tout completed by stepping past TOUT";
718 retval =
"initial condition calculation completed successfully";
722 retval =
std::string (
"a large amount of work has been expended (t =")
727 retval =
"the error tolerances are too stringent";
731 retval =
std::string (
"error weight became zero during problem. (t = ")
733 +
"; solution component i vanished, and atol or atol(i) == 0)";
737 retval =
std::string (
"repeated error test failures on the last attempted step (t = ")
742 retval =
std::string (
"the corrector could not converge (t = ")
747 retval =
std::string (
"the matrix of partial derivatives is singular (t = ")
752 retval =
std::string (
"the corrector could not converge (t = ")
753 + t_curr +
"; repeated test failures)";
757 retval =
std::string (
"corrector could not converge because IRES was -1 (t = ")
762 retval =
std::string (
"return requested in user-supplied function (t = ")
767 retval =
"failed to compute consistent initial conditions";
771 retval =
std::string (
"unrecoverable error encountered inside user's PSOL function (t = ")
776 retval =
std::string (
"the Krylov linear system solver failed to converge (t = ")
781 retval =
"unrecoverable error (see printed message)";
785 retval =
"unknown error state";
Array< double > absolute_tolerance(void) const
bool is_empty(void) const
static DAEFunc::DAERHSFunc user_fun
octave_idx_type size(void) const
static octave_idx_type ddaspk_f(const double &time, const double *state, const double *deriv, const double &, double *delta, octave_idx_type &ires, double *, octave_idx_type *)
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)
octave_idx_type numel(void) const
Number of elements in the array.
identity matrix If supplied two scalar respectively For allows like xample val
static octave_idx_type nn
subroutine ddaspk(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
static DAEFunc::DAEJacFunc user_jac
DAERHSFunc function(void) const
octave_idx_type print_initial_condition_info(void) const
Matrix(* DAEJacFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, double cj)
T & elem(octave_idx_type n)
static octave_idx_type ddaspk_psol(const octave_idx_type &, const double &, const double *, const double *, const double *, const double &, const double *, double *, octave_idx_type *, double *, const double &, octave_idx_type &, double *, octave_idx_type *)
octave_idx_type use_initial_condition_heuristics(void) const
#define F77_XFCN(f, F, args)
static octave_idx_type ddaspk_j(const double &time, const double *state, const double *deriv, double *pd, const double &cj, double *, octave_idx_type *)
DAEJacFunc jacobian_function(void) const
double maximum_step_size(void) const
Array< octave_idx_type > iwork
virtual void force_restart(void)
octave_idx_type(* daspk_fcn_ptr)(const double &, const double *, const double *, const double &, double *, octave_idx_type &, double *, octave_idx_type *)
void set_stop_time(double tt)
double initial_step_size(void) const
void resize(const dim_vector &dv, const T &rfv)
octave_idx_type maximum_order(void) const
Array< octave_idx_type > algebraic_variables(void) const
octave_idx_type compute_consistent_initial_condition(void) const
octave_idx_type enforce_inequality_constraints(void) const
octave_idx_type exclude_algebraic_variables_from_error_test(void) const
Array< octave_idx_type > inequality_constraint_types(void) const
void clear_stop_time(void)
static uint32_t state[624]
std::string error_message(void) const
F77_RET_T F77_FUNC(ddaspk, DDASPK)(daspk_fcn_ptr
=val(i)}if ode{val(i)}occurs in table i
octave_idx_type nint_big(double x)
Array< double > relative_tolerance(void) const
octave_idx_type(* daspk_psol_ptr)(const octave_idx_type &, const double &, const double *, const double *, const double *, const double &, const double *, double *, octave_idx_type *, double *, const double &, octave_idx_type &, double *, octave_idx_type *)
Array< octave_idx_type > info
ColumnVector(* DAERHSFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, octave_idx_type &ires)
const T * fortran_vec(void) const
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
octave_idx_type(* daspk_jac_ptr)(const double &, const double *, const double *, double *, const double &, double *, octave_idx_type *)
Array< double > initial_condition_heuristics(void) const