1 #include "MUQ/Modeling/ODEBase.h"
2
3 #include <cvodes/cvodes.h> // prototypes for CVODE fcts. and consts.
4 #include <cvodes/cvodes_spgmr.h> // prototypes & constants for CVSPGMR solver
5 #include <cvodes/cvodes_spbcgs.h> // prototypes & constants for CVSPBCG solver
6 #include <cvodes/cvodes_sptfqmr.h> // prototypes & constants for SPTFQMR solver
7 #include <cvodes/cvodes_dense.h> // prototype for CVDense
8
9 #include <sundials/sundials_types.h> // definition of type
10 #include <sundials/sundials_math.h> // contains the macros ABS, SQR, and EXP
11
12 #if MUQ_HAS_PARCER==1
13 #include <nvector/nvector_parallel.h>
14 #endif
15
16 namespace pt = boost::property_tree;
17 using namespace muq::Modeling;
18
19 #if MUQ_HAS_PARCER==1
ODEBase(std::shared_ptr<ModPiece> const & rhs,Eigen::VectorXi const & inputSizes,Eigen::VectorXi const & outputSizes,pt::ptree const & pt,std::shared_ptr<parcer::Communicator> const & comm)20 ODEBase::ODEBase(std::shared_ptr<ModPiece> const& rhs, Eigen::VectorXi const& inputSizes, Eigen::VectorXi const& outputSizes, pt::ptree const& pt, std::shared_ptr<parcer::Communicator> const& comm) :
21 #else
22 ODEBase::ODEBase(std::shared_ptr<ModPiece> const& rhs, Eigen::VectorXi const& inputSizes, Eigen::VectorXi const& outputSizes, pt::ptree const& pt) :
23 #endif
24 ModPiece(inputSizes, outputSizes),
25 rhs(rhs),
26 reltol(pt.get<double>("RelativeTolerance", 1.0e-8)),
27 abstol(pt.get<double>("AbsoluteTolerance", 1.0e-8)), maxStepSize(pt.get<double>("MaxStepSize", 1.0)),
28 maxNumSteps(pt.get<unsigned int>("MaxNumSteps", 500)),
29 autonomous(pt.get<bool>("Autonomous", true)),
30 checkPtGap(pt.get<unsigned int>("CheckPointGap", 50))
31 #if MUQ_HAS_PARCER==1
32 , globalSize(pt.get<unsigned int>("GlobalSize", std::numeric_limits<unsigned int>::quiet_NaN())),
33 comm(comm)
34 #endif
35 {
36 // we must know the number of inputs for the rhs and it must have at least one (the state)
37 assert(rhs->numInputs>0);
38
39 // determine the multistep method and the nonlinear solver
40 const std::string& multiStepMethod = pt.get<std::string>("MultistepMethod", "BDF");
41 assert(multiStepMethod.compare("Adams")==0 || multiStepMethod.compare("BDF")==0); // Adams or BDF
42 multiStep = (multiStepMethod.compare("BDF")==0) ? CV_BDF : CV_ADAMS;
43 const std::string& nonlinearSolver = pt.get<std::string>("NonlinearSolver", "Newton");
44 assert(nonlinearSolver.compare("Iter")==0 || nonlinearSolver.compare("Newton")==0); // Iter or Newton
45 solveMethod = (nonlinearSolver.compare("Newton")==0) ? CV_NEWTON : CV_FUNCTIONAL;
46
47 // determine the method of thelinear solver
48 const std::string& linearSolver = pt.get<std::string>("LinearSolver", "Dense");
49 if( linearSolver.compare("Dense")==0 ) {
50 slvr = LinearSolver::Dense;
51 } else if( linearSolver.compare("SPGMR")==0 ) {
52 slvr = LinearSolver::SPGMR;
53 } else if( linearSolver.compare("SPBCG")==0 ) {
54 slvr = LinearSolver::SPBCG;
55 } else if( linearSolver.compare("SPTFQMR")==0 ) {
56 slvr = LinearSolver::SPTFQMR;
57 } else {
58 std::cerr << "\nInvalid CVODES linear solver type. Options are Dense, SPGMR, SPBCG, or SPTFQMR\n\n";
59 assert(false);
60 }
61 }
62
~ODEBase()63 ODEBase::~ODEBase() {}
64
CheckFlag(void * flagvalue,std::string const & funcname,unsigned int const opt) const65 bool ODEBase::CheckFlag(void* flagvalue, std::string const& funcname, unsigned int const opt) const {
66 // there are only two options
67 assert(opt==0 || opt==1);
68
69 // check if Sundials function returned nullptr pointer - no memory allocated
70 if( opt==0 && flagvalue==nullptr ) {
71 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", funcname.c_str());
72
73 // return failure
74 return false;
75 }
76
77 // check if flag<0
78 if( opt==1 ) {
79 // get int value
80 int *errflag = (int *) flagvalue;
81
82 if( *errflag<0 ) { // negative indicates Sundials error
83 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", funcname.c_str(), *errflag);
84
85 // return failure
86 return false;
87 }
88 }
89
90 // return success
91 return true;
92 }
93
ErrorHandler(int error_code,const char * module,const char * function,char * msg,void * user_data)94 void ODEBase::ErrorHandler(int error_code, const char *module, const char *function, char *msg, void *user_data) {}
95
CreateSolverMemoryB(void * cvode_mem,double const timeFinal,N_Vector const & lambda,N_Vector const & nvGrad,std::shared_ptr<ODEData> data) const96 int ODEBase::CreateSolverMemoryB(void* cvode_mem, double const timeFinal, N_Vector const& lambda, N_Vector const& nvGrad, std::shared_ptr<ODEData> data) const {
97 int indexB, flag;
98
99 // initialize the adjoint solver
100 flag = CVodeAdjInit(cvode_mem, checkPtGap, CV_HERMITE);
101 assert(CheckFlag(&flag, "CVodeAdjInit", 1));
102
103 // creat adjoint solver
104 flag = CVodeCreateB(cvode_mem, multiStep, solveMethod, &indexB);
105 assert(CheckFlag(&flag, "CVodeCreateB", 1));
106
107 // set the pointer to user-defined data
108 assert(data);
109 flag = CVodeSetUserDataB(cvode_mem, indexB, data.get());
110 assert(CheckFlag(&flag, "CVodeSetUserDataB", 1));
111
112 // set the adjoint right hand side
113 flag = CVodeInitB(cvode_mem, indexB, AdjointRHS, timeFinal, lambda);
114 assert(CheckFlag(&flag, "CVodeInitB", 1));
115
116 // specify the relative and absolute tolerances
117 flag = CVodeSStolerancesB(cvode_mem, indexB, reltol, abstol);
118 assert(CheckFlag(&flag, "CVodeSStolerancesB", 1));
119
120 // set the maximum time step size
121 flag = CVodeSetMaxStepB(cvode_mem, indexB, maxStepSize);
122 assert(CheckFlag(&flag, "CVodeSetMaxStepB", 1));
123
124 flag = CVodeSetMaxNumStepsB(cvode_mem, indexB, 10000000);
125 assert(CheckFlag(&flag, "CVodesSetMaxNumStepsB", 1));
126
127 // determine which linear solver to use
128 if( slvr==LinearSolver::Dense ) { // dense linear solver
129 // specify the dense linear solver
130 #if MUQ_HAS_PARCER==1
131 flag = CVDenseB(cvode_mem, indexB, comm? NV_GLOBLENGTH_P(lambda) : NV_LENGTH_S(lambda));
132 #else
133 flag = CVDenseB(cvode_mem, indexB, NV_LENGTH_S(lambda));
134 #endif
135 assert(CheckFlag(&flag, "CVDenseB", 1));
136
137 // set the Jacobian routine to Jac (user-supplied)
138 flag = CVDlsSetDenseJacFnB(cvode_mem, indexB, AdjointJacobian);
139 assert(CheckFlag(&flag, "CVDlsSetDenseJacFnB", 1));
140 } else { // sparse linear solver
141 if( slvr==LinearSolver::SPGMR ) {
142 flag = CVSpgmrB(cvode_mem, indexB, 0, 0);
143 assert(CheckFlag(&flag, "CVSpgmrB", 1));
144 } else if( slvr==LinearSolver::SPBCG ) {
145 flag = CVSpbcgB(cvode_mem, indexB, 0, 0);
146 assert(CheckFlag(&flag, "CVSpbcgB", 1));
147 } else if( slvr==LinearSolver::SPTFQMR ) {
148 flag = CVSptfqmrB(cvode_mem, indexB, 0, 0);
149 assert(CheckFlag(&flag, "CVSptfqmrB", 1));
150 } else {
151 std::cerr << "\nInvalid CVODES linear solver type. Options are Dense, SPGMR, SPBCG, or SPTFQMR\n\n";
152 assert(false);
153 }
154
155 // set the Jacobian-times-vector function
156 flag = CVSpilsSetJacTimesVecFnB(cvode_mem, indexB, AdjointJacobianAction);
157 assert(CheckFlag(&flag, "CVSpilsSetJacTimesVecFnB", 1));
158 }
159
160 flag = CVodeQuadInitB(cvode_mem, indexB, AdjointQuad, nvGrad);
161 assert(CheckFlag(&flag, "CVodeQuadInitB", 1));
162
163 flag = CVodeSetQuadErrCon(cvode_mem, true);
164 assert(CheckFlag(&flag, "CVodeSetQuadErrCon", 1));
165
166 flag = CVodeQuadSStolerancesB(cvode_mem, indexB, reltol, abstol);
167 assert(CheckFlag(&flag, "CVodeQuadSStolerancesB", 1));
168
169 return indexB;
170 }
171
CreateSolverMemory(void * cvode_mem,N_Vector const & state,std::shared_ptr<ODEData> data) const172 void ODEBase::CreateSolverMemory(void* cvode_mem, N_Vector const& state, std::shared_ptr<ODEData> data) const {
173 // a flag used for error checking
174 int flag;
175
176 // set the pointer to user-defined data
177 assert(data);
178 flag = CVodeSetUserData(cvode_mem, data.get());
179 assert(CheckFlag(&flag, "CVodeSetUserData", 1));
180
181 // tell the solver how to deal with errors
182 flag = CVodeSetErrHandlerFn(cvode_mem, ErrorHandler, data.get());
183 assert(CheckFlag(&flag, "CVodeSetErrHandlerFun", 1));
184
185 // tell the solver how to evaluate the rhs
186 flag = CVodeInit(cvode_mem, EvaluateRHS, 0.0, state);
187 assert(CheckFlag(&flag, "CVodeInit", 1));
188
189 // specify the relative and absolute tolerances
190 flag = CVodeSStolerances(cvode_mem, reltol, abstol);
191 assert(CheckFlag(&flag, "CVodeSStolerances", 1));
192
193 // set the maximum time step size
194 flag = CVodeSetMaxStep(cvode_mem, maxStepSize);
195 assert(CheckFlag(&flag, "CVodeSetMaxStep", 1));
196
197 flag = CVodeSetMaxNumSteps(cvode_mem, maxNumSteps);
198 assert(CheckFlag(&flag, "CVodesSetMaxNumSteps", 1));
199
200 // determine which linear solver to use
201 if( slvr==LinearSolver::Dense ) { // dense linear solver
202 // specify the dense linear solver
203 #if MUQ_HAS_PARCER==1
204 flag = CVDense(cvode_mem, comm? NV_GLOBLENGTH_P(state) : NV_LENGTH_S(state));
205 #else
206 flag = CVDense(cvode_mem, NV_LENGTH_S(state));
207 #endif
208 assert(CheckFlag(&flag, "CVDense", 1));
209
210 // set the Jacobian routine to Jac (user-supplied)
211 flag = CVDlsSetDenseJacFn(cvode_mem, RHSJacobian);
212 assert(CheckFlag(&flag, "CVDlsSetDenseJacFn", 1));
213 } else { // sparse linear solver
214 if( slvr==LinearSolver::SPGMR ) {
215 flag = CVSpgmr(cvode_mem, 0, 0);
216 assert(CheckFlag(&flag, "CVSpgmr", 1));
217 } else if( slvr==LinearSolver::SPBCG ) {
218 flag = CVSpbcg(cvode_mem, 0, 0);
219 assert(CheckFlag(&flag, "CVSpbcg", 1));
220 } else if( slvr==LinearSolver::SPTFQMR ) {
221 flag = CVSptfqmr(cvode_mem, 0, 0);
222 assert(CheckFlag(&flag, "CVSptfqmr", 1));
223 } else {
224 std::cerr << "\nInvalid CVODES linear solver type. Options are Dense, SPGMR, SPBCG, or SPTFQMR\n\n";
225 assert(false);
226 }
227
228 // set the Jacobian-times-vector function
229 flag = CVSpilsSetJacTimesVecFn(cvode_mem, RHSJacobianAction);
230 assert(CheckFlag(&flag, "CVSpilsSetJacTimesVecFn", 1));
231 }
232 }
233
234 #define NVectorMap(eigenmap, comm, vec) Eigen::Map<Eigen::VectorXd> eigenmap( \
235 comm? NV_DATA_P(vec) : NV_DATA_S(vec), \
236 comm? NV_LOCLENGTH_P(vec) : NV_LENGTH_S(vec));
237
EvaluateRHS(realtype time,N_Vector state,N_Vector statedot,void * user_data)238 int ODEBase::EvaluateRHS(realtype time, N_Vector state, N_Vector statedot, void *user_data) {
239 // get the data type
240 ODEData* data = (ODEData*)user_data;
241 assert(data);
242 assert(data->rhs);
243
244 // set the state input
245 #if MUQ_HAS_PARCER==1
246 NVectorMap(stateref, data->comm, state);
247 #else
248 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
249 #endif
250 data->UpdateInputs(stateref, time);
251
252 // evaluate the rhs
253 #if MUQ_HAS_PARCER==1
254 NVectorMap(statedotref, data->comm, statedot);
255 #else
256 Eigen::Map<Eigen::VectorXd> statedotref(NV_DATA_S(statedot), NV_LENGTH_S(statedot));
257 #endif
258 statedotref = data->rhs->Evaluate(data->inputs) [0];
259
260 return 0;
261 }
262
AdjointRHS(realtype time,N_Vector state,N_Vector lambda,N_Vector deriv,void * user_data)263 int ODEBase::AdjointRHS(realtype time, N_Vector state, N_Vector lambda, N_Vector deriv, void *user_data) {
264 // get the data type
265 ODEData* data = (ODEData*)user_data;
266 assert(data);
267 assert(data->rhs);
268
269 // set the state input
270 #if MUQ_HAS_PARCER==1
271 NVectorMap(stateref, data->comm, state);
272 #else
273 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
274 #endif
275 data->UpdateInputs(stateref, time);
276
277 #if MUQ_HAS_PARCER==1
278 NVectorMap(lambdaMap, data->comm, lambda);
279 #else
280 Eigen::Map<Eigen::VectorXd> lambdaMap(NV_DATA_S(lambda), NV_LENGTH_S(lambda));
281 #endif
282 const Eigen::VectorXd& lam = lambdaMap;
283
284 #if MUQ_HAS_PARCER==1
285 NVectorMap(derivMap, data->comm, deriv);
286 #else
287 Eigen::Map<Eigen::VectorXd> derivMap(NV_DATA_S(deriv), NV_LENGTH_S(deriv));
288 #endif
289
290 #if MUQ_HAS_PARCER==1
291 int globalind = 0;
292 int localSize = 0;
293 if( data->comm ) {
294 localSize = NV_LOCLENGTH_P(state);
295 for( unsigned int i=0; i<data->comm->GetSize(); ++i ) {
296 int temp = localSize;
297 data->comm->Bcast(temp, i);
298 globalind += i<data->comm->GetRank()? temp : 0;
299 }
300 } else {
301 localSize = NV_LENGTH_S(state);
302 }
303 #else
304 int globalind = 0;
305 int localSize = NV_LENGTH_S(state);
306 #endif
307 derivMap = -1.0*data->rhs->Gradient(0, data->autonomous? 0 : 1, data->inputs, lam).segment(globalind, localSize);
308
309 return 0;
310 }
311
RHSJacobianAction(N_Vector v,N_Vector Jv,realtype time,N_Vector state,N_Vector rhs,void * user_data,N_Vector tmp)312 int ODEBase::RHSJacobianAction(N_Vector v, N_Vector Jv, realtype time, N_Vector state, N_Vector rhs, void *user_data, N_Vector tmp) {
313 // get the data type
314 ODEData* data = (ODEData*)user_data;
315 assert(data);
316 assert(data->rhs);
317
318 // set the state input
319 #if MUQ_HAS_PARCER==1
320 NVectorMap(stateref, data->comm, state);
321 #else
322 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
323 #endif
324 data->UpdateInputs(stateref, time);
325
326 #if MUQ_HAS_PARCER==1
327 NVectorMap(vref, data->comm, v);
328 #else
329 Eigen::Map<Eigen::VectorXd> vref(NV_DATA_S(v), NV_LENGTH_S(v));
330 #endif
331 const Eigen::VectorXd& veig = vref;
332
333 // compute the jacobain wrt the state
334 #if MUQ_HAS_PARCER==1
335 NVectorMap(Jvref, data->comm, Jv);
336 #else
337 Eigen::Map<Eigen::VectorXd> Jvref(NV_DATA_S(Jv), NV_LENGTH_S(Jv));
338 #endif
339 Jvref = data->rhs->ApplyJacobian(0, data->autonomous? 0 : 1, data->inputs, veig);
340
341 return 0;
342 }
343
AdjointJacobianAction(N_Vector target,N_Vector output,realtype time,N_Vector state,N_Vector lambda,N_Vector adjRhs,void * user_data,N_Vector tmp)344 int ODEBase::AdjointJacobianAction(N_Vector target, N_Vector output, realtype time, N_Vector state, N_Vector lambda, N_Vector adjRhs, void *user_data, N_Vector tmp) {
345 // get the data type
346 ODEData* data = (ODEData*)user_data;
347 assert(data);
348 assert(data->rhs);
349
350 // set the state input
351 #if MUQ_HAS_PARCER==1
352 NVectorMap(stateref, data->comm, state);
353 #else
354 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
355 #endif
356 data->UpdateInputs(stateref, time);
357
358 #if MUQ_HAS_PARCER==1
359 NVectorMap(targetMap, data->comm, target);
360 #else
361 Eigen::Map<Eigen::VectorXd> targetMap(NV_DATA_S(target), NV_LENGTH_S(target));
362 #endif
363 const Eigen::VectorXd& targ = targetMap;
364 #if MUQ_HAS_PARCER==1
365 NVectorMap(outputMap, data->comm, output);
366 #else
367 Eigen::Map<Eigen::VectorXd> outputMap(NV_DATA_S(output), NV_LENGTH_S(output));
368 #endif
369
370 outputMap = -1.0*data->rhs->Gradient(0, data->autonomous? 0 : 1, data->inputs, targ);
371
372 return 0;
373 }
374
RHSJacobian(long int N,realtype time,N_Vector state,N_Vector rhs,DlsMat jac,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)375 int ODEBase::RHSJacobian(long int N, realtype time, N_Vector state, N_Vector rhs, DlsMat jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
376 // get the data type
377 ODEData* data = (ODEData*)user_data;
378 assert(data);
379 assert(data->rhs);
380
381 // set the state input
382 #if MUQ_HAS_PARCER==1
383 NVectorMap(stateref, data->comm, state);
384 #else
385 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
386 #endif
387 data->UpdateInputs(stateref, time);
388
389 // evaluate the jacobian
390 Eigen::Map<Eigen::MatrixXd> jacref(jac->data, jac->M, jac->N);
391 jacref = data->rhs->Jacobian(0, data->autonomous? 0 : 1, data->inputs);
392
393 return 0;
394 }
395
AdjointJacobian(long int N,realtype time,N_Vector state,N_Vector lambda,N_Vector rhs,DlsMat jac,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)396 int ODEBase::AdjointJacobian(long int N, realtype time, N_Vector state, N_Vector lambda, N_Vector rhs, DlsMat jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
397 // get the data type
398 ODEData* data = (ODEData*)user_data;
399 assert(data);
400 assert(data->rhs);
401
402 // set the state input
403 #if MUQ_HAS_PARCER==1
404 NVectorMap(stateref, data->comm, state);
405 #else
406 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
407 #endif
408 data->UpdateInputs(stateref, time);
409
410 // evaluate the jacobian
411 Eigen::Map<Eigen::MatrixXd> jacref(jac->data, jac->M, jac->N);
412 jacref = -1.0*data->rhs->Jacobian(0, data->autonomous? 0 : 1, data->inputs).transpose();
413
414 return 0;
415 }
416
AdjointQuad(realtype time,N_Vector state,N_Vector lambda,N_Vector quadRhs,void * user_data)417 int ODEBase::AdjointQuad(realtype time, N_Vector state, N_Vector lambda, N_Vector quadRhs, void *user_data) {
418 // get the data type
419 ODEData* data = (ODEData*)user_data;
420 assert(data);
421 assert(data->rhs);
422
423 // set the state input
424 #if MUQ_HAS_PARCER==1
425 NVectorMap(stateref, data->comm, state);
426 #else
427 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(state), NV_LENGTH_S(state));
428 #endif
429 data->UpdateInputs(stateref, time);
430
431 #if MUQ_HAS_PARCER==1
432 NVectorMap(lambdaMap, data->comm, lambda);
433 #else
434 Eigen::Map<Eigen::VectorXd> lambdaMap(NV_DATA_S(lambda), NV_LENGTH_S(lambda));
435 #endif
436 const Eigen::VectorXd& lam = lambdaMap;
437
438 #if MUQ_HAS_PARCER==1
439 NVectorMap(nvQuadMap, data->comm, quadRhs);
440 #else
441 Eigen::Map<Eigen::VectorXd> nvQuadMap(NV_DATA_S(quadRhs), NV_LENGTH_S(quadRhs));
442 #endif
443
444 nvQuadMap = data->rhs->Gradient(0, data->autonomous? data->wrtIn : data->wrtIn+1, data->inputs, lam);
445
446 return 0;
447 }
448
SetUpSensitivity(void * cvode_mem,unsigned int const paramSize,N_Vector * sensState) const449 void ODEBase::SetUpSensitivity(void *cvode_mem, unsigned int const paramSize, N_Vector *sensState) const {
450 // initialze the forward sensitivity solver
451 int flag = CVodeSensInit(cvode_mem, paramSize, CV_SIMULTANEOUS, ForwardSensitivityRHS, sensState);
452 assert(CheckFlag(&flag, "CVodeSensInit", 1));
453
454 // set sensitivity tolerances
455 Eigen::VectorXd absTolVec = Eigen::VectorXd::Constant(paramSize, abstol);
456 flag = CVodeSensSStolerances(cvode_mem, reltol, absTolVec.data());
457 assert(CheckFlag(&flag, "CVodeSensSStolerances", 1));
458
459 // error control strategy should test the sensitivity variables
460 flag = CVodeSetSensErrCon(cvode_mem, true);
461 assert(CheckFlag(&flag, "CVodeSetSensErrCon", 1));
462 }
463
ForwardSensitivityRHS(int Ns,realtype time,N_Vector y,N_Vector ydot,N_Vector * ys,N_Vector * ySdot,void * user_data,N_Vector tmp1,N_Vector tmp2)464 int ODEBase::ForwardSensitivityRHS(int Ns, realtype time, N_Vector y, N_Vector ydot, N_Vector *ys, N_Vector *ySdot, void *user_data, N_Vector tmp1, N_Vector tmp2) {
465 // get the data type
466 ODEData* data = (ODEData*)user_data;
467 assert(data);
468 assert(data->rhs);
469
470 // set the state input
471 #if MUQ_HAS_PARCER==1
472 NVectorMap(stateref, data->comm, y);
473 #else
474 Eigen::Map<Eigen::VectorXd> stateref(NV_DATA_S(y), NV_LENGTH_S(y));
475 #endif
476 data->UpdateInputs(stateref, time);
477
478 for( unsigned int i=0; i<Ns; ++i ) {
479 #if MUQ_HAS_PARCER==1
480 NVectorMap(rhsMap, data->comm, ySdot[i]);
481 NVectorMap(sensMap, data->comm, ys[i]);
482 #else
483 Eigen::Map<Eigen::VectorXd> rhsMap(NV_DATA_S(ySdot[i]), stateref.size());
484 Eigen::Map<Eigen::VectorXd> sensMap(NV_DATA_S(ys[i]), stateref.size());
485 #endif
486
487 rhsMap = data->rhs->ApplyJacobian(0, data->autonomous? 0 : 1, data->inputs, (Eigen::VectorXd)sensMap);
488 }
489
490 // the derivative of the rhs wrt the parameter
491 const Eigen::MatrixXd& dfdp = data->rhs->Jacobian(0, data->autonomous? data->wrtIn : data->wrtIn+1, data->inputs);
492
493 // now, loop through and fill in the rhs vectors stored in ySdot
494 for( unsigned int i=0; i<Ns; ++i ) {
495 #if MUQ_HAS_PARCER==1
496 Eigen::Map<Eigen::VectorXd> rhsMap(data->comm? NV_DATA_P(ySdot[i]) : NV_DATA_S(ySdot[i]), stateref.size());
497 Eigen::Map<Eigen::VectorXd> sensMap(data->comm? NV_DATA_P(ys[i]) : NV_DATA_S(ys[i]), stateref.size());
498 #else
499 Eigen::Map<Eigen::VectorXd> rhsMap(NV_DATA_S(ySdot[i]), stateref.size());
500 Eigen::Map<Eigen::VectorXd> sensMap(NV_DATA_S(ys[i]), stateref.size());
501 #endif
502 rhsMap += dfdp.col(i);
503 }
504
505 return 0;
506 }
507