| 193 | **********************************************************************/ |
| 194 | |
| 195 | static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred, |
| 196 | N_Vector fpred, boole *jcurPtr, N_Vector vtemp1, |
| 197 | N_Vector vtemp2, N_Vector vtemp3) |
| 198 | { |
| 199 | real r; |
| 200 | N_Vector ftemp, y; |
| 201 | boole invOK; |
| 202 | CVDiagMem cvdiag_mem; |
| 203 | |
| 204 | cvdiag_mem = (CVDiagMem) lmem; |
| 205 | |
| 206 | /* Rename work vectors for use as temporary values of y and f */ |
| 207 | ftemp = vtemp1; |
| 208 | y = vtemp2; |
| 209 | |
| 210 | /* Form y with perturbation = FRACT*(func. iter. correction) */ |
| 211 | r = FRACT * rl1; |
| 212 | N_VLinearSum(h, fpred, -ONE, zn[1], ftemp); |
| 213 | N_VLinearSum(r, ftemp, ONE, ypred, y); |
| 214 | |
| 215 | /* Evaluate f at perturbed y */ |
| 216 | f(N, tn, y, M, f_data); |
| 217 | nfe++; |
| 218 | |
| 219 | /* Construct M = I - gamma*J with J = diag(deltaf_i/deltay_i) */ |
| 220 | N_VLinearSum(ONE, M, -ONE, fpred, M); |
| 221 | N_VLinearSum(FRACT, ftemp, -h, M, M); |
| 222 | N_VProd(ftemp, ewt, y); |
| 223 | /* Protect against deltay_i being at roundoff level */ |
| 224 | N_VCompare(uround, y, bit); |
| 225 | N_VAddConst(bit, -ONE, bitcomp); |
| 226 | N_VProd(ftemp, bit, y); |
| 227 | N_VLinearSum(FRACT, y, -ONE, bitcomp, y); |
| 228 | N_VDiv(M, y, M); |
| 229 | N_VProd(M, bit, M); |
| 230 | N_VLinearSum(ONE, M, -ONE, bitcomp, M); |
| 231 | |
| 232 | /* Invert M with test for zero components */ |
| 233 | invOK = N_VInvTest(M, M); |
| 234 | if (!invOK) return(1); |
| 235 | |
| 236 | /* Set jcur = TRUE, save gamma in gammasv, and return */ |
| 237 | *jcurPtr = TRUE; |
| 238 | gammasv = gamma; |
| 239 | return(0); |
| 240 | } |
| 241 | |
| 242 | /*************** CVDiagSolve ***************************************** |
| 243 |
nothing calls this directly
no test coverage detected