diff --git a/CHANGELOG.md b/CHANGELOG.md index 22216ca69a..f30f2da4c1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,8 @@ stage is computed (`ARKodeSetPreprocessStepFn`, `ARKodeSetPostprocessStepFn`, should treat the state vector as read-only, otherwise all theoretical guarantees of solution accuracy and stability will be lost. +Removed extraneous copy of output vector when using ARKODE in ``ARK_ONE_STEP`` mode. + ### Bug Fixes Fixed a CMake bug where the SuperLU_MT interface would not be built and diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index 9194c08ced..749f045051 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -15,6 +15,8 @@ These are considered **advanced** functions, as they should treat the state vect read-only, otherwise all theoretical guarantees of solution accuracy and stability will be lost. +Removed extraneous copy of output vector when using ARKODE in ``ARK_ONE_STEP`` mode. + **Bug Fixes** Fixed a CMake bug where the SuperLU_MT interface would not be built and diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index 2793ac89f9..36d5759d45 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -1148,13 +1148,12 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout, break; } - /* In ONE_STEP mode, copy y and exit loop */ + /* In ONE_STEP mode, exit loop (arkCompleteStep already copied yn to ycur, an alias to yout) */ if (itask == ARK_ONE_STEP) { istate = ARK_SUCCESS; ark_mem->tretlast = *tret = ark_mem->tcur; - N_VScale(ONE, ark_mem->yn, yout); - ark_mem->next_h = ark_mem->hprime; + ark_mem->next_h = ark_mem->hprime; break; } diff --git a/src/arkode/arkode_root.c b/src/arkode/arkode_root.c index 5d621987ad..5714e87de0 100644 --- a/src/arkode/arkode_root.c +++ b/src/arkode/arkode_root.c @@ -471,8 +471,9 @@ int arkRootCheck1(void* arkode_mem) hratio = SUNMAX(rootmem->ttol / SUNRabs(ark_mem->h), TENTH); smallh = hratio * ark_mem->h; tplus = rootmem->tlo + smallh; - N_VLinearSum(ONE, ark_mem->yn, smallh, ark_mem->fn, ark_mem->ycur); - retval = rootmem->gfun(tplus, ark_mem->ycur, rootmem->ghi, rootmem->root_data); + N_VLinearSum(ONE, ark_mem->yn, smallh, ark_mem->fn, ark_mem->tempv4); + retval = rootmem->gfun(tplus, ark_mem->tempv4, rootmem->ghi, + rootmem->root_data); rootmem->nge++; if (retval != 0) { @@ -533,11 +534,11 @@ int arkRootCheck2(void* arkode_mem) /* return if no roots in previous step */ if (rootmem->irfnd == 0) { return (ARK_SUCCESS); } - /* Set ark_ycur = y(tlo) */ - (void)ARKodeGetDky(ark_mem, rootmem->tlo, 0, ark_mem->ycur); + /* Set tempv4 = y(tlo) */ + (void)ARKodeGetDky(ark_mem, rootmem->tlo, 0, ark_mem->tempv4); /* Evaluate root-finding function: glo = g(tlo, y(tlo)) */ - retval = rootmem->gfun(rootmem->tlo, ark_mem->ycur, rootmem->glo, + retval = rootmem->gfun(rootmem->tlo, ark_mem->tempv4, rootmem->glo, rootmem->root_data); rootmem->nge++; if (retval != 0) { return (ARK_RTFUNC_FAIL); } @@ -569,7 +570,7 @@ int arkRootCheck2(void* arkode_mem) if ((tplus - ark_mem->tcur) * ark_mem->h >= ZERO) { /* hratio = smallh/ark_mem->h; */ - N_VLinearSum(ONE, ark_mem->ycur, smallh, ark_mem->fn, ark_mem->ycur); + N_VLinearSum(ONE, ark_mem->tempv4, smallh, ark_mem->fn, ark_mem->ycur); } else { @@ -632,24 +633,24 @@ int arkRootCheck3(void* arkode_mem, sunrealtype tout, int itask) if (itask == ARK_ONE_STEP) { rootmem->thi = ark_mem->tcur; - N_VScale(ONE, ark_mem->yn, ark_mem->ycur); + N_VScale(ONE, ark_mem->yn, ark_mem->tempv4); } if (itask == ARK_NORMAL) { if ((tout - ark_mem->tcur) * ark_mem->h >= ZERO) { rootmem->thi = ark_mem->tcur; - N_VScale(ONE, ark_mem->yn, ark_mem->ycur); + N_VScale(ONE, ark_mem->yn, ark_mem->tempv4); } else { rootmem->thi = tout; - (void)ARKodeGetDky(ark_mem, rootmem->thi, 0, ark_mem->ycur); + (void)ARKodeGetDky(ark_mem, rootmem->thi, 0, ark_mem->tempv4); } } /* Set rootmem->ghi = g(thi) and call arkRootfind to search (tlo,thi) for roots. */ - retval = rootmem->gfun(rootmem->thi, ark_mem->ycur, rootmem->ghi, + retval = rootmem->gfun(rootmem->thi, ark_mem->tempv4, rootmem->ghi, rootmem->root_data); rootmem->nge++; if (retval != 0) { return (ARK_RTFUNC_FAIL); } @@ -861,8 +862,8 @@ int arkRootfind(void* arkode_mem) tmid = rootmem->thi - fracsub * (rootmem->thi - rootmem->tlo); } - (void)ARKodeGetDky(ark_mem, tmid, 0, ark_mem->ycur); - retval = rootmem->gfun(tmid, ark_mem->ycur, rootmem->grout, + (void)ARKodeGetDky(ark_mem, tmid, 0, ark_mem->tempv4); + retval = rootmem->gfun(tmid, ark_mem->tempv4, rootmem->grout, rootmem->root_data); rootmem->nge++; if (retval != 0) { return (ARK_RTFUNC_FAIL); }