| 1209 | *****************************************************************/ |
| 1210 | |
| 1211 | static boole CVHin(CVodeMem cv_mem, real tout) |
| 1212 | { |
| 1213 | int sign, count; |
| 1214 | real tdiff, tdist, tround, hlb, hub; |
| 1215 | real hg, hgs, hnew, hrat, h0, yddnrm; |
| 1216 | |
| 1217 | /* Test for tout too close to tn */ |
| 1218 | |
| 1219 | if ((tdiff = tout-tn) == ZERO) return(FALSE); |
| 1220 | |
| 1221 | sign = (tdiff > ZERO) ? 1 : -1; |
| 1222 | tdist = ABS(tdiff); |
| 1223 | tround = uround * MAX(ABS(tn), ABS(tout)); |
| 1224 | if (tdist < TWO*tround) return(FALSE); |
| 1225 | |
| 1226 | /* Set lower and upper bounds on h0, and take geometric mean |
| 1227 | Exit with this value if the bounds cross each other */ |
| 1228 | |
| 1229 | hlb = HLB_FACTOR * tround; |
| 1230 | hub = CVUpperBoundH0(cv_mem, tdist); |
| 1231 | hg = RSqrt(hlb*hub); |
| 1232 | if (hub < hlb) { |
| 1233 | if (sign == -1) hg = -hg; |
| 1234 | h = hg; |
| 1235 | return(TRUE); |
| 1236 | } |
| 1237 | |
| 1238 | /* Loop up to MAX_ITERS times to find h0. |
| 1239 | Stop if new and previous values differ by a factor < 2. |
| 1240 | Stop if hnew/hg > 2 after one iteration, as this probably means |
| 1241 | that the ydd value is bad because of cancellation error. */ |
| 1242 | |
| 1243 | count = 0; |
| 1244 | loop { |
| 1245 | hgs = hg*sign; |
| 1246 | yddnrm = CVYddNorm(cv_mem, hgs); |
| 1247 | hnew = (yddnrm*hub*hub > TWO) ? RSqrt(TWO/yddnrm) : RSqrt(hg*hub); |
| 1248 | count++; |
| 1249 | if (count >= MAX_ITERS) break; |
| 1250 | hrat = hnew/hg; |
| 1251 | if ((hrat > HALF) && (hrat < TWO)) break; |
| 1252 | if ((count >= 2) && (hrat > TWO)) { |
| 1253 | hnew = hg; |
| 1254 | break; |
| 1255 | } |
| 1256 | hg = hnew; |
| 1257 | } |
| 1258 | |
| 1259 | /* Apply bounds, bias factor, and attach sign */ |
| 1260 |
no test coverage detected