1 /* -----------------------------------------------------------------------------
2 * Programmer(s): David J. Gardner @ LLNL
3 * -----------------------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2021, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 * -----------------------------------------------------------------------------
14 * This the implementation file for the CVODES nonlinear solver interface.
15 * ---------------------------------------------------------------------------*/
16
17 /*
18 * When sensitivities are computed using the CV_SIMULTANEOUS approach and the
19 * Newton solver is selected the iteraiton is a quasi-Newton method on the
20 * combined system (by approximating the Jacobian matrix by its block diagonal)
21 * and thus only solve linear systems with multiple right hand sides (all
22 * sharing the same coefficient matrix - whatever iteration matrix we decide on)
23 * we set-up the linear solver to handle N equations at a time.
24 */
25
26 #include "cvodes_impl.h"
27 #include "sundials/sundials_math.h"
28 #include "sundials/sundials_nvector_senswrapper.h"
29
30 /* constant macros */
31 #define ONE RCONST(1.0)
32
33 /* private functions */
34 static int cvNlsResidualSensSim(N_Vector ycorSim, N_Vector resSim,
35 void* cvode_mem);
36 static int cvNlsFPFunctionSensSim(N_Vector ycorSim, N_Vector resSim,
37 void* cvode_mem);
38
39 static int cvNlsLSetupSensSim(booleantype jbad, booleantype* jcur,
40 void* cvode_mem);
41 static int cvNlsLSolveSensSim(N_Vector deltaSim, void* cvode_mem);
42 static int cvNlsConvTestSensSim(SUNNonlinearSolver NLS,
43 N_Vector ycorSim, N_Vector delSim,
44 realtype tol, N_Vector ewtSim, void* cvode_mem);
45
46 /* -----------------------------------------------------------------------------
47 * Exported functions
48 * ---------------------------------------------------------------------------*/
49
CVodeSetNonlinearSolverSensSim(void * cvode_mem,SUNNonlinearSolver NLS)50 int CVodeSetNonlinearSolverSensSim(void *cvode_mem, SUNNonlinearSolver NLS)
51 {
52 CVodeMem cv_mem;
53 int retval, is;
54
55 /* Return immediately if CVode memory is NULL */
56 if (cvode_mem == NULL) {
57 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
58 "CVodeSetNonlinearSolverSensSim", MSGCV_NO_MEM);
59 return(CV_MEM_NULL);
60 }
61 cv_mem = (CVodeMem) cvode_mem;
62
63 /* Return immediately if NLS memory is NULL */
64 if (NLS == NULL) {
65 cvProcessError(NULL, CV_ILL_INPUT, "CVODES",
66 "CVodeSetNonlinearSolverSensSim",
67 "NLS must be non-NULL");
68 return (CV_ILL_INPUT);
69 }
70
71 /* check for required nonlinear solver functions */
72 if ( NLS->ops->gettype == NULL ||
73 NLS->ops->solve == NULL ||
74 NLS->ops->setsysfn == NULL ) {
75 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
76 "CVodeSetNonlinearSolverSensSim",
77 "NLS does not support required operations");
78 return(CV_ILL_INPUT);
79 }
80
81 /* check that sensitivities were initialized */
82 if (!(cv_mem->cv_sensi)) {
83 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
84 "CVodeSetNonlinearSolverSensSim",
85 MSGCV_NO_SENSI);
86 return(CV_ILL_INPUT);
87 }
88
89 /* check that simultaneous corrector was selected */
90 if (cv_mem->cv_ism != CV_SIMULTANEOUS) {
91 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
92 "CVodeSetNonlinearSolverSensStg",
93 "Sensitivity solution method is not CV_SIMULTANEOUS");
94 return(CV_ILL_INPUT);
95 }
96
97 /* free any existing nonlinear solver */
98 if ((cv_mem->NLSsim != NULL) && (cv_mem->ownNLSsim))
99 retval = SUNNonlinSolFree(cv_mem->NLSsim);
100
101 /* set SUNNonlinearSolver pointer */
102 cv_mem->NLSsim = NLS;
103
104 /* Set NLS ownership flag. If this function was called to attach the default
105 NLS, CVODE will set the flag to SUNTRUE after this function returns. */
106 cv_mem->ownNLSsim = SUNFALSE;
107
108 /* set the nonlinear system function */
109 if (SUNNonlinSolGetType(NLS) == SUNNONLINEARSOLVER_ROOTFIND) {
110 retval = SUNNonlinSolSetSysFn(cv_mem->NLSsim, cvNlsResidualSensSim);
111 } else if (SUNNonlinSolGetType(NLS) == SUNNONLINEARSOLVER_FIXEDPOINT) {
112 retval = SUNNonlinSolSetSysFn(cv_mem->NLSsim, cvNlsFPFunctionSensSim);
113 } else {
114 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
115 "CVodeSetNonlinearSolverSensSim",
116 "Invalid nonlinear solver type");
117 return(CV_ILL_INPUT);
118 }
119
120 if (retval != CV_SUCCESS) {
121 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
122 "CVodeSetNonlinearSolverSensSim",
123 "Setting nonlinear system function failed");
124 return(CV_ILL_INPUT);
125 }
126
127 /* set convergence test function */
128 retval = SUNNonlinSolSetConvTestFn(cv_mem->NLSsim, cvNlsConvTestSensSim,
129 cvode_mem);
130 if (retval != CV_SUCCESS) {
131 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
132 "CVodeSetNonlinearSolverSensSim",
133 "Setting convergence test function failed");
134 return(CV_ILL_INPUT);
135 }
136
137 /* set max allowed nonlinear iterations */
138 retval = SUNNonlinSolSetMaxIters(cv_mem->NLSsim, NLS_MAXCOR);
139 if (retval != CV_SUCCESS) {
140 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",
141 "CVodeSetNonlinearSolverSensSim",
142 "Setting maximum number of nonlinear iterations failed");
143 return(CV_ILL_INPUT);
144 }
145
146 /* create vector wrappers if necessary */
147 if (cv_mem->simMallocDone == SUNFALSE) {
148
149 cv_mem->zn0Sim = N_VNewEmpty_SensWrapper(cv_mem->cv_Ns+1);
150 if (cv_mem->zn0Sim == NULL) {
151 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES",
152 "CVodeSetNonlinearSolverSensSim", MSGCV_MEM_FAIL);
153 return(CV_MEM_FAIL);
154 }
155
156 cv_mem->ycorSim = N_VNewEmpty_SensWrapper(cv_mem->cv_Ns+1);
157 if (cv_mem->ycorSim == NULL) {
158 N_VDestroy(cv_mem->zn0Sim);
159 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES",
160 "CVodeSetNonlinearSolverSensSim", MSGCV_MEM_FAIL);
161 return(CV_MEM_FAIL);
162 }
163
164 cv_mem->ewtSim = N_VNewEmpty_SensWrapper(cv_mem->cv_Ns+1);
165 if (cv_mem->ewtSim == NULL) {
166 N_VDestroy(cv_mem->zn0Sim);
167 N_VDestroy(cv_mem->ycorSim);
168 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES",
169 "CVodeSetNonlinearSolverSensSim", MSGCV_MEM_FAIL);
170 return(CV_MEM_FAIL);
171 }
172
173 cv_mem->simMallocDone = SUNTRUE;
174 }
175
176 /* attach vectors to vector wrappers */
177 NV_VEC_SW(cv_mem->zn0Sim, 0) = cv_mem->cv_zn[0];
178 NV_VEC_SW(cv_mem->ycorSim, 0) = cv_mem->cv_acor;
179 NV_VEC_SW(cv_mem->ewtSim, 0) = cv_mem->cv_ewt;
180
181 for (is=0; is < cv_mem->cv_Ns; is++) {
182 NV_VEC_SW(cv_mem->zn0Sim, is+1) = cv_mem->cv_znS[0][is];
183 NV_VEC_SW(cv_mem->ycorSim, is+1) = cv_mem->cv_acorS[is];
184 NV_VEC_SW(cv_mem->ewtSim, is+1) = cv_mem->cv_ewtS[is];
185 }
186
187 /* Reset the acnrmcur flag to SUNFALSE */
188 cv_mem->cv_acnrmcur = SUNFALSE;
189
190 return(CV_SUCCESS);
191 }
192
193
194 /*---------------------------------------------------------------
195 CVodeGetNonlinearSystemDataSens:
196
197 This routine provides access to the relevant data needed to
198 compute the nonlinear system function.
199 ---------------------------------------------------------------*/
CVodeGetNonlinearSystemDataSens(void * cvode_mem,realtype * tcur,N_Vector ** ySpred,N_Vector ** ySn,realtype * gamma,realtype * rl1,N_Vector ** znS1,void ** user_data)200 int CVodeGetNonlinearSystemDataSens(void *cvode_mem, realtype *tcur,
201 N_Vector **ySpred, N_Vector **ySn,
202 realtype *gamma, realtype *rl1,
203 N_Vector **znS1, void **user_data)
204 {
205 CVodeMem cv_mem;
206
207 if (cvode_mem == NULL) {
208 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
209 "CVodeGetNonlinearSystemDataSens", MSGCV_NO_MEM);
210 return(CV_MEM_NULL);
211 }
212
213 cv_mem = (CVodeMem) cvode_mem;
214
215 *tcur = cv_mem->cv_tn;
216 *ySpred = cv_mem->cv_znS[0];
217 *ySn = cv_mem->cv_yS;
218 *gamma = cv_mem->cv_gamma;
219 *rl1 = cv_mem->cv_rl1;
220 *znS1 = cv_mem->cv_znS[1];
221 *user_data = cv_mem->cv_user_data;
222
223 return(CV_SUCCESS);
224 }
225
226
227 /* -----------------------------------------------------------------------------
228 * Private functions
229 * ---------------------------------------------------------------------------*/
230
231
cvNlsInitSensSim(CVodeMem cvode_mem)232 int cvNlsInitSensSim(CVodeMem cvode_mem)
233 {
234 int retval;
235
236 /* set the linear solver setup wrapper function */
237 if (cvode_mem->cv_lsetup)
238 retval = SUNNonlinSolSetLSetupFn(cvode_mem->NLSsim, cvNlsLSetupSensSim);
239 else
240 retval = SUNNonlinSolSetLSetupFn(cvode_mem->NLSsim, NULL);
241
242 if (retval != CV_SUCCESS) {
243 cvProcessError(cvode_mem, CV_ILL_INPUT, "CVODES", "cvNlsInitSensSim",
244 "Setting the linear solver setup function failed");
245 return(CV_NLS_INIT_FAIL);
246 }
247
248 /* set the linear solver solve wrapper function */
249 if (cvode_mem->cv_lsolve)
250 retval = SUNNonlinSolSetLSolveFn(cvode_mem->NLSsim, cvNlsLSolveSensSim);
251 else
252 retval = SUNNonlinSolSetLSolveFn(cvode_mem->NLSsim, NULL);
253
254 if (retval != CV_SUCCESS) {
255 cvProcessError(cvode_mem, CV_ILL_INPUT, "CVODES", "cvNlsInitSensSim",
256 "Setting linear solver solve function failed");
257 return(CV_NLS_INIT_FAIL);
258 }
259
260 /* initialize nonlinear solver */
261 retval = SUNNonlinSolInitialize(cvode_mem->NLSsim);
262
263 if (retval != CV_SUCCESS) {
264 cvProcessError(cvode_mem, CV_ILL_INPUT, "CVODES", "cvNlsInitSensSim",
265 MSGCV_NLS_INIT_FAIL);
266 return(CV_NLS_INIT_FAIL);
267 }
268
269 return(CV_SUCCESS);
270 }
271
272
cvNlsLSetupSensSim(booleantype jbad,booleantype * jcur,void * cvode_mem)273 static int cvNlsLSetupSensSim(booleantype jbad, booleantype* jcur,
274 void* cvode_mem)
275 {
276 CVodeMem cv_mem;
277 int retval;
278
279 if (cvode_mem == NULL) {
280 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
281 "cvNlsLSetupSensSim", MSGCV_NO_MEM);
282 return(CV_MEM_NULL);
283 }
284 cv_mem = (CVodeMem) cvode_mem;
285
286 /* if the nonlinear solver marked the Jacobian as bad update convfail */
287 if (jbad)
288 cv_mem->convfail = CV_FAIL_BAD_J;
289
290 /* setup the linear solver */
291 retval = cv_mem->cv_lsetup(cv_mem, cv_mem->convfail, cv_mem->cv_y,
292 cv_mem->cv_ftemp, &(cv_mem->cv_jcur),
293 cv_mem->cv_vtemp1, cv_mem->cv_vtemp2,
294 cv_mem->cv_vtemp3);
295 cv_mem->cv_nsetups++;
296
297 /* update Jacobian status */
298 *jcur = cv_mem->cv_jcur;
299
300 cv_mem->cv_forceSetup = SUNFALSE;
301 cv_mem->cv_gamrat = ONE;
302 cv_mem->cv_gammap = cv_mem->cv_gamma;
303 cv_mem->cv_crate = ONE;
304 cv_mem->cv_crateS = ONE;
305 cv_mem->cv_nstlp = cv_mem->cv_nst;
306
307 if (retval < 0) return(CV_LSETUP_FAIL);
308 if (retval > 0) return(SUN_NLS_CONV_RECVR);
309
310 return(CV_SUCCESS);
311 }
312
313
cvNlsLSolveSensSim(N_Vector deltaSim,void * cvode_mem)314 static int cvNlsLSolveSensSim(N_Vector deltaSim, void* cvode_mem)
315 {
316 CVodeMem cv_mem;
317 int retval, is;
318 N_Vector delta;
319 N_Vector *deltaS;
320
321 if (cvode_mem == NULL) {
322 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
323 "cvNlsLSolveSensSim", MSGCV_NO_MEM);
324 return(CV_MEM_NULL);
325 }
326 cv_mem = (CVodeMem) cvode_mem;
327
328 /* extract state delta from the vector wrapper */
329 delta = NV_VEC_SW(deltaSim,0);
330
331 /* solve the state linear system */
332 retval = cv_mem->cv_lsolve(cv_mem, delta, cv_mem->cv_ewt, cv_mem->cv_y,
333 cv_mem->cv_ftemp);
334
335 if (retval < 0) return(CV_LSOLVE_FAIL);
336 if (retval > 0) return(SUN_NLS_CONV_RECVR);
337
338 /* extract sensitivity deltas from the vector wrapper */
339 deltaS = NV_VECS_SW(deltaSim)+1;
340
341 /* solve the sensitivity linear systems */
342 for (is=0; is<cv_mem->cv_Ns; is++) {
343 retval = cv_mem->cv_lsolve(cv_mem, deltaS[is], cv_mem->cv_ewtS[is],
344 cv_mem->cv_y, cv_mem->cv_ftemp);
345
346 if (retval < 0) return(CV_LSOLVE_FAIL);
347 if (retval > 0) return(SUN_NLS_CONV_RECVR);
348 }
349
350 return(CV_SUCCESS);
351 }
352
353
cvNlsConvTestSensSim(SUNNonlinearSolver NLS,N_Vector ycorSim,N_Vector deltaSim,realtype tol,N_Vector ewtSim,void * cvode_mem)354 static int cvNlsConvTestSensSim(SUNNonlinearSolver NLS,
355 N_Vector ycorSim, N_Vector deltaSim,
356 realtype tol, N_Vector ewtSim, void* cvode_mem)
357 {
358 CVodeMem cv_mem;
359 int m, retval;
360 realtype del, delS, Del;
361 realtype dcon;
362 N_Vector ycor, delta, ewt;
363 N_Vector *deltaS, *ewtS;
364
365 if (cvode_mem == NULL) {
366 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
367 "cvNlsConvTestSensSim", MSGCV_NO_MEM);
368 return(CV_MEM_NULL);
369 }
370 cv_mem = (CVodeMem) cvode_mem;
371
372 /* extract the current state and sensitivity corrections */
373 ycor = NV_VEC_SW(ycorSim,0);
374
375 /* extract state and sensitivity deltas */
376 delta = NV_VEC_SW(deltaSim,0);
377 deltaS = NV_VECS_SW(deltaSim)+1;
378
379 /* extract state and sensitivity error weights */
380 ewt = NV_VEC_SW(ewtSim,0);
381 ewtS = NV_VECS_SW(ewtSim)+1;
382
383 /* compute the norm of the state and sensitivity corrections */
384 del = N_VWrmsNorm(delta, ewt);
385 delS = cvSensUpdateNorm(cv_mem, del, deltaS, ewtS);
386
387 /* norm used in error test */
388 Del = delS;
389
390 /* get the current nonlinear solver iteration count */
391 retval = SUNNonlinSolGetCurIter(NLS, &m);
392 if (retval != CV_SUCCESS) return(CV_MEM_NULL);
393
394 /* Test for convergence. If m > 0, an estimate of the convergence
395 rate constant is stored in crate, and used in the test.
396
397 Recall that, even when errconS=SUNFALSE, all variables are used in the
398 convergence test. Hence, we use Del (and not del). However, acnrm is used
399 in the error test and thus it has different forms depending on errconS
400 (and this explains why we have to carry around del and delS).
401 */
402 if (m > 0) {
403 cv_mem->cv_crate = SUNMAX(CRDOWN * cv_mem->cv_crate, Del/cv_mem->cv_delp);
404 }
405 dcon = Del * SUNMIN(ONE, cv_mem->cv_crate) / tol;
406
407 /* check if nonlinear system was solved successfully */
408 if (dcon <= ONE) {
409 if (m == 0) {
410 cv_mem->cv_acnrm = (cv_mem->cv_errconS) ? delS : del;
411 } else {
412 cv_mem->cv_acnrm = (cv_mem->cv_errconS) ?
413 N_VWrmsNorm(ycorSim, ewtSim) : N_VWrmsNorm(ycor, ewt);
414 }
415 cv_mem->cv_acnrmcur = SUNTRUE;
416 return(CV_SUCCESS);
417 }
418
419 /* check if the iteration seems to be diverging */
420 if ((m >= 1) && (Del > RDIV*cv_mem->cv_delp)) return(SUN_NLS_CONV_RECVR);
421
422 /* Save norm of correction and loop again */
423 cv_mem->cv_delp = Del;
424
425 /* Not yet converged */
426 return(SUN_NLS_CONTINUE);
427 }
428
429
cvNlsResidualSensSim(N_Vector ycorSim,N_Vector resSim,void * cvode_mem)430 static int cvNlsResidualSensSim(N_Vector ycorSim, N_Vector resSim, void* cvode_mem)
431 {
432 CVodeMem cv_mem;
433 int retval;
434 N_Vector ycor, res;
435 N_Vector *ycorS, *resS;
436 realtype cvals[3];
437 N_Vector* XXvecs[3];
438
439 if (cvode_mem == NULL) {
440 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
441 "cvNlsResidualSensSim", MSGCV_NO_MEM);
442 return(CV_MEM_NULL);
443 }
444 cv_mem = (CVodeMem) cvode_mem;
445
446 /* extract state and residual vectors from the vector wrapper */
447 ycor = NV_VEC_SW(ycorSim,0);
448 res = NV_VEC_SW(resSim,0);
449
450 /* update the state based on the current correction */
451 N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, ycor, cv_mem->cv_y);
452
453 /* evaluate the rhs function */
454 retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y, cv_mem->cv_ftemp,
455 cv_mem->cv_user_data);
456 cv_mem->cv_nfe++;
457 if (retval < 0) return(CV_RHSFUNC_FAIL);
458 if (retval > 0) return(RHSFUNC_RECVR);
459
460 /* compute the resiudal */
461 N_VLinearSum(cv_mem->cv_rl1, cv_mem->cv_zn[1], ONE, ycor, res);
462 N_VLinearSum(-cv_mem->cv_gamma, cv_mem->cv_ftemp, ONE, res, res);
463
464 /* extract sensitivity and residual vectors from the vector wrapper */
465 ycorS = NV_VECS_SW(ycorSim)+1;
466 resS = NV_VECS_SW(resSim)+1;
467
468 /* update sensitivities based on the current correction */
469 retval = N_VLinearSumVectorArray(cv_mem->cv_Ns,
470 ONE, cv_mem->cv_znS[0],
471 ONE, ycorS, cv_mem->cv_yS);
472 if (retval != CV_SUCCESS) return(CV_VECTOROP_ERR);
473
474 /* evaluate the sensitivity rhs function */
475 retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn,
476 cv_mem->cv_y, cv_mem->cv_ftemp,
477 cv_mem->cv_yS, cv_mem->cv_ftempS,
478 cv_mem->cv_vtemp1, cv_mem->cv_vtemp2);
479
480 if (retval < 0) return(CV_SRHSFUNC_FAIL);
481 if (retval > 0) return(SRHSFUNC_RECVR);
482
483 /* compute the sensitivity resiudal */
484 cvals[0] = cv_mem->cv_rl1; XXvecs[0] = cv_mem->cv_znS[1];
485 cvals[1] = ONE; XXvecs[1] = ycorS;
486 cvals[2] = -cv_mem->cv_gamma; XXvecs[2] = cv_mem->cv_ftempS;
487
488 retval = N_VLinearCombinationVectorArray(cv_mem->cv_Ns,
489 3, cvals, XXvecs, resS);
490 if (retval != CV_SUCCESS) return(CV_VECTOROP_ERR);
491
492 return(CV_SUCCESS);
493 }
494
495
cvNlsFPFunctionSensSim(N_Vector ycorSim,N_Vector resSim,void * cvode_mem)496 static int cvNlsFPFunctionSensSim(N_Vector ycorSim, N_Vector resSim, void* cvode_mem)
497 {
498 CVodeMem cv_mem;
499 int retval, is;
500 N_Vector ycor, res;
501 N_Vector *ycorS, *resS;
502
503 if (cvode_mem == NULL) {
504 cvProcessError(NULL, CV_MEM_NULL, "CVODES",
505 "cvNlsFPFunctionSensSim", MSGCV_NO_MEM);
506 return(CV_MEM_NULL);
507 }
508 cv_mem = (CVodeMem) cvode_mem;
509
510 /* extract state and residual vectors from the vector wrapper */
511 ycor = NV_VEC_SW(ycorSim,0);
512 res = NV_VEC_SW(resSim,0);
513
514 /* update the state based on the current correction */
515 N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, ycor, cv_mem->cv_y);
516
517 /* evaluate the rhs function */
518 retval = cv_mem->cv_f(cv_mem->cv_tn, cv_mem->cv_y, res,
519 cv_mem->cv_user_data);
520 cv_mem->cv_nfe++;
521 if (retval < 0) return(CV_RHSFUNC_FAIL);
522 if (retval > 0) return(RHSFUNC_RECVR);
523
524 /* evaluate fixed point function */
525 N_VLinearSum(cv_mem->cv_h, res, -ONE, cv_mem->cv_zn[1], res);
526 N_VScale(cv_mem->cv_rl1, res, res);
527
528 /* extract sensitivity and residual vectors from the vector wrapper */
529 ycorS = NV_VECS_SW(ycorSim)+1;
530 resS = NV_VECS_SW(resSim)+1;
531
532 /* update the sensitivities based on the current correction */
533 N_VLinearSumVectorArray(cv_mem->cv_Ns,
534 ONE, cv_mem->cv_znS[0],
535 ONE, ycorS, cv_mem->cv_yS);
536
537 /* evaluate the sensitivity rhs function */
538 retval = cvSensRhsWrapper(cv_mem, cv_mem->cv_tn,
539 cv_mem->cv_y, res,
540 cv_mem->cv_yS, resS,
541 cv_mem->cv_vtemp1, cv_mem->cv_vtemp2);
542
543 if (retval < 0) return(CV_SRHSFUNC_FAIL);
544 if (retval > 0) return(SRHSFUNC_RECVR);
545
546 /* evaluate sensitivity fixed point function */
547 for (is=0; is<cv_mem->cv_Ns; is++) {
548 N_VLinearSum(cv_mem->cv_h, resS[is], -ONE, cv_mem->cv_znS[1][is], resS[is]);
549 N_VScale(cv_mem->cv_rl1, resS[is], resS[is]);
550 }
551
552 return(CV_SUCCESS);
553 }
554