Cantera  2.0
IDA_Solver.cpp
Go to the documentation of this file.
1 /**
2  * @file IDA_Solver.cpp
3  *
4  */
5 
6 // Copyright 2006 California Institute of Technology
7 
8 #include "cantera/numerics/IDA_Solver.h"
9 #include "cantera/base/stringUtils.h"
10 
11 #include <iostream>
12 
13 #if SUNDIALS_VERSION >= 24
14 #include "sundials/sundials_types.h"
15 #include "sundials/sundials_math.h"
16 #include "ida/ida.h"
17 #include "ida/ida_dense.h"
18 #include "ida/ida_spgmr.h"
19 #include "ida/ida_band.h"
20 #include "nvector/nvector_serial.h"
21 
22 using namespace std;
23 
24 #if SUNDIALS_VERSION < 25
25 typedef int sd_size_t;
26 #else
27 typedef long int sd_size_t;
28 #endif
29 
30 
31 inline static N_Vector nv(void* x)
32 {
33  return reinterpret_cast<N_Vector>(x);
34 }
35 
36 namespace Cantera
37 {
38 
39 /**
40  * A simple class to hold an array of parameter values and a pointer to
41  * an instance of a subclass of ResidEval.
42  */
43 class ResidData
44 {
45 
46 public:
47 
48  ResidData(ResidJacEval* f, IDA_Solver* s, int npar = 0) {
49  m_func = f;
50  m_solver = s;
51  }
52 
53  virtual ~ResidData() {
54  }
55 
56  ResidJacEval* m_func;
57  IDA_Solver* m_solver;
58 };
59 }
60 
61 //======================================================================================================================
62 extern "C" {
63  //! Function called by IDA to evaluate the residual, given y and ydot.
64  /*!
65  * IDA allows passing in a void* pointer to access external data. Instead of requiring the user to provide a
66  * residual function directly to IDA (which would require using
67  * the sundials data types N_Vector, etc.), we define this function as the single function that IDA always calls. The
68  * real evaluation of the residual is done by an instance of a subclass of ResidEval, passed in to this
69  * function as a pointer in the parameters.
70  *
71  * FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its residual routines ------
72  * A IDAResFn res should return a value of 0 if successful, a positive
73  * value if a recoverable error occured (e.g. yy has an illegal value),
74  * or a negative value if a nonrecoverable error occured. In the latter
75  * case, the program halts. If a recoverable error occured, the integrator
76  * will attempt to correct and retry.
77  */
78  static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data)
79  {
80  double* ydata = NV_DATA_S(y);
81  double* ydotdata = NV_DATA_S(ydot);
82  double* rdata = NV_DATA_S(r);
83  Cantera::ResidData* d = (Cantera::ResidData*) f_data;
84  Cantera::ResidJacEval* f = d->m_func;
85  Cantera::IDA_Solver* s = d->m_solver;
86  double delta_t = s->getCurrentStepFromIDA();
87  // TODO evaluate evalType. Assumed to be Base_ResidEval
88  int retn = 0;
89  int flag = f->evalResidNJ(t, delta_t, ydata, ydotdata, rdata);
90  if (flag < 0) {
91  // This signals to IDA that a nonrecoverable error has occurred.
92  retn = flag;
93  }
94  return retn;
95  }
96 
97  //! Function called by by IDA to evaluate the Jacobian, given y and ydot.
98  /*!
99  *
100  *
101  * typedef int (*IDADlsDenseJacFn)(sd_size_t N, realtype t, realtype c_j,
102  * N_Vector y, N_Vector yp, N_Vector r,
103  * DlsMat Jac, void *user_data,
104  * N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
105  *
106  * A IDADlsDenseJacFn should return
107  * 0 if successful,
108  * a positive int if a recoverable error occurred, or
109  * a negative int if a nonrecoverable error occurred.
110  * In the case of a recoverable error return, the integrator will
111  * attempt to recover by reducing the stepsize (which changes cj).
112  */
113  static int ida_jacobian(sd_size_t nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r,
114  DlsMat Jac, void* f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
115  {
116  doublereal* ydata = NV_DATA_S(y);
117  doublereal* ydotdata = NV_DATA_S(ydot);
118  doublereal* rdata = NV_DATA_S(r);
119  Cantera::ResidData* d = (Cantera::ResidData*) f_data;
120  Cantera::ResidJacEval* f = d->m_func;
121  doublereal* const* colPts = Jac->cols;
122  Cantera::IDA_Solver* s = d->m_solver;
123  double delta_t = s->getCurrentStepFromIDA();
124  // printf(" delta_t = %g 1/cj = %g\n", delta_t, 1.0/c_j);
125  f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata);
126  return 0;
127  }
128 
129 
130 }
131 
132 namespace Cantera
133 {
134 
135 //====================================================================================================================
136 /*
137  * Constructor. Default settings: dense jacobian, no user-supplied
138  * Jacobian function, Newton iteration.
139  */
140 IDA_Solver::IDA_Solver(ResidJacEval& f) :
141  DAE_Solver(f),
142  m_ida_mem(0),
143  m_t0(0.0),
144  m_y(0),
145  m_ydot(0),
146  m_id(0),
147  m_constraints(0),
148  m_abstol(0),
149  m_type(0),
150  m_itol(IDA_SS),
151  m_iter(0),
152  m_reltol(1.e-9),
153  m_abstols(1.e-15),
154  m_nabs(0),
155  m_hmax(0.0),
156  m_hmin(0.0),
157  m_h0(0.0),
158  m_maxsteps(20000),
159  m_maxord(0),
160  m_formJac(0),
161  m_tstop(0.0),
162  m_told_old(0.0),
163  m_told(0.0),
164  m_tcurrent(0.0),
165  m_deltat(0.0),
166  m_maxErrTestFails(-1),
167  m_maxNonlinIters(0),
168  m_maxNonlinConvFails(-1),
169  m_setSuppressAlg(0),
170  m_fdata(0),
171  m_mupper(0),
172  m_mlower(0)
173 {
174 }
175 //====================================================================================================================
176 IDA_Solver::~IDA_Solver()
177 {
178  if (m_ida_mem) {
179  IDAFree(&m_ida_mem);
180  }
181  if (m_y) {
182  N_VDestroy_Serial(nv(m_y));
183  }
184  if (m_ydot) {
185  N_VDestroy_Serial(nv(m_ydot));
186  }
187  if (m_abstol) {
188  N_VDestroy_Serial(nv(m_abstol));
189  }
190  if (m_constraints) {
191  N_VDestroy_Serial(nv(m_constraints));
192  }
193  delete m_fdata;
194 }
195 //====================================================================================================================
196 doublereal IDA_Solver::solution(int k) const
197 {
198  return NV_Ith_S(nv(m_y),k);
199 }
200 //====================================================================================================================
201 const doublereal* IDA_Solver::solutionVector() const
202 {
203  return NV_DATA_S(nv(m_y));
204 }
205 //====================================================================================================================
206 doublereal IDA_Solver::derivative(int k) const
207 {
208  return NV_Ith_S(nv(m_ydot),k);
209 }
210 //====================================================================================================================
211 const doublereal* IDA_Solver::derivativeVector() const
212 {
213  return NV_DATA_S(nv(m_ydot));
214 }
215 //====================================================================================================================
216 
217 void IDA_Solver::setTolerances(double reltol, double* abstol)
218 {
219  m_itol = IDA_SV;
220  if (!m_abstol) {
221  m_abstol = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
222  }
223  for (int i = 0; i < m_neq; i++) {
224  NV_Ith_S(nv(m_abstol), i) = abstol[i];
225  }
226  m_reltol = reltol;
227  if (m_ida_mem) {
228  int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
229  if (flag != IDA_SUCCESS) {
230  throw IDA_Err("Memory allocation failed.");
231  }
232  }
233 }
234 //====================================================================================================================
235 void IDA_Solver::setTolerances(doublereal reltol, doublereal abstol)
236 {
237  m_itol = IDA_SS;
238  m_reltol = reltol;
239  m_abstols = abstol;
240  if (m_ida_mem) {
241  int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
242  if (flag != IDA_SUCCESS) {
243  throw IDA_Err("Memory allocation failed.");
244  }
245  }
246 }
247 //====================================================================================================================
248 void IDA_Solver::setLinearSolverType(int solverType)
249 {
250  m_type = solverType;
251 }
252 //====================================================================================================================
253 void IDA_Solver::setDenseLinearSolver()
254 {
255  setLinearSolverType(0);
256 }
257 //====================================================================================================================
258 void IDA_Solver::setBandedLinearSolver(int m_upper, int m_lower)
259 {
260  m_type = 2;
261  m_upper = m_mupper;
262  m_mlower = m_lower;
263 }
264 //====================================================================================================================
265 void IDA_Solver::setMaxOrder(int n)
266 {
267  m_maxord = n;
268 }
269 //====================================================================================================================
270 void IDA_Solver::setMaxNumSteps(int n)
271 {
272  m_maxsteps = n;
273 }
274 //====================================================================================================================
275 void IDA_Solver::setInitialStepSize(doublereal h0)
276 {
277  m_h0 = h0;
278 }
279 //====================================================================================================================
280 void IDA_Solver::setStopTime(doublereal tstop)
281 {
282  m_tstop = tstop;
283 }
284 //====================================================================================================================
285 doublereal IDA_Solver::getCurrentStepFromIDA()
286 {
287  doublereal hcur;
288  IDAGetCurrentStep(m_ida_mem, &hcur);
289  return hcur;
290 }
291 //====================================================================================================================
292 void IDA_Solver::setJacobianType(int formJac)
293 {
294  m_formJac = formJac;
295  if (m_ida_mem) {
296  if (m_formJac == 1) {
297  int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
298  if (flag != IDA_SUCCESS) {
299  throw IDA_Err("IDADlsSetDenseJacFn failed.");
300  }
301  }
302  }
303 }
304 //====================================================================================================================
305 void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails)
306 {
307  m_maxErrTestFails = maxErrTestFails;
308 }
309 //====================================================================================================================
310 void IDA_Solver::setMaxNonlinIterations(int n)
311 {
312  m_maxNonlinIters = n;
313 }
314 //====================================================================================================================
315 void IDA_Solver::setMaxNonlinConvFailures(int n)
316 {
317  m_maxNonlinConvFails = n;
318 }
319 //====================================================================================================================
320 void IDA_Solver::inclAlgebraicInErrorTest(bool yesno)
321 {
322  if (yesno) {
323  m_setSuppressAlg = 0;
324  } else {
325  m_setSuppressAlg = 1;
326  }
327 }
328 
329 //====================================================================================================================
330 void IDA_Solver::init(doublereal t0)
331 {
332 
333  m_t0 = t0;
334  m_told = t0;
335  m_told_old = t0;
336  m_tcurrent = t0;
337  if (m_y) {
338  N_VDestroy_Serial(nv(m_y));
339  }
340  if (m_ydot) {
341  N_VDestroy_Serial(nv(m_ydot));
342  }
343  if (m_id) {
344  N_VDestroy_Serial(nv(m_id));
345  }
346  if (m_constraints) {
347  N_VDestroy_Serial(nv(m_constraints));
348  }
349 
350  m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
351  m_ydot = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
352  m_constraints = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
353 
354  for (int i=0; i<m_neq; i++) {
355  NV_Ith_S(nv(m_y), i) = 0.0;
356  NV_Ith_S(nv(m_ydot), i) = 0.0;
357  NV_Ith_S(nv(m_constraints), i) = 0.0;
358  }
359 
360  // get the initial conditions
361  m_resid.getInitialConditions(m_t0, NV_DATA_S(nv(m_y)), NV_DATA_S(nv(m_ydot)));
362 
363  if (m_ida_mem) {
364  IDAFree(&m_ida_mem);
365  }
366 
367  /* Call IDACreate */
368  m_ida_mem = IDACreate();
369 
370  int flag = 0;
371 
372 
373 
374  if (m_itol == IDA_SV) {
375 #if SUNDIALS_VERSION <= 23
376  // vector atol
377  flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
378  m_itol, m_reltol, nv(m_abstol));
379  if (flag != IDA_SUCCESS) {
380  if (flag == IDA_MEM_FAIL) {
381  throw IDA_Err("Memory allocation failed.");
382  } else if (flag == IDA_ILL_INPUT) {
383  throw IDA_Err("Illegal value for IDAMalloc input argument.");
384  } else {
385  throw IDA_Err("IDAMalloc failed.");
386  }
387  }
388 
389 #elif SUNDIALS_VERSION >= 24
390  flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
391  if (flag != IDA_SUCCESS) {
392  if (flag == IDA_MEM_FAIL) {
393  throw IDA_Err("Memory allocation failed.");
394  } else if (flag == IDA_ILL_INPUT) {
395  throw IDA_Err("Illegal value for IDAMalloc input argument.");
396  } else {
397  throw IDA_Err("IDAMalloc failed.");
398  }
399  }
400  flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
401  if (flag != IDA_SUCCESS) {
402  throw IDA_Err("Memory allocation failed.");
403  }
404 #endif
405  } else {
406 #if SUNDIALS_VERSION <= 23
407  // scalar atol
408  flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
409  m_itol, m_reltol, &m_abstols);
410  if (flag != IDA_SUCCESS) {
411  if (flag == IDA_MEM_FAIL) {
412  throw IDA_Err("Memory allocation failed.");
413  } else if (flag == IDA_ILL_INPUT) {
414  throw IDA_Err("Illegal value for IDAMalloc input argument.");
415  } else {
416  throw IDA_Err("IDAMalloc failed.");
417  }
418  }
419 
420 #elif SUNDIALS_VERSION >= 24
421  flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
422  if (flag != IDA_SUCCESS) {
423  if (flag == IDA_MEM_FAIL) {
424  throw IDA_Err("Memory allocation failed.");
425  } else if (flag == IDA_ILL_INPUT) {
426  throw IDA_Err("Illegal value for IDAMalloc input argument.");
427  } else {
428  throw IDA_Err("IDAMalloc failed.");
429  }
430  }
431  flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
432  if (flag != IDA_SUCCESS) {
433  throw IDA_Err("Memory allocation failed.");
434  }
435 #endif
436  }
437 
438  //-----------------------------------
439  // set the linear solver type
440  //-----------------------------------
441 
442  if (m_type == 1 || m_type == 0) {
443  long int N = m_neq;
444  flag = IDADense(m_ida_mem, N);
445  if (flag) {
446  throw IDA_Err("IDADense failed");
447  }
448  } else if (m_type == 2) {
449  long int N = m_neq;
450  long int nu = m_mupper;
451  long int nl = m_mlower;
452  IDABand(m_ida_mem, N, nu, nl);
453  } else {
454  throw IDA_Err("unsupported linear solver type");
455  }
456 
457  if (m_formJac == 1) {
458  flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
459  if (flag != IDA_SUCCESS) {
460  throw IDA_Err("IDADlsSetDenseJacFn failed.");
461  }
462  }
463 
464  // pass a pointer to func in m_data
465  m_fdata = new ResidData(&m_resid, this, m_resid.nparams());
466 #if SUNDIALS_VERSION <= 23
467  flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
468  if (flag != IDA_SUCCESS) {
469  throw IDA_Err("IDASetRdata failed.");
470  }
471 #elif SUNDIALS_VERSION >= 24
472  flag = IDASetUserData(m_ida_mem, (void*)m_fdata);
473  if (flag != IDA_SUCCESS) {
474  throw IDA_Err("IDASetUserData failed.");
475  }
476 #endif
477 
478  // set options
479  if (m_maxord > 0) {
480  flag = IDASetMaxOrd(m_ida_mem, m_maxord);
481  if (flag != IDA_SUCCESS) {
482  throw IDA_Err("IDASetMaxOrd failed.");
483  }
484  }
485  if (m_maxsteps > 0) {
486  flag = IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
487  if (flag != IDA_SUCCESS) {
488  throw IDA_Err("IDASetMaxNumSteps failed.");
489  }
490  }
491  if (m_h0 > 0.0) {
492  flag = IDASetInitStep(m_ida_mem, m_h0);
493  if (flag != IDA_SUCCESS) {
494  throw IDA_Err("IDASetInitStep failed.");
495  }
496  }
497  if (m_tstop > 0.0) {
498  flag = IDASetStopTime(m_ida_mem, m_tstop);
499  if (flag != IDA_SUCCESS) {
500  throw IDA_Err("IDASetStopTime failed.");
501  }
502  }
503  if (m_maxErrTestFails >= 0) {
504  flag = IDASetMaxErrTestFails(m_ida_mem, m_maxErrTestFails);
505  if (flag != IDA_SUCCESS) {
506  throw IDA_Err("IDASetMaxErrTestFails failed.");
507  }
508  }
509  if (m_maxNonlinIters > 0) {
510  flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
511  if (flag != IDA_SUCCESS) {
512  throw IDA_Err("IDASetmaxNonlinIters failed.");
513  }
514  }
515  if (m_maxNonlinConvFails >= 0) {
516  flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
517  if (flag != IDA_SUCCESS) {
518  throw IDA_Err("IDASetMaxConvFails failed.");
519  }
520  }
521  if (m_setSuppressAlg != 0) {
522  flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
523  if (flag != IDA_SUCCESS) {
524  throw IDA_Err("IDASetSuppressAlg failed.");
525  }
526  }
527 
528 
529 
530 }
531 //====================================================================================================================
532 // Calculate consistent value of the starting solution given the starting solution derivatives
533 /*
534  * This method may be called if the initial conditions do not
535  * satisfy the residual equation F = 0. Given the derivatives
536  * of all variables, this method computes the initial y
537  * values.
538  */
539 void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout)
540 {
541  int icopt = IDA_Y_INIT;
542  doublereal tout1 = tout;
543  if (tout == 0.0) {
544  double h0 = 1.0E-5;
545  if (m_h0 > 0.0) {
546  h0 = m_h0;
547  }
548  tout1 = m_t0 + h0;
549  }
550 
551  int flag = IDACalcIC(m_ida_mem, icopt, tout1);
552  if (flag != IDA_SUCCESS) {
553  throw IDA_Err("IDACalcIC failed: error = " + int2str(flag));
554  }
555 
556 
557  flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot));
558  if (flag != IDA_SUCCESS) {
559  throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
560  }
561  doublereal* yy = NV_DATA_S(nv(m_y));
562  doublereal* yyp = NV_DATA_S(nv(m_ydot));
563 
564  for (int i = 0; i < m_neq; i++) {
565  y[i] = yy[i];
566  yp[i] = yyp[i];
567  }
568 }
569 //====================================================================================================================
570 /*
571  * This method may be called if the initial conditions do not
572  * satisfy the residual equation F = 0. Given the initial
573  * values of all differential variables, it computes the
574  * initial values of all algebraic variables and the initial
575  * derivatives of all differential variables.
576  *
577  * @param y Calculated value of the solution vector after the procedure ends
578  * @param yp Calculated value of the solution derivative after the procedure
579  * @param The first value of t at which a soluton will be
580  * requested (from IDASolve). (This is needed here to
581  * determine the direction of integration and rough scale
582  * in the independent variable t.
583  */
584 void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout)
585 {
586 
587  int icopt = IDA_YA_YDP_INIT;
588  doublereal tout1 = tout;
589  if (tout == 0.0) {
590  double h0 = 1.0E-5;
591  if (m_h0 > 0.0) {
592  h0 = m_h0;
593  }
594  tout1 = m_t0 + h0;
595  }
596 
597  int flag = IDACalcIC(m_ida_mem, icopt, tout1);
598  if (flag != IDA_SUCCESS) {
599  throw IDA_Err("IDACalcIC failed: error = " + int2str(flag));
600  }
601 
602 
603  flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot));
604  if (flag != IDA_SUCCESS) {
605  throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
606  }
607  doublereal* yy = NV_DATA_S(nv(m_y));
608  doublereal* yyp = NV_DATA_S(nv(m_ydot));
609 
610  for (int i = 0; i < m_neq; i++) {
611  y[i] = yy[i];
612  yp[i] = yyp[i];
613  }
614 }
615 //====================================================================================================================
616 int IDA_Solver::solve(double tout)
617 {
618  double tretn;
619  int flag;
620  flag = IDASetStopTime(m_ida_mem, tout);
621  if (flag != IDA_SUCCESS) {
622  throw IDA_Err(" IDA error encountered.");
623  }
624  do {
625  if (tout <= m_tcurrent) {
626  throw IDA_Err(" tout <= tcurrent");
627  }
628  m_told_old = m_told;
629  m_told = m_tcurrent;
630  flag = IDASolve(m_ida_mem, tout, &tretn, nv(m_y), nv(m_ydot), IDA_ONE_STEP);
631  if (flag < 0) {
632  throw IDA_Err(" IDA error encountered.");
633  } else if (flag == IDA_TSTOP_RETURN) {
634  // we've reached our goal, and have actually integrated past it
635  } else if (flag == IDA_ROOT_RETURN) {
636  // not sure what to do with this yet
637  } else if (flag == IDA_WARNING) {
638  throw IDA_Err(" IDA Warning encountered.");
639  }
640  m_tcurrent = tretn;
641  m_deltat = m_tcurrent - m_told;
642  } while (tretn < tout);
643 
644  if (flag != IDA_SUCCESS && flag != IDA_TSTOP_RETURN) {
645  throw IDA_Err(" IDA error encountered.");
646  }
647  return flag;
648 }
649 //====================================================================================================================
650 double IDA_Solver::step(double tout)
651 {
652  double t;
653  int flag;
654  if (tout <= m_tcurrent) {
655  throw IDA_Err(" tout <= tcurrent");
656  }
657  m_told_old = m_told;
658  m_told = m_tcurrent;
659  flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP);
660  if (flag < 0) {
661  throw IDA_Err(" IDA error encountered.");
662  } else if (flag == IDA_TSTOP_RETURN) {
663  // we've reached our goal, and have actually integrated past it
664  } else if (flag == IDA_ROOT_RETURN) {
665  // not sure what to do with this yet
666  } else if (flag == IDA_WARNING) {
667  throw IDA_Err(" IDA Warning encountered.");
668  }
669  m_tcurrent = t;
670  m_deltat = m_tcurrent - m_told;
671  return t;
672 }
673 //====================================================================================================================
674 doublereal IDA_Solver::getOutputParameter(int flag) const
675 {
676  long int lenrw, leniw;
677  switch (flag) {
678  case REAL_WORKSPACE_SIZE:
679  flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
680  return doublereal(lenrw);
681  break;
682  }
683  return 0.0;
684 }
685 //====================================================================================================================
686 
687 }
688 #endif
Generated by   doxygen 1.8.2