diff --git a/.ci-pipelines/ci-common-defs.sh b/.ci-pipelines/ci-common-defs.sh index 8f2c730..070861a 100755 --- a/.ci-pipelines/ci-common-defs.sh +++ b/.ci-pipelines/ci-common-defs.sh @@ -16,6 +16,7 @@ F90_feuler F90_lsode F90_radau F90_rk +F90_rkadj F90_rktlm F90_ros F90_rosadj @@ -27,7 +28,9 @@ F90_rostlm F90_ros_upcase F90_saprc_2006 F90_sd +F90_sd4 F90_sdadj +F90_sdtlm F90_seulex F90_small_strato " diff --git a/CHANGELOG.md b/CHANGELOG.md index 4589ef7..4dc0e47 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added GitHub Action to run C-I tests with GCC compilers v9, v10, v11, v12, and v13 - Added "Lint" GitHub Action to check other actions for security issues +- Added new example files: `rkadj.kpp`, `sd4.kpp`, `sdtlm.kpp` +- Added new C-I tests: `F90_rkadj`, `F90_sd4`, `F90_sdtlm` ### Changed - Updated ReadTheDocs documentation to reflect that C-I tests are now done as a GitHub Action @@ -23,9 +25,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Moved the `which kpp` instruction to the end of the "Build the KPP executable" section in the installation guide on ReadTheDocs - Updated rules to ignore files in `.gitignore` and updated comments accordingly - Fixed a bug that prevented `.ci-pipelines/ci-cleanup-script.sh` from removing KPP-generated files for MCM mechanisms +- Fixed typo in error message in `int/rosenbrock_autoreduce.f90` ### Removed - Removed C-I tests on Microsoft Azure Dev Pipelines +- Replaced BLAS functions (`WAXPY`, `WCOPY`, `WSCAL`, `WADD`, `WLAMCH`) with pure F90 code from `int/*.f90` integrators (thanks to AI for the help) ## [3.3.0] - 2025-07-17 ### Added diff --git a/ci-tests/F90_rkadj/F90_rkadj.kpp b/ci-tests/F90_rkadj/F90_rkadj.kpp new file mode 120000 index 0000000..b3bf812 --- /dev/null +++ b/ci-tests/F90_rkadj/F90_rkadj.kpp @@ -0,0 +1 @@ +../../examples/rkadj.kpp \ No newline at end of file diff --git a/ci-tests/F90_sd4/F90_sd4.kpp b/ci-tests/F90_sd4/F90_sd4.kpp new file mode 120000 index 0000000..f11975e --- /dev/null +++ b/ci-tests/F90_sd4/F90_sd4.kpp @@ -0,0 +1 @@ +../../examples/sd4.kpp \ No newline at end of file diff --git a/ci-tests/F90_sdtlm/F90_sdtlm.kpp b/ci-tests/F90_sdtlm/F90_sdtlm.kpp new file mode 120000 index 0000000..eb0b88c --- /dev/null +++ b/ci-tests/F90_sdtlm/F90_sdtlm.kpp @@ -0,0 +1 @@ +../../examples/sdtlm.kpp \ No newline at end of file diff --git a/docs/source/tech_info/06_info_for_kpp_developers.rst b/docs/source/tech_info/06_info_for_kpp_developers.rst index cb6798e..7a1f2cc 100644 --- a/docs/source/tech_info/06_info_for_kpp_developers.rst +++ b/docs/source/tech_info/06_info_for_kpp_developers.rst @@ -415,6 +415,10 @@ List of continuous integration tests - Fortran90 - small_strato - runge_kutta + * - :code:`F90_rkadj` + - Fortran90 + - small_strato + - runge_kutta_adj * - :code:`F90_rktlm` - Fortran90 - small_strato @@ -455,6 +459,10 @@ List of continuous integration tests - Fortran90 - saprcnov - rosenbrock + * - :code:`F90_sd4` + - Fortran90 + - small_strato + - sdirk4 * - :code:`F90_sd` - Fortran90 - small_strato @@ -463,6 +471,10 @@ List of continuous integration tests - Fortran90 - small_strato - sdirk_adj + * - :code:`F90_sdtlm` + - Fortran90 + - small_strato + - sdirk_tlm * - :code:`F90_seulex` - Fortran90 - saprcnov diff --git a/examples/rkadj.kpp b/examples/rkadj.kpp new file mode 100644 index 0000000..3e28e80 --- /dev/null +++ b/examples/rkadj.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR runge_kutta_adj +#LANGUAGE Fortran90 +#DRIVER general_adj diff --git a/examples/sd4.kpp b/examples/sd4.kpp new file mode 100644 index 0000000..7dc9ea9 --- /dev/null +++ b/examples/sd4.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR sdirk4 +#LANGUAGE Fortran90 +#DRIVER general diff --git a/examples/sdtlm.kpp b/examples/sdtlm.kpp new file mode 100644 index 0000000..9865d6b --- /dev/null +++ b/examples/sdtlm.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR sdirk_tlm +#LANGUAGE Fortran90 +#DRIVER general_tlm diff --git a/int/dvode.f90 b/int/dvode.f90 index 3e5f7ec..64430e4 100644 --- a/int/dvode.f90 +++ b/int/dvode.f90 @@ -4,7 +4,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP - USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, Set2zero, WLAMCH + USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve IMPLICIT NONE PUBLIC diff --git a/int/lsode.f90 b/int/lsode.f90 index 93a4bb8..3b6d73f 100644 --- a/int/lsode.f90 +++ b/int/lsode.f90 @@ -11,7 +11,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, WLAMCH + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC diff --git a/int/radau5.f90 b/int/radau5.f90 index 89a4125..f409dde 100644 --- a/int/radau5.f90 +++ b/int/radau5.f90 @@ -412,7 +412,7 @@ SUBROUTINE RADAU5(N,T,Tend,Y,H,RelTol,AbsTol, & !~~~> Roundoff SMALLEST NUMBER SATISFYING 1.0d0+Roundoff>1.0d0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> RCNTRL(1) = Hmin - not used Hmin = ZERO @@ -642,12 +642,12 @@ SUBROUTINE RAD_Integrator( & ! STARTING VALUES FOR NEWTON ITERATION !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) - CALL Set2zero(N,F1) - CALL Set2zero(N,F2) - CALL Set2zero(N,F3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp + F1(1:N) = 0.0_dp + F2(1:N) = 0.0_dp + F3(1:N) = 0.0_dp ELSE C3Q=H/Hold C1Q=rkC(1)*C3Q @@ -726,9 +726,9 @@ SUBROUTINE RAD_Integrator( & END IF END IF DYNOLD=MAX(DYNO,Roundoff) - CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1 - CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2 - CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3 + F1(1:N) = F1(1:N) + Z1(1:N) ! F1 <- F1 + Z1 + F2(1:N) = F2(1:N) + Z2(1:N) ! F2 <- F2 + Z2 + F3(1:N) = F3(1:N) + Z3(1:N) ! F3 <- F3 + Z3 ! Z(1,2,3) = Transf x F(1,2,3) CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3) NewtonDone = (FacConv*DYNO <= TolNewton) diff --git a/int/rosenbrock.f90 b/int/rosenbrock.f90 index e5ae740..2460c8d 100644 --- a/int/rosenbrock.f90 +++ b/int/rosenbrock.f90 @@ -335,7 +335,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -519,9 +519,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -581,55 +578,49 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & RETURN END IF -!~~~> Compute the stages -Stage: DO istage = 1, ros_S + !~~~> Compute the stages + Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N) - ioffset = N*(istage-1) + ioffset = N*(istage-1) ! For the 1st istage the function has been computed previously - IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) + IF ( istage == 1 ) THEN Fcn(1:N) = Fcn0(1:N) - ! istage>1 and a new function evaluation is needed at the current istage - ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) + ! istage>1 and a new function evaluation is needed at the current istage + ELSEIF ( ros_NewF(istage) ) THEN Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) & + + ros_A((istage-1)*(istage-2)/2+j) * K(N*(j-1)+1:N*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) - K(ioffset+1:ioffset+N) = Fcn(1:N) - DO j = 1, istage-1 + END IF ! if istage == 1 elseif ros_NewF(istage) + K(ioffset+1:ioffset+N) = Fcn(1:N) + DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) - END DO - IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HC * K(N*(j-1)+1:N*j) + END DO + IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) - END IF - CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HG * dFdT(1:N) + END IF + CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) END DO Stage - !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -645,7 +636,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -732,8 +722,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -781,16 +771,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_adj.f90 b/int/rosenbrock_adj.f90 index 90c6e2b..aa644f9 100644 --- a/int/rosenbrock_adj.f90 +++ b/int/rosenbrock_adj.f90 @@ -438,7 +438,7 @@ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, & !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -747,8 +747,6 @@ SUBROUTINE ros_DPush( S, T, H, Ystage, K, E, P ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - !CALL WCOPY(NVAR*S,Ystage,1,chk_Y(1,stack_ptr),1) - !CALL WCOPY(NVAR*S,K,1,chk_K(1,stack_ptr),1) chk_Y(1:NVAR*S,stack_ptr) = Ystage(1:NVAR*S) chk_K(1:NVAR*S,stack_ptr) = K(1:NVAR*S) IF (SaveLU) THEN @@ -784,11 +782,8 @@ SUBROUTINE ros_DPop( S, T, H, Ystage, K, E, P ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - !CALL WCOPY(NVAR*S,chk_Y(1,stack_ptr),1,Ystage,1) - !CALL WCOPY(NVAR*S,chk_K(1,stack_ptr),1,K,1) Ystage(1:NVAR*S) = chk_Y(1:NVAR*S,stack_ptr) K(1:NVAR*S) = chk_K(1:NVAR*S,stack_ptr) - !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1) IF (SaveLU) THEN #ifdef FULL_ALGEBRA E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr) @@ -817,9 +812,6 @@ SUBROUTINE ros_CPush( T, H, Y, dY, d2Y ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - !CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - !CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1) - !CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1) chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR) chk_d2Y(1:NVAR,stack_ptr) = d2Y(1:NVAR) @@ -840,9 +832,6 @@ SUBROUTINE ros_CPop( T, H, Y, dY, d2Y ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - !CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) - !CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1) - !CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr) d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr) @@ -1004,60 +993,59 @@ SUBROUTINE ros_FwdInt ( Y, & RETURN END IF -!~~~> Compute the stages -Stage: DO istage = 1, ros_S + !~~~> Compute the stages + Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR) - ioffset = NVAR*(istage-1) + ioffset = NVAR*(istage-1) ! For the 1st istage the function has been computed previously - IF ( istage == 1 ) THEN - CALL WCOPY(NVAR,Fcn0,1,Fcn,1) + IF ( istage == 1 ) THEN + Fcn(1:NVAR) = Fcn0(1:NVAR) IF (AdjointType == Adj_discrete) THEN ! Save stage solution - ! CALL WCOPY(NVAR,Y,1,Ystage(1),1) Ystage(1:NVAR) = Y(1:NVAR) - CALL WCOPY(NVAR,Y,1,Ynew,1) + Ynew(1:NVAR) = Y(1:NVAR) END IF - ! istage>1 and a new function evaluation is needed at the current istage - ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR,Y,1,Ynew,1) + ! istage>1 and a new function evaluation is needed at the current istage + ELSEIF ( ros_NewF(istage) ) THEN + Ynew(1:NVAR) = Y(1:NVAR) DO j = 1, istage-1 - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1),1,Ynew,1) + Ynew(1:NVAR) = Ynew(1:NVAR) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - END IF ! if istage == 1 elseif ros_NewF(istage) + END IF ! if istage == 1 elseif ros_NewF(istage) ! save stage solution every time even if ynew is not updated - IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN - ! CALL WCOPY(NVAR,Ynew,1,Ystage(ioffset+1),1) + IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR) - END IF - CALL WCOPY(NVAR,Fcn,1,K(ioffset+1),1) - DO j = 1, istage-1 + END IF + K(ioffset+1:ioffset+NVAR) = Fcn(1:NVAR) + DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1) - END DO - IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN + K(ioffset+1:ioffset+NVAR) = K(ioffset+1:ioffset+NVAR) & + + HC * K(NVAR*(j-1)+1:NVAR*j) + END DO + IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1) - END IF - CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1)) + K(ioffset+1:ioffset+NVAR) = K(ioffset+1:ioffset+NVAR) + HG * dFdT(1:NVAR) + END IF + CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1)) END DO Stage - !~~~> Compute the new solution - CALL WCOPY(NVAR,Y,1,Ynew,1) + Ynew(1:NVAR) = Y(1:NVAR) DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1) + Ynew(1:NVAR) = Ynew(1:NVAR) + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j) END DO !~~~> Compute the error estimation - CALL WSCAL(NVAR,ZERO,Yerr,1) + Yerr(1:NVAR) = ZERO DO j=1,ros_S - CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1) + Yerr(1:NVAR) = Yerr(1:NVAR) + ros_E(j) * K(NVAR*(j-1)+1:NVAR*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -1079,11 +1067,11 @@ SUBROUTINE ros_FwdInt ( Y, & CALL Jac_SP_Vec( Jac0, Fcn0, K(1) ) #endif IF (.NOT. Autonomous) THEN - CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) + K(1:NVAR) = K(1:NVAR) + dFdT(1:NVAR) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) END IF - CALL WCOPY(NVAR,Ynew,1,Y,1) + Y(1:NVAR) = Ynew(1:NVAR) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1126,7 +1114,7 @@ SUBROUTINE ros_FwdInt ( Y, & #endif IF (.NOT. Autonomous) THEN CALL ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) - CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) + K(1:NVAR) = K(1:NVAR) + dFdT(1:NVAR) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) !~~~> Deallocate stage buffer: only needed for discrete adjoint @@ -1212,7 +1200,7 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & Ghimj(i,i) = Ghimj(i,i)+Tau END DO #else - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau END DO @@ -1231,16 +1219,18 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & !~~~> Compute U DO m = 1,NADJ - CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1) - CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1) + U(istart:istart+NVAR-1,m) = Lambda(1:NVAR,m) + U(istart:istart+NVAR-1,m) = ros_M(istage) * U(istart:istart+NVAR-1,m) END DO ! m=1:NADJ DO j = istage+1, ros_S jstart = NVAR*(j-1) + 1 HA = ros_A((j-1)*(j-2)/2+istage) HC = ros_C((j-1)*(j-2)/2+istage)/(Direction*H) DO m = 1,NADJ - CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1) - CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1) + U(istart:istart+NVAR-1,m) = U(istart:istart+NVAR-1,m) & + + HA * V(jstart:jstart+NVAR-1,m) + U(istart:istart+NVAR-1,m) = U(istart:istart+NVAR-1,m) & + + HC * U(jstart:jstart+NVAR-1,m) END DO ! m=1:NADJ END DO DO m = 1,NADJ @@ -1273,10 +1263,10 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & istart = NVAR*(istage-1) + 1 DO m = 1,NADJ ! Add V_i - CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + V(istart:istart+NVAR-1,m) ! Add (H0xK_i)^T * U_i CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp ) - CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + Tmp(1:NVAR) END DO ! m=1:NADJ END DO ! Add H * dJac_dT_0^T * \sum(gamma_i U_i) @@ -1286,14 +1276,15 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & Tmp(1:NVAR) = ZERO DO istage = 1, ros_S istart = NVAR*(istage-1) + 1 - CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1) + Tmp(1:NVAR) = Tmp(1:NVAR) & + + ros_Gamma(istage) * U(istart:istart+NVAR-1,m) END DO #ifdef FULL_ALGEBRA Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp) #else CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2) #endif - CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + H * Tmp2(1:NVAR) END DO ! m=1:NADJ END IF ! .NOT.Autonomous @@ -1404,7 +1395,7 @@ SUBROUTINE ros_CadjInt ( & #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif - CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) + dFdT(1:NVAR,iadj) = -dFdT(1:NVAR,iadj) END DO END IF @@ -1412,7 +1403,7 @@ SUBROUTINE ros_CadjInt ( & #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) + Jac0(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1441,16 +1432,17 @@ SUBROUTINE ros_CadjInt ( & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) + Fcn(1:NVAR,iadj) = Fcn0(1:NVAR,iadj) END DO ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) + Ynew(1:NVAR,1:NADJ) = Y(1:NVAR,1:NADJ) DO j = 1, istage-1 - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) - END DO + DO iadj = 1, NADJ + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO Tau = T + ros_Alpha(istage)*Direction*H ! CALL FunTemplate(Tau,Ynew,Fcn) @@ -1461,7 +1453,7 @@ SUBROUTINE ros_CadjInt ( & #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) + Jac(1:LU_NONZERO) = -Jac(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1469,25 +1461,25 @@ SUBROUTINE ros_CadjInt ( & #else CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj)) #endif - !CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1) END DO END IF ! if istage == 1 elseif ros_NewF(istage) DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) + K(ioffset+1:ioffset+NVAR,iadj) = Fcn(1:NVAR,iadj) END DO DO j = 1, istage-1 - HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & - K(ioffset+1,iadj),1) - END DO + HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HC * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN - HG = Direction*H*ros_Gamma(istage) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) - END DO + HG = Direction*H*ros_Gamma(istage) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HG * dFdT(1:NVAR,iadj) + END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) @@ -1498,18 +1490,20 @@ SUBROUTINE ros_CadjInt ( & !~~~> Compute the new solution DO iadj = 1, NADJ - CALL WCOPY(NVAR,Y(1,iadj),1,Ynew(1,iadj),1) + Ynew(1:NVAR,iadj) = Y(1:NVAR,iadj) DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) END DO END DO !~~~> Compute the error estimation - CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1) + Yerr(1:NVAR,1:NADJ) = ZERO DO j=1,ros_S - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1,iadj),1,Yerr(1,iadj),1) - END DO + DO iadj = 1, NADJ + Yerr(1:NVAR,iadj) = Yerr(1:NVAR,iadj) & + + ros_E(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO !~~~> Max error among all adjoint components iadj = 1 @@ -1524,7 +1518,7 @@ SUBROUTINE ros_CadjInt ( & ! ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1) + Y(1:NVAR,1:NADJ) = Ynew(1:NVAR,1:NADJ) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1598,9 +1592,6 @@ SUBROUTINE ros_SimpleCadjInt ( & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1617,7 +1608,6 @@ SUBROUTINE ros_SimpleCadjInt ( & T = chk_T(istack) H = chk_H(istack-1) - !CALL WCOPY(NVAR,chk_Y(1,istack),1,Y0,1) Y0(1:NVAR) = chk_Y(1:NVAR,istack) !~~~> Compute the Jacobian at current time @@ -1633,7 +1623,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif - CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) + dFdT(1:NVAR,iadj) = -dFdT(1:NVAR,iadj) END DO END IF @@ -1641,7 +1631,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) + Jac0(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1659,8 +1649,7 @@ SUBROUTINE ros_SimpleCadjInt ( & Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv END DO @@ -1681,17 +1670,16 @@ SUBROUTINE ros_SimpleCadjInt ( & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) - END DO + Fcn(1:NVAR,1:NADJ) = Fcn0(1:NVAR,1:NADJ) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) + Ynew(1:NVAR,1:NADJ) = Y(1:NVAR,1:NADJ) DO j = 1, istage-1 - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) - END DO + DO iadj = 1, NADJ + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO Tau = T + ros_Alpha(istage)*Direction*H CALL ros_Hermite3( chk_T(istack-1), chk_T(istack), Tau, & @@ -1702,7 +1690,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) + Jac(1:LU_NONZERO) = -Jac(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1712,22 +1700,20 @@ SUBROUTINE ros_SimpleCadjInt ( & #endif END DO END IF ! if istage == 1 elseif ros_NewF(istage) - - DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) - END DO + K(ioffset+1:ioffset+NVAR,1:NADJ) = Fcn(1:NVAR,1:NADJ) DO j = 1, istage-1 - HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & - K(ioffset+1,iadj),1) - END DO + HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HC * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN - HG = Direction*H*ros_Gamma(istage) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) - END DO + HG = Direction*H*ros_Gamma(istage) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HG * dFdT(1:NVAR,iadj) + END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) @@ -1739,7 +1725,8 @@ SUBROUTINE ros_SimpleCadjInt ( & !~~~> Compute the new solution DO iadj = 1, NADJ DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Y(1,iadj),1) + Y(1:NVAR,iadj) = Y(1:NVAR,iadj) & + + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) END DO END DO @@ -1801,8 +1788,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(NVAR,(ONE/Delta),dFdT,1) + dFdT(1:NVAR) = dFdT(1:NVAR) - Fcn0(1:NVAR) + dFdT(1:NVAR) = dFdT(1:NVAR) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1830,11 +1817,11 @@ SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, Jac0, dJdT ) CALL JacTemplate( T+Delta, Y, dJdT ) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA - CALL WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1) + dJdT(1:NVAR,1:NVAR) = dJdT(1:NVAR,1:NVAR) - Jac0(1:NVAR,1:NVAR) + dJdT(1:NVAR,1:NVAR) = dJdT(1:NVAR,1:NVAR) * (ONE/Delta) #else - CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) - Jac0(1:LU_NONZERO) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) * (ONE/Delta) #endif END SUBROUTINE ros_JacTimeDerivative @@ -1884,15 +1871,15 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - CALL WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1) - CALL WSCAL(NVAR*NVAR,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv @@ -2058,30 +2045,30 @@ SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y ) amb(i) = amb(i-1)*amb(1) END DO - ! c(1) = ya; - CALL WCOPY(NVAR,Ya,1,C(1,1),1) + C(1:NVAR,1) = Ya(1:NVAR) + ! c(2) = ja; - CALL WCOPY(NVAR,Ja,1,C(1,2),1) + C(1:NVAR,2) = Ja(1:NVAR) + ! c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb ; - CALL WCOPY(NVAR,Ya,1,C(1,3),1) - CALL WSCAL(NVAR,-3.0*amb(2),C(1,3),1) - CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1,3),1) - CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1,3),1) - CALL WAXPY(NVAR,amb(1),Jb,1,C(1,3),1) + C(1:NVAR,3) = 2.0*amb(1) * Ja(1:NVAR) & + + amb(1) * Jb(1:NVAR) & + - 3.0*amb(2) * Ya(1:NVAR) & + + 3.0*amb(2) * Yb(1:NVAR) + ! c(4) = 1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb ; - CALL WCOPY(NVAR,Ya,1,C(1,4),1) - CALL WSCAL(NVAR,-2.0*amb(3),C(1,4),1) - CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1,4),1) - CALL WAXPY(NVAR,amb(2),Ja,1,C(1,4),1) - CALL WAXPY(NVAR,amb(2),Jb,1,C(1,4),1) + C(1:NVAR,4) = amb(2) * Ja(1:NVAR) & + + amb(2) * Jb(1:NVAR) & + - 2.0*amb(3) * Ya(1:NVAR) & + + 2.0*amb(3) * Yb(1:NVAR) +! Unrolled loop: Y = Tau^3*c(4) + Tau^2*c(3) + Tau*c(2) + c(1) Tau = T - a - CALL WCOPY(NVAR,C(1,4),1,Y,1) - CALL WSCAL(NVAR,Tau**3,Y,1) - DO j = 3,1,-1 - CALL WAXPY(NVAR,TAU**(j-1),C(1,j),1,Y,1) - END DO + Y(1:NVAR) = Tau**3 * C(1:NVAR,4) & + + Tau**2 * C(1:NVAR,3) & + + Tau * C(1:NVAR,2) & + + C(1:NVAR,1) END SUBROUTINE ros_Hermite3 @@ -2109,45 +2096,43 @@ SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y ) END DO ! c(1) = ya; - CALL WCOPY(NVAR,Ya,1,C(1,1),1) + C(1:NVAR,1) = Ya(1:NVAR) + ! c(2) = ja; - CALL WCOPY(NVAR,Ja,1,C(1,2),1) + C(1:NVAR,2) = Ja(1:NVAR) + ! c(3) = ha/2; - CALL WCOPY(NVAR,Ha,1,C(1,3),1) - CALL WSCAL(NVAR,HALF,C(1,3),1) + C(1:NVAR,3) = Ha(1:NVAR) * HALF ! c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb + 1.5*amb(1)*ha - 0.5*amb(1)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,4),1) - CALL WSCAL(NVAR,10.0*amb(3),C(1,4),1) - CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1,4),1) - CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1,4),1) - CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1,4),1) - CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1,4),1) - CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1,4),1) + C(1:NVAR,4) = 10.0*amb(3) * Ya(1:NVAR) & + - 10.0*amb(3) * Yb(1:NVAR) & + - 6.0*amb(2) * Ja(1:NVAR) & + - 4.0*amb(2) * Jb(1:NVAR) & + + 1.5*amb(1) * Ha(1:NVAR) & + - 0.5*amb(1) * Hb(1:NVAR) ! c(5) = 15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb + 1.5*amb(2)*ha - 1*amb(2)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,5),1) - CALL WSCAL(NVAR, 15.0*amb(4),C(1,5),1) - CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1,5),1) - CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1,5),1) - CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1,5),1) - CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1,5),1) - CALL WAXPY(NVAR,-amb(2),Hb,1,C(1,5),1) + C(1:NVAR,5) = 15.0*amb(4) * Ya(1:NVAR) & + - 15.0*amb(4) * Yb(1:NVAR) & + - 8.0*amb(3) * Ja(1:NVAR) & + - 7.0*amb(3) * Jb(1:NVAR) & + + 1.5*amb(2) * Ha(1:NVAR) & + - amb(2) * Hb(1:NVAR) ! c(6) = 6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb + 0.5*amb(3)*ha -0.5*amb(3)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,6),1) - CALL WSCAL(NVAR, 6.0*amb(5),C(1,6),1) - CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1,6),1) - CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1,6),1) - CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1,6),1) - CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1,6),1) - CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1,6),1) + C(1:NVAR,6) = 6.0*amb(5) * Ya(1:NVAR) & + - 6.0*amb(5) * Yb(1:NVAR) & + - 3.0*amb(4) * Ja(1:NVAR) & + - 3.0*amb(4) * Jb(1:NVAR) & + + 0.5*amb(3) * Ha(1:NVAR) & + - 0.5*amb(3) * Hb(1:NVAR) + Tau = T - a - CALL WCOPY(NVAR,C(1,6),1,Y,1) + Y(1:NVAR) = C(1:NVAR,6) DO j = 5,1,-1 - CALL WSCAL(NVAR,Tau,Y,1) - CALL WAXPY(NVAR,ONE,C(1,j),1,Y,1) + Y(1:NVAR) = Tau * Y(1:NVAR) + C(1:NVAR,j) END DO END SUBROUTINE ros_Hermite5 diff --git a/int/rosenbrock_autoreduce.f90 b/int/rosenbrock_autoreduce.f90 index a7a8e07..b6e2821 100644 --- a/int/rosenbrock_autoreduce.f90 +++ b/int/rosenbrock_autoreduce.f90 @@ -307,7 +307,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ISTATUS(1:8) = 0 RSTATUS(1:4) = ZERO -!~~~> Autonomous or time dependent ODE. Default is time dependent. +!~~~> Autonomous (1) or time dependent ODE (0). Default is time dependent. Autonomous = .NOT.(ICNTRL(1) == 0) !~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1) @@ -360,7 +360,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & AR_target_spc = ICNTRL(14) !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -528,7 +528,7 @@ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR) CASE (-6) PRINT * , '--> No of steps exceeds maximum bound' CASE (-7) - PRINT * , '--> Step size too small: T + 10*H = T', & + PRINT * , '--> Step size too small: T + 0.1*H = T', & ' or H < Roundoff' CASE (-8) PRINT * , '--> Matrix is repeatedly singular' @@ -587,9 +587,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -661,29 +658,28 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + & + ros_A((istage-1)*(istage-2)/2+j) * K(N*(j-1)+1:N*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = Fcn(1:N) DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HC * K(N*(j-1)+1:N*j) END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HG * dFdT(1:N) END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) @@ -691,17 +687,15 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -717,7 +711,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -1018,13 +1011,11 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - call WCOPY(N,Fcn0,1,Fcn,1) - ! Fcn(1:N) = Fcn0(1:N) + Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ! K = 0.0_dp ! is this fix needed? hplin 14:04 -- not. 3 hours wiser later ELSEIF ( ros_NewF(istage) ) THEN - call WCOPY(N,Y,1,Ynew,1) - ! Ynew(1:N) = Y(1:N) + Ynew(1:N) = Y(1:N) DO j = 1, istage-1 ! In full vector space. Just use WAXPY as normal ! other entries in K are set to 1 previously. @@ -1090,7 +1081,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HC * K(N*(j-1)+SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HC * K(N*(j-1)+SPC_MAP(i)) ENDDO ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP) ! loop unrolling is consistently slower here. 18:58 @@ -1102,7 +1094,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HG * dFdT(SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HG * dFdT(SPC_MAP(i)) ENDDO ENDIF @@ -1153,8 +1146,7 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - !Y(1:N) = Ynew(1:N) + Y(1:N) = Ynew(1:N) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1482,13 +1474,11 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - call WCOPY(N,Fcn0,1,Fcn,1) - ! Fcn(1:N) = Fcn0(1:N) + Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ! K = 0.0_dp ! is this fix needed? hplin 14:04 -- not. 3 hours wiser later ELSEIF ( ros_NewF(istage) ) THEN - call WCOPY(N,Y,1,Ynew,1) - ! Ynew(1:N) = Y(1:N) + Ynew(1:N) = Y(1:N) DO j = 1, istage-1 ! In full vector space. Just use WAXPY as normal ! other entries in K are set to 1 previously. @@ -1554,7 +1544,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HC * K(N*(j-1)+SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HC * K(N*(j-1)+SPC_MAP(i)) ENDDO ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP) ! loop unrolling is consistently slower here. 18:58 @@ -1566,7 +1557,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HG * dFdT(SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HG * dFdT(SPC_MAP(i)) ENDDO ENDIF @@ -1616,8 +1608,7 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - !Y(1:N) = Ynew(1:N) + Y(1:N) = Ynew(1:N) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1729,8 +1720,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1778,16 +1769,12 @@ SUBROUTINE ros_cPrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,rNVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:cNONZERO) = -Jac0(JVS_MAP(1:cNONZERO)) ghinv = ONE/(Direction*H*gam) DO i=1,rNVAR @@ -1873,8 +1860,6 @@ SUBROUTINE ros_cSolve( A, Pivot, b, map1, map2 ) Btmp = 0.d0 Atmp(map1(1:cNONZERO)) = A btmp(map2(1:rNVAR)) = b -! call cWCOPY(cNONZERO,LU_NONZERO,A,1,Atmp,1,map1) -! call cWCOPY(rNVAR,NVAR,B,1,Btmp,1,map2) CALL KppSolve( Atmp, btmp ) b = btmp(map2(1:rNVAR)) #endif @@ -1926,16 +1911,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_h211b_qssa.f90 b/int/rosenbrock_h211b_qssa.f90 index e925687..2d209ca 100644 --- a/int/rosenbrock_h211b_qssa.f90 +++ b/int/rosenbrock_h211b_qssa.f90 @@ -409,7 +409,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ENDIF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -662,9 +662,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -837,31 +834,25 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - ! CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - ! K(N*(j-1)+1),1,Ynew,1) - Ynew(1:N) = Ynew(1:N) + K(N*(j-1)+1:N*j) * ros_A((istage-1)*(istage-2)/2+j) + Ynew(1:N) = Ynew(1:N) + & + K(N*(j-1)+1:N*j) * ros_A((istage-1)*(istage-2)/2+j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = Fcn(1:N) DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - ! CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + K(N*(j-1)+1:N*j) * HC END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) -! CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + dFdT(1:N) * HG END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) @@ -870,18 +861,14 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - ! CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) Ynew(1:N) = Ynew(1:N) + K(N*(j-1)+1:N*j) * ros_m(j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - ! CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) Yerr(1:N) = Yerr(1:N) + K(N*(j-1)+1:N*j) * ros_E(j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -908,7 +895,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -995,9 +981,7 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - ! CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) dFdt(1:N) = dFdt(1:N) - ONE * FcN0(1:N) - ! CALL WSCAL(N,(ONE/Delta),dFdT,1) dFDT(1:N) = dFDT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1045,16 +1029,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_tlm.f90 b/int/rosenbrock_tlm.f90 index 4710de3..de5f3d4 100644 --- a/int/rosenbrock_tlm.f90 +++ b/int/rosenbrock_tlm.f90 @@ -372,7 +372,7 @@ SUBROUTINE RosenbrockTLM(N,Y,NTLM,Y_tlm, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -545,9 +545,6 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & LOGICAL :: RejectLastH, RejectMoreH, Singular !~~~> Local parameters KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -625,21 +622,23 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & ioffset = N*(istage-1) ! Initialize stage solution - CALL WCOPY(N,Y,1,Ynew,1) - CALL WCOPY(N*NTLM,Y_tlm,1,Ynew_tlm,1) + Ynew(1:N) = Y(1:N) + Ynew_tlm(1:N,1:NTLM) = Y_tlm(1:N,1:NTLM) ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - CALL WCOPY(N,Fcn0,1,Fcn,1) - CALL WCOPY(N*NTLM,Fcn0_tlm,1,Fcn_tlm,1) + Fcn(1:N) = Fcn0(1:N) + Fcn_tlm(1:N,1:NTLM) = Fcn0_tlm(1:N,1:NTLM) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(N*(j-1)+1:N*j) DO itlm=1,NTLM - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K_tlm(N*(j-1)+1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Ynew_tlm(1:N,itlm) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO Tau = T + ros_Alpha(istage)*Direction*H @@ -651,29 +650,33 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & CALL Jac_SP_Vec ( Jac, Ynew_tlm(1,itlm), Fcn_tlm(1,itlm) ) END DO END IF ! if istage == 1 elseif ros_NewF(istage) - CALL WCOPY(N,Fcn,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = Fcn(1:N) DO itlm=1,NTLM - CALL WCOPY(N,Fcn_tlm(1,itlm),1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = Fcn_tlm(1:N,itlm) END DO DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HC * K(N*(j-1)+1:N*j) DO itlm=1,NTLM - CALL WAXPY(N,HC,K_tlm(N*(j-1)+1,itlm),1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + HC * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HG * dFdT(1:N) DO itlm=1,NTLM CALL Jac_SP_Vec ( dJdT, Ynew_tlm(1,itlm), Tmp ) - CALL WAXPY(N,HG,Tmp,1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + HG * Tmp(1:N) END DO END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) DO itlm=1,NTLM CALL Hess_Vec ( Hes0, K(ioffset+1), Y_tlm(1,itlm), Tmp ) - CALL WAXPY(N,ONE,Tmp,1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + Tmp(1:N) CALL ros_Solve(Ghimj, Pivot, K_tlm(ioffset+1,itlm)) END DO @@ -681,29 +684,31 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & !~~~> Compute the new solution - CALL WCOPY(N,Y,1,Ynew,1) + Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO DO itlm=1,NTLM - CALL WCOPY(N,Y_tlm(1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Y_tlm(1:N,itlm) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K_tlm(N*(j-1)+1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Ynew_tlm(1:N,itlm) & + + ros_M(j) * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO !~~~> Compute the error estimation - CALL Set2zero(N,Yerr) + Yerr(1:N) = 0.0_dp DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) IF (TLMtruncErr) THEN Err1 = 0.0d0 - CALL Set2zero(N*NTLM,Yerr_tlm) + Yerr_tlm(1:N,1:NTLM) = 0.0_dp DO itlm=1,NTLM DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K_tlm(N*(j-1)+1,itlm),1,Yerr_tlm(1,itlm),1) + Yerr_tlm(1:N,itlm) = Yerr_tlm(1:N,itlm) & + + ros_E(j) * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO Err = ros_ErrorNorm_tlm(Y_tlm,Ynew_tlm,Yerr_tlm,AbsTol_tlm,RelTol_tlm,Err,VectorTol) @@ -717,8 +722,8 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - CALL WCOPY(N*NTLM,Ynew_tlm,1,Y_tlm,1) + Y(1:N) = Ynew(1:N) + Y_tlm(1:N,1:NTLM) = Ynew_tlm(1:N,1:NTLM) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -832,8 +837,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -856,8 +861,8 @@ SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, Jac0, dJdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL JacTemplate( T+Delta, Y, dJdT ) ISTATUS(Njac) = ISTATUS(Njac) + 1 - CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) - Jac0(1:LU_NONZERO) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) * (ONE/Delta) END SUBROUTINE ros_JacTimeDerivative @@ -895,8 +900,8 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & DO WHILE (Singular) !~~~> Construct Ghimj = 1/(H*ham) - Jac0 - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv diff --git a/int/runge_kutta.f90 b/int/runge_kutta.f90 index 27bbab7..73f2eae 100644 --- a/int/runge_kutta.f90 +++ b/int/runge_kutta.f90 @@ -19,8 +19,9 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Precision USE KPP_ROOT_Parameters USE KPP_ROOT_Global - USE KPP_ROOT_Jacobian, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra + USE KPP_ROOT_Jacobian, ONLY : LU_DIAG + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppDecompCmplx, KppSolveCmplx IMPLICIT NONE PUBLIC @@ -391,7 +392,7 @@ SUBROUTINE RungeKutta( N,T,Tend,Y,RelTol,AbsTol, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -591,9 +592,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -640,9 +641,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -670,11 +671,11 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Prepare the loop-independent part of the right-hand side ! G = H*rkBgam(0)*F0 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - IF (rkMethod/=L3A) CALL WAXPY(N,rkBgam(0)*H, F0,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + IF (rkMethod/=L3A) G(1:N) = G(1:N) + rkBgam(0)*H * F0(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -683,12 +684,12 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -723,8 +724,8 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) - + Z4(1:N) = Z4(1:N) + DZ4(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT SDNewtonLoop @@ -742,21 +743,21 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN - CALL Set2Zero(N, DZ4) + DZ4(1:N) = 0.0_dp IF (rkMethod==L3A) THEN DZ4(1:N) = H*rkF(0)*F0(1:N) - IF (rkF(1) /= ZERO) CALL WAXPY(N, rkF(1), Z1, 1, DZ4, 1) - IF (rkF(2) /= ZERO) CALL WAXPY(N, rkF(2), Z2, 1, DZ4, 1) - IF (rkF(3) /= ZERO) CALL WAXPY(N, rkF(3), Z3, 1, DZ4, 1) + IF (rkF(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(1) * Z1(1:N) + IF (rkF(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(2) * Z2(1:N) + IF (rkF(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(3) * Z3(1:N) TMP = Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ1) - CALL WAXPY(N, H*rkBgam(4), DZ1, 1, DZ4, 1) + DZ4(1:N) = DZ4(1:N) + H*rkBgam(4) * DZ1(1:N) ELSE ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) END IF Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE @@ -790,9 +791,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum(d_i Z_i) - IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1) /= ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) @@ -978,33 +979,33 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,F0,Z1,Z2,Z3,R1,R2,R3) KPP_REAL :: T, H KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F0,F,R1,R2,R3,TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 IF (rkMethod==L3A) THEN - CALL WAXPY(N,-H*rkA(1,0),F0,1,R1,1) ! R1 <- R1 - h*A_10*F0 - CALL WAXPY(N,-H*rkA(2,0),F0,1,R2,1) ! R2 <- R2 - h*A_20*F0 - CALL WAXPY(N,-H*rkA(3,0),F0,1,R3,1) ! R3 <- R3 - h*A_30*F0 + R1(1:N) = R1(1:N) - H*rkA(1,0) * F0(1:N) ! R1 <- R1 - h*A_10*F0 + R2(1:N) = R2(1:N) - H*rkA(2,0) * F0(1:N) ! R2 <- R2 - h*A_20*F0 + R3(1:N) = R3(1:N) - H*rkA(3,0) * F0(1:N) ! R3 <- R3 - h*A_30*F0 END IF - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS diff --git a/int/runge_kutta_adj.f90 b/int/runge_kutta_adj.f90 index d7008ce..1074f8c 100644 --- a/int/runge_kutta_adj.f90 +++ b/int/runge_kutta_adj.f90 @@ -20,8 +20,10 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters USE KPP_ROOT_Global USE KPP_ROOT_Jacobian - USE KPP_ROOT_LinearAlgebra - + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppSolveTR, KppSolveCmplx, & + KppSolveTRCmplx, KppDecompCmplx, & + WGEFA, WGESL IMPLICIT NONE PUBLIC SAVE @@ -476,7 +478,7 @@ SUBROUTINE RungeKuttaADJ( N, Y, NADJ, Lambda,Tstart,Tend, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff = WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -810,15 +812,11 @@ SUBROUTINE rk_DPush( T, H, Y, Zstage, NewIt, E1, E2 )!, Jcb ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - ! CALL WCOPY(NVAR*3,Zstage,1,chk_Z(1,stack_ptr),1) chk_Y(1:N,stack_ptr) = Y(1:N) chk_Z(1:3*N,stack_ptr) = Zstage(1:3*N) chk_NiT( stack_ptr ) = NewIt IF (SaveLU) THEN #ifdef FULL_ALGEBRA - ! CALL WCOPY(NVAR*NVAR, E1,1,chk_E1(1,stack_ptr),1) - ! CALL WCOPYCmplx(NVAR*NVAR, E2,1,chk_E2(1,stack_ptr),1) DO j=1,NVAR DO i=1,NVAR chk_E1(NVAR*(j-1)+i,stack_ptr) = E1(i,j) @@ -826,13 +824,10 @@ SUBROUTINE rk_DPush( T, H, Y, Zstage, NewIt, E1, E2 )!, Jcb ) END DO END DO #else - ! CALL WCOPY(LU_NONZERO, E1,1,chk_E1(1,stack_ptr),1) - ! CALL WCOPYCmplx(LU_NONZERO, E2,1,chk_E2(1,stack_ptr),1) chk_E1(1:LU_NONZERO,stack_ptr) = E1(1:LU_NONZERO) chk_E2(1:LU_NONZERO,stack_ptr) = E2(1:LU_NONZERO) #endif END IF - !CALL WCOPY(LU_NONZERO,Jcb,1,chk_J(1,stack_ptr),1) END SUBROUTINE rk_DPush @@ -860,15 +855,11 @@ SUBROUTINE rk_DPop( T, H, Y, Zstage, NewIt, E1, E2 ) !, Jcb ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) - ! CALL WCOPY(NVAR*3,chk_Z(1,stack_ptr),1,Zstage,1) Zstage(1:3*NVAR) = chk_Z(1:3*NVAR,stack_ptr) NewIt = chk_NiT( stack_ptr ) IF (SaveLU) THEN #ifdef FULL_ALGEBRA - ! CALL WCOPY(NVAR*NVAR,chk_E1(1,stack_ptr),1, E1,1) - ! CALL WCOPYCmplx(NVAR*NVAR,chk_E2(1,stack_ptr),1, E2,1) DO j=1,NVAR DO i=1,NVAR E1(i,j) = chk_E1(NVAR*(j-1)+i,stack_ptr) @@ -876,13 +867,10 @@ SUBROUTINE rk_DPop( T, H, Y, Zstage, NewIt, E1, E2 ) !, Jcb ) END DO END DO #else - ! CALL WCOPY(LU_NONZERO,chk_E1(1,stack_ptr),1, E1,1) - ! CALL WCOPYCmplx(LU_NONZERO,chk_E2(1,stack_ptr),1, E2,1) E1(1:LU_NONZERO) = chk_E1(1:LU_NONZERO,stack_ptr) E2(1:LU_NONZERO) = chk_E2(1:LU_NONZERO,stack_ptr) #endif END IF - !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1) stack_ptr = stack_ptr - 1 @@ -903,9 +891,6 @@ SUBROUTINE rk_CPush(T, H, Y, dY, d2Y ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - ! CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1) - ! CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1) chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR) chk_d2Y(1:NVAR,stack_ptr)= d2Y(1:NVAR) @@ -927,9 +912,6 @@ SUBROUTINE rk_CPop( T, H, Y, dY, d2Y ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) - ! CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1) - ! CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr) d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr) @@ -1034,9 +1016,9 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -1084,9 +1066,9 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -1122,11 +1104,11 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! G = H*rkBgam(0)*DZ4 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - CALL WAXPY(N,rkBgam(0)*H, DZ4,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + G(1:N) = G(1:N) + rkBgam(0)*H * DZ4(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -1135,12 +1117,13 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) + !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -1175,8 +1158,8 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) - + Z4(1:N) = Z4(1:N) + DZ4(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT SDNewtonLoop @@ -1196,12 +1179,12 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - CALL Set2Zero(N, DZ4) - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) - Err = RK_ErrorNorm(N,SCAL,DZ4) + DZ4(1:N) = 0.0_dp + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) + Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE CALL RK_ErrorEstimate(N,H,Y,T, & E1,IP1,Z1,Z2,Z3,SCAL,Err,FirstStep,Reject) @@ -1217,9 +1200,6 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ accept:IF (Err < ONE) THEN !~~~> STEP IS ACCEPTED IF (AdjointType == KPP_ROOT_discrete) THEN ! Save stage solution - ! CALL WCOPY(N,Z1,1,Zstage(1),1) - ! CALL WCOPY(N,Z2,1,Zstage(N+1),1) - ! CALL WCOPY(N,Z3,1,Zstage(2*N+1),1) Zstage(1:N) = Z1(1:N) Zstage(N+1:2*N) = Z2(1:N) Zstage(2*N+1:3*N) = Z3(1:N) @@ -1243,9 +1223,10 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum (d_i Z_i) - IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1) /= ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) + ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) @@ -1367,12 +1348,12 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF !~~~> Jacobian values at stage vectors - CALL WADD(N,Y,Zstage(1),TMP) ! TMP <- Y + Z1 - CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) - CALL WADD(N,Y,Zstage(1+N),TMP) ! TMP <- Y + Z2 - CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) - CALL WADD(N,Y,Zstage(1+2*N),TMP) ! TMP <- Y + Z3 - CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) + TMP(1:N) = Y(1:N) + Zstage(1:N) ! TMP <- Y + Z1 + CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) + TMP(1:N) = Y(1:N) + Zstage(N+1:2*N) ! TMP <- Y + Z2 + CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) + TMP(1:N) = Y(1:N) + Zstage(2*N+1:3*N) ! TMP <- Y + Z3 + CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) END IF ! .not.Reject @@ -1398,7 +1379,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) Jbig(i,i) = Jbig(i,i) + ONE END DO CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,ISING) - ! CALL WGEFA(3*N,Jbig,IPbig,ISING) IF (ISING /= 0) THEN PRINT*,'Big guy is singular'; STOP END IF @@ -1435,7 +1415,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) DO i=1, 3*N Jbig(i,i) = ONE + Jbig(i,i) END DO - ! CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,ISING) CALL WGEFA(3*N,Jbig,IPbig,ISING) IF (ISING /= 0) THEN PRINT*,'Big guy is singular'; STOP @@ -1448,13 +1427,10 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Adj:DO iadj = 1, NADJ !~~~> Starting values for Newton iteration - ! CALL WCOPY(N,Lambda(1,iadj),1,U1(1,iadj),1) - ! CALL WCOPY(N,Lambda(1,iadj),1,U2(1,iadj),1) - ! CALL WCOPY(N,Lambda(1,iadj),1,U3(1,iadj),1) - CALL Set2Zero(N,U1(1,iadj)) - CALL Set2Zero(N,U2(1,iadj)) - CALL Set2Zero(N,U3(1,iadj)) - + U1(1:N,iadj) = 0.0_dp + U2(1:N,iadj) = 0.0_dp + U3(1:N,iadj) = 0.0_dp + !~~~> Initializations for Newton iteration NewtonDone = .FALSE. !~~~> Right Hand Side - part G for all Newton iterations @@ -1508,9 +1484,9 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF ! (AdjointSolve == Solve_adaptive) ! Update solution - CALL WAXPY(N,-ONE,DU1,1,U1(1,iadj),1) ! U1 <- U1 - DU1 - CALL WAXPY(N,-ONE,DU2,1,U2(1,iadj),1) ! U2 <- U2 - DU2 - CALL WAXPY(N,-ONE,DU3,1,U3(1,iadj),1) ! U3 <- U3 - DU3 + U1(1:N,iadj) = U1(1:N,iadj) - DU1(1:N) ! U1 <- U1 - DU1 + U2(1:N,iadj) = U2(1:N,iadj) - DU2(1:N) ! U2 <- U2 - DU2 + U3(1:N,iadj) = U3(1:N,iadj) - DU3(1:N) ! U3 <- U3 - DU3 IF (AdjointSolve == Solve_adaptive) THEN ! When performing an adaptive number of iterations @@ -1531,9 +1507,9 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF ! Update adjoint solution: Y_adj <- Y_adj + sum (U_i) - CALL WAXPY(N,ONE,U1(1,iadj),1,Lambda(1,iadj),1) - CALL WAXPY(N,ONE,U2(1,iadj),1,Lambda(1,iadj),1) - CALL WAXPY(N,ONE,U3(1,iadj),1,Lambda(1,iadj),1) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U1(1:N,iadj) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U2(1:N,iadj) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U3(1:N,iadj) ELSE ! NewtonConverge = .false. @@ -1542,7 +1518,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) X(N+1:2*N) = -G2(1:N) X(2*N+1:3*N) = -G3(1:N) CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,ISING) - ! CALL WGESL('T',3*N,Jbig,IPbig,X) Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N) #else ! Commented lines for sparse big algebra: @@ -1555,7 +1530,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) X(1:N) = -G1(1:N) X(N+1:2*N) = -G2(1:N) X(2*N+1:3*N) = -G3(1:N) - ! CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,ISING) CALL WGESL('T',3*N,Jbig,IPbig,X) Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N) #endif @@ -1758,28 +1732,28 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,R1,R2,R3) KPP_REAL :: T,H KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F,R1,R2,R3,TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 - - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 - + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 + + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP, F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 + END SUBROUTINE RK_PrepareRHS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1802,46 +1776,48 @@ SUBROUTINE RK_PrepareRHS_Adj(N,H,Jac1,Jac2,Jac3,Lambda, & KPP_REAL, DIMENSION(N) :: F,TMP - CALL WCOPY(N,G1,1,R1,1) ! R1 <- G1 - CALL WCOPY(N,G2,1,R2,1) ! R2 <- G2 - CALL WCOPY(N,G3,1,R3,1) ! R3 <- G3 + R1(1:N) = G1(1:N) ! R1 <- G1 + R2(1:N) = G2(1:N) ! R2 <- G2 + R3(1:N) = G3(1:N) ! R3 <- G3 + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,1),U1,1,F,1) ! F1 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,1),U2,1,F,1) ! F1 <- F1 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,1),U3,1,F,1) ! F1 <- F1 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac1),F) #else - CALL JacTR_SP_Vec ( Jac1, F, TMP ) ! R1 <- -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + CALL JacTR_SP_Vec ( Jac1, F, TMP ) ! R1 <- -Jac(Y+Z1)^t*h*sum(A_j1*U_j) #endif - CALL WAXPY(N,ONE,U1,1,TMP,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) - CALL WAXPY(N,ONE,TMP,1,R1,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + TMP(1:N) = TMP(1:N) + U1(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R1(1:N) = R1(1:N) + TMP(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,2),U1,1,F,1) ! F2 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,2),U2,1,F,1) ! F2 <- F2 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,2),U3,1,F,1) ! F2 <- F2 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac2),F) #else - CALL JacTR_SP_Vec ( Jac2, F, TMP ) ! R2 <- -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + CALL JacTR_SP_Vec ( Jac2, F, TMP ) ! R2 <- -Jac(Y+Z2)^t*h*sum(A_j2*U_j) #endif - CALL WAXPY(N,ONE,U2,1,TMP,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) - CALL WAXPY(N,ONE,TMP,1,R2,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + TMP(1:N) = TMP(1:N) + U2(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R2(1:N) = R2(1:N) + TMP(1:N) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,3),U1,1,F,1) ! F3 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,3),U2,1,F,1) ! F3 <- F3 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,3),U3,1,F,1) ! F3 <- F3 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac3),F) #else CALL JacTR_SP_Vec ( Jac3, F, TMP ) ! R3 <- -Jac(Y+Z3)^t*h*sum(A_j3*U_j) #endif - CALL WAXPY(N,ONE,U3,1,TMP,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j) - CALL WAXPY(N,ONE,TMP,1,R3,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j) - + TMP(1:N) = TMP(1:N) + U3(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R3(1:N) = R3(1:N) + TMP(1:N) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) END SUBROUTINE RK_PrepareRHS_Adj @@ -1864,29 +1840,29 @@ SUBROUTINE RK_PrepareRHS_G(N,H,Jac1,Jac2,Jac3,Lambda, & #endif KPP_REAL, DIMENSION(N) :: F - CALL SET2ZERO(N,G1) - CALL SET2ZERO(N,G2) - CALL SET2ZERO(N,G3) + G1(1:N) = 0.0_dp + G2(1:N) = 0.0_dp + G3(1:N) = 0.0_dp #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac1),Lambda) #else CALL JacTR_SP_Vec ( Jac1, Lambda, F ) ! F1 <- Jac(Y+Z1)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(1),F,1,G1,1) ! R1 <- R1 - h*B_1*F1 + G1(1:N) = G1(1:N) - H*rkB(1) * F(1:N) ! R1 <- R1 - h*B_1*F1 #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac2),Lambda) #else CALL JacTR_SP_Vec ( Jac2, Lambda, F ) ! F2 <- Jac(Y+Z2)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(2),F,1,G2,1) ! R2 <- R2 - h*B_2*F2 + G2(1:N) = G2(1:N) - H*rkB(2) * F(1:N) ! R2 <- R2 - h*B_2*F2 #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac3),Lambda) #else CALL JacTR_SP_Vec ( Jac3, Lambda, F ) ! F3 <- Jac(Y+Z3)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(3),F,1,G3,1) ! R3 <- R3 - h*B_3*F3 + G3(1:N) = G3(1:N) - H*rkB(3) * F(1:N) ! R3 <- R3 - h*B_3*F3 END SUBROUTINE RK_PrepareRHS_G diff --git a/int/runge_kutta_tlm.f90 b/int/runge_kutta_tlm.f90 index d0afc1f..b175444 100644 --- a/int/runge_kutta_tlm.f90 +++ b/int/runge_kutta_tlm.f90 @@ -21,6 +21,8 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Jacobian USE KPP_ROOT_LinearAlgebra + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppDecompCmplx, KppSolveCmplx IMPLICIT NONE PUBLIC @@ -436,7 +438,7 @@ SUBROUTINE RungeKuttaTLM( N, NTLM, Y, Y_tlm, T, Tend,RelTol,AbsTol, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -648,9 +650,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -697,9 +699,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -734,11 +736,11 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! G = H*rkBgam(0)*DZ4 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - CALL WAXPY(N,rkBgam(0)*H, DZ4,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + G(1:N) = G(1:N) + rkBgam(0)*H * DZ4(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -747,12 +749,13 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) + !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -786,7 +789,7 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) + Z4(1:N) = Z4(1:N) + DZ4(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -807,11 +810,12 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - CALL Set2Zero(N, DZ4) - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) + DZ4(1:N) = 0.0_dp + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) + Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE CALL RK_ErrorEstimate(N,H,Y,T, & @@ -824,11 +828,11 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) IF (Err < ONE) THEN !~~~> Jacobian values - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) TLMDIR: IF (TLMDirect) THEN @@ -853,7 +857,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Jbig(i,i) = ONE + Jbig(i,i) END DO !~~~> Solve the big system - ! CALL DGETRF(3*NVAR,3*NVAR,Jbig,3*NVAR,IPbig,j) CALL WGEFA(3*N,Jbig,IPbig,info) IF (info /= 0) THEN PRINT*,'Big big guy is singular'; STOP @@ -865,7 +868,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Zbig(2*NVAR+j) = Y_tlm(j,itlm) END DO Zbig = MATMUL(Ebig,Zbig) - !CALL DGETRS ('N',3*NVAR,1,Jbig,3*NVAR,IPbig,Zbig,3*NVAR,ISING) CALL WGESL('N',3*N,Jbig,IPbig,Zbig) DO j=1,NVAR Z1_tlm(j,itlm) = Zbig(j) @@ -907,7 +909,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) DO i=1, 3*N Jbig(i,i) = ONE + Jbig(i,i) END DO - ! CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,info) CALL WGEFA(3*N,Jbig,IPbig,info) IF (info /= 0) THEN PRINT*,'Big guy is singular'; STOP @@ -927,7 +928,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) ! Compute RHS CALL RK_PrepareRHS_TLMdirect(N,H,Jac1,Jac2,Jac3,Y_tlm(1,itlm),Zbig) ! Solve the system - ! CALL DGETRS('N',3*N,1,Jbig,3*N,IPbig,Zbig,3*N,ISING) CALL WGESL('N',3*N,Jbig,IPbig,Zbig) Z1_tlm(1:NVAR,itlm) = Zbig(1:NVAR) Z2_tlm(1:NVAR,itlm) = Zbig(NVAR+1:2*NVAR) @@ -944,9 +944,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Tlm:DO itlm = 1, NTLM !~~~> Starting values for Newton iteration - CALL Set2zero(N,Z1_tlm(1,itlm)) - CALL Set2zero(N,Z2_tlm(1,itlm)) - CALL Set2zero(N,Z3_tlm(1,itlm)) + Z1_tlm(1:N,itlm) = 0.0_dp + Z2_tlm(1:N,itlm) = 0.0_dp + Z3_tlm(1:N,itlm) = 0.0_dp !~~~> Initializations for Newton iteration IF (TLMNewtonEst) THEN @@ -998,9 +998,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) END IF !(TLMNewtonEst) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1_tlm(1,itlm),1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2_tlm(1,itlm),1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3_tlm(1,itlm),1) ! Z3 <- Z3 - DZ3 + Z1_tlm(1:N,itlm) = Z1_tlm(1:N,itlm) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2_tlm(1:N,itlm) = Z2_tlm(1:N,itlm) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3_tlm(1:N,itlm) = Z3_tlm(1:N,itlm) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations IF (TLMNewtonEst) THEN @@ -1061,14 +1061,15 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum(d_i Z_i) - IF (rkD(1)/=ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2)/=ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3)/=ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1)/=ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2)/=ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3)/=ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) + ! Update TLM solution: Y <- Y + sum(d_i*Z_i_tlm) DO itlm = 1,NTLM - IF (rkD(1)/=ZERO) CALL WAXPY(N,rkD(1),Z1_tlm(1,itlm),1,Y_tlm(1,itlm),1) - IF (rkD(2)/=ZERO) CALL WAXPY(N,rkD(2),Z2_tlm(1,itlm),1,Y_tlm(1,itlm),1) - IF (rkD(3)/=ZERO) CALL WAXPY(N,rkD(3),Z3_tlm(1,itlm),1,Y_tlm(1,itlm),1) + IF (rkD(1)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(1) * Z1_tlm(1:N,itlm) + IF (rkD(2)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(2) * Z2_tlm(1:N,itlm) + IF (rkD(3)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(3) * Z3_tlm(1:N,itlm) END DO ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) @@ -1264,27 +1265,27 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,R1,R2,R3) KPP_REAL, INTENT(INOUT), DIMENSION(N) :: R1,R2,R3 KPP_REAL, DIMENSION(N) :: F, TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 - - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 + + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS @@ -1308,39 +1309,40 @@ SUBROUTINE RK_PrepareRHS_TLM(N,H,Jac1,Jac2,Jac3,Y_tlm, & #endif KPP_REAL, DIMENSION(N) :: F, TMP - CALL WCOPY(N,Z1_tlm,1,R1,1) ! R1 <- Z1_tlm - CALL WCOPY(N,Z2_tlm,1,R2,1) ! R2 <- Z2_tlm - CALL WCOPY(N,Z3_tlm,1,R3,1) ! R3 <- Z3_tlm + R1(1:N) = Z1_tlm(1:N) ! R1 <- Z1_tlm + R2(1:N) = Z2_tlm(1:N) ! R2 <- Z2_tlm + R3(1:N) = Z3_tlm(1:N) ! R3 <- Z3_tlm - CALL WADD(N,Y_tlm,Z1_tlm,TMP) ! TMP <- Y + Z1 + TMP(1:N) = Y_tlm(1:N) + Z1_tlm(1:N) ! TMP <- Y + Z1 #ifdef FULL_ALGEBRA F = MATMUL(Jac1,TMP) #else - CALL Jac_SP_Vec ( Jac1, TMP, F ) ! F1 <- Jac(Y+Z1)*(Y_tlm+Z1_tlm) + CALL Jac_SP_Vec ( Jac1, TMP, F ) ! F1 <- Jac(Y+Z1)*(Y_tlm+Z1_tlm) #endif - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + - CALL WADD(N,Y_tlm,Z2_tlm,TMP) ! TMP <- Y + Z2 + TMP(1:N) = Y_tlm(1:N) + Z2_tlm(1:N) ! TMP <- Y + Z2 #ifdef FULL_ALGEBRA F = MATMUL(Jac2,TMP) #else - CALL Jac_SP_Vec ( Jac2, TMP, F ) ! F2 <- Jac(Y+Z2)*(Y_tlm+Z2_tlm) + CALL Jac_SP_Vec ( Jac2, TMP, F ) ! F2 <- Jac(Y+Z2)*(Y_tlm+Z2_tlm) #endif - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 - CALL WADD(N,Y_tlm,Z3_tlm,TMP) ! TMP <- Y + Z3 + TMP(1:N) = Y_tlm(1:N) + Z3_tlm(1:N) ! TMP <- Y + Z3 #ifdef FULL_ALGEBRA F = MATMUL(Jac3,TMP) #else - CALL Jac_SP_Vec ( Jac3, TMP, F ) ! F3 <- Jac(Y+Z3)*(Y_tlm+Z3_tlm) + CALL Jac_SP_Vec ( Jac3, TMP, F ) ! F3 <- Jac(Y+Z3)*(Y_tlm+Z3_tlm) #endif - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS_TLM diff --git a/int/sdirk.f90 b/int/sdirk.f90 index 490586f..9269541 100644 --- a/int/sdirk.f90 +++ b/int/sdirk.f90 @@ -22,9 +22,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WCOPY, WAXPY, & - WSCAL, WADD + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -366,7 +364,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -591,17 +589,16 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) - + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp !~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) IF (istage > 1) THEN DO j = 1, istage-1 ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) IF (StartNewton) THEN - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END IF END DO END IF @@ -613,13 +610,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) ! RHS <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N) - CALL WSCAL(N, H*rkGamma, RHS, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1) - CALL WAXPY (N, ONE, G,1, RHS,1) + RHS(1:N) = RHS(1:N) * (H * rkGamma) + RHS(1:N) = RHS(1:N) - Z(1:N,istage) + RHS(1:N) = RHS(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, RHS ) @@ -648,7 +645,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -675,9 +672,9 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF (sdMethod /= BEL) THEN ! All methods but Backward Euler - CALL Set2zero(N,TMP) + TMP(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) + IF (rkE(i)/=ZERO) TMP(1:N) = TMP(1:N) + rkE(i) * Z(1:N,i) END DO CALL SDIRK_Solve( H, N, E, IP, IER, TMP ) @@ -704,7 +701,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) DO i = 1,rkS - IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) + IF (rkD(i)/=ZERO) Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) END DO !~~~> Update scaling coefficients @@ -918,10 +915,10 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) KPP_REAL, INTENT(IN) :: E(LU_NONZERO) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) - KPP_REAL :: HGammaInv - HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + ! NOTE: This line reproduces the results of the + ! previous WAXPY call (@yantosca, 16 Oct 2025) + RHS(1:N) = RHS(1:N) * (ONE / (H * rkGamma)) #ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else diff --git a/int/sdirk4.f90 b/int/sdirk4.f90 index 7e1c04e..435cbe0 100644 --- a/int/sdirk4.f90 +++ b/int/sdirk4.f90 @@ -12,8 +12,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WAXPY, WCOPY + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -320,7 +319,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -553,11 +552,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & NewtonFactor(istage) = MAX(NewtonFactor(istage),Roundoff)**0.8d0 !~~~> STARTING VALUES FOR NEWTON ITERATION - CALL Set2zero(N,G) - CALL Set2zero(N,Z(1,istage)) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp IF (istage==1) THEN IF (FIRST.OR.NewtonReject) THEN - CALL Set2zero(N,Z(1,istage)) + Z(1:N,istage) = 0.0_dp ELSE W=ONE+rkGamma*H/Hold DO i=1,N @@ -567,12 +566,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & ELSE DO j = 1, istage-1 ! Gj(:) = sum_j Beta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkBeta(istage,j),Z(1,j),1,G,1) - ! CALL WAXPY(N,H*rkA(istage,j),FV(1,j),1,G,1) + G(1:N) = G(1:N) + rkBeta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END DO - IF (istage==5) CALL WCOPY(N,Z(1,istage),1,Yhat,1) ! Yhat(:) <- Z5(:) + IF (istage==5) Yhat(1:N) = Z(1:N,istage) ! Yhat(:) <- Z5(:) END IF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -597,7 +595,8 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & TMP(1:N) = Y(1:N) + Z(1:N,istage) CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) TMP(1:N) = G(1:N) - Z(1:N,istage) - CALL WAXPY(N,HGammaInv,TMP,1,RHS,1) ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:)) + ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:)) + RHS(1:N) = RHS(1:N) + HGammaInv * TMP(1:N) !~~~> SOLVE THE LINEAR SYSTEMS #ifdef FULL_ALGEBRA @@ -630,8 +629,9 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & END IF END IF NewtonErrOld = NewtonErr - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) ! Z(:) <-- Z(:)+RHS(:) - + ! Z(:) <-- Z(:)+RHS(:) + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) + END DO Newton !~~> END OF SIMPLIFIED NEWTON ITERATION @@ -671,13 +671,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & Hold=H !~~~> COEFFICIENTS FOR CONTINUOUS SOLUTION - CALL WAXPY(N,ONE,Z(1,5),1,Y,1) ! Y(:) <-- Y(:)+Z5(:) - CALL WCOPY(N,Z(1,5),1,Yhat,1) ! Yhat <-- Z5 + Y(1:N) = Y(1:N) + Z(1:N,5) ! Y(:) <-- Y(:)+Z5(:) + Yhat(1:N) = Z(1:N,5) ! Yhat <-- Z5 DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj - CALL Set2zero(N,CONT(1,i)) + CONT(1:N,i) = 0.0_dp DO j = 1,5 - CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) + CONT(1:N,i) = CONT(1:N,i) + rkD(i,j) * Z(1:N,j) END DO END DO @@ -997,11 +997,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JS) - DO j=1,NVAR - DO j=1,NVAR - JV(i,j) = 0.0D0 - END DO - END DO + JV(1:NVAR,1:NVAR) = 0.0d0 DO i=1,LU_NONZERO JV(LU_IROW(i),LU_ICOL(i)) = JS(i) END DO diff --git a/int/sdirk_adj.f90 b/int/sdirk_adj.f90 index 8f58fa0..59f9b53 100644 --- a/int/sdirk_adj.f90 +++ b/int/sdirk_adj.f90 @@ -20,14 +20,13 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG USE KPP_ROOT_Jacobian, ONLY: Jac_SP_Vec, JacTR_SP_Vec - USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, & - KppSolveTR, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD - + USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, KppSolveTR + IMPLICIT NONE PUBLIC SAVE -!~~~> Flags to determine if we should call the UPDATE_* routines from within +!~~~> Flags to determine if we should call the UPDATE_* routines from within !~~~> the integrator. If using KPP in an external model, you might want to !~~~> disable these calls (via ICNTRL(15)) to avoid excess computations. LOGICAL, PRIVATE :: Do_Update_RCONST @@ -38,7 +37,7 @@ MODULE KPP_ROOT_Integrator INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, & Nrej=5, Ndec=6, Nsol=7, Nsng=8, & Ntexit=1, Nhexit=2, Nhnew=3 - + CONTAINS SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & @@ -56,7 +55,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & !~~~> NADJ - No. of cost functionals for which adjoints ! are evaluated simultaneously ! If single cost functional is considered (like in -! most applications) simply set NADJ = 1 +! most applications) simply set NADJ = 1 INTEGER :: NADJ !~~~> Lambda - Sensitivities w.r.t. concentrations ! Note: Lambda (1:NVAR,j) contains sensitivities of @@ -64,7 +63,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & KPP_REAL :: Lambda(NVAR,NADJ) !~~~> Tolerances for adjoint calculations ! (used for full continuous adjoint, and for controlling -! iterations when used to solve the discrete adjoint) +! iterations when used to solve the discrete adjoint) KPP_REAL, INTENT(IN) :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ) KPP_REAL, INTENT(IN) :: TIN ! Start Time KPP_REAL, INTENT(IN) :: TOUT ! End Time @@ -89,7 +88,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & ICNTRL(5) = 8 ! Max no. of Newton iterations ICNTRL(7) = 1 ! Adjoint solution by: 0=Newton, 1=direct ICNTRL(8) = 1 ! Save fwd LU factorization: 0 = do *not* save, 1 = save - ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. + ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. !~~~> if optional parameters are given, and if they are /= 0, !~~~> then use them to overwrite default settings @@ -99,7 +98,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & IF ( PRESENT( RCNTRL_U ) ) THEN WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U ENDIF - + !~~~> Determine the settings of the Do_Update_* flags, which determine !~~~> whether or not we need to call Update_* routines in the integrator !~~~> (or not, if we are calling them from a higher-level) @@ -109,7 +108,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & ! = 2 ! Call Update_PHOTO from within the integrator ! = 3 ! Call Update_RCONST and Update_PHOTO from w/in the int. ! = 4 ! Call Update_SUN from within the integrator - ! = 5 ! Call Update_SUN and Update_RCONST from within the int. + ! = 5 ! Call Update_SUN and Update_RCONST from within the int. ! = 6 ! Call Update_SUN and Update_PHOTO from within the int. ! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int. CALL Integrator_Update_Options( ICNTRL(15), & @@ -136,12 +135,12 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & !~~~> Debug option: number of steps ! Ntotal = Ntotal + ISTATUS(Nstp) ! WRITE(6,777) ISTATUS(Nstp),Ntotal,VAR(ind_O3),VAR(ind_NO2) - ! 777 FORMAT('NSTEPS=',I5,' (',I5,') O3=',E24.14,' NO2=',E24.14) + ! 777 FORMAT('NSTEPS=',I5,' (',I5,') O3=',E24.14,' NO2=',E24.14) IF (Ierr < 0) THEN PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (Ierr=',Ierr,')' ENDIF - + ! if optional parameters are given for output ! use them to store information in them IF ( PRESENT( ISTATUS_U ) ) ISTATUS_U = ISTATUS @@ -168,14 +167,14 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ! This code is based on the SDIRK4 routine in the above book. ! ! Methods: -! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 -! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant -! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 +! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 +! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant +! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! ! (C) Adrian Sandu, July 2005 ! Virginia Polytechnic Institute and State University ! Contact: sandu@cs.vt.edu -! Revised by Philipp Miehe and Adrian Sandu, May 2006 +! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! This implementation is part of KPP - the Kinetic PreProcessor !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! @@ -210,7 +209,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. -!~~~> +!~~~> ! ICNTRL(1) = not used ! ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors @@ -299,16 +298,16 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE -! Arguments +! Arguments INTEGER, INTENT(IN) :: N, NADJ, ICNTRL(20) KPP_REAL, INTENT(INOUT) :: Y(NVAR), Lambda(NVAR,NADJ) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, & RelTol(NVAR), AbsTol(NVAR), RCNTRL(20), & RelTol_adj(NVAR,NADJ), AbsTol_adj(NVAR,NADJ) INTEGER, INTENT(OUT) :: Ierr - INTEGER, INTENT(INOUT) :: ISTATUS(20) + INTEGER, INTENT(INOUT) :: ISTATUS(20) KPP_REAL, INTENT(OUT) :: RSTATUS(20) - + !~~~> SDIRK method coefficients, up to 5 stages INTEGER, PARAMETER :: Smax = 5 INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5 @@ -328,16 +327,16 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & #else KPP_REAL, DIMENSION(:,:), POINTER :: chk_J #endif -! Local variables +! Local variables KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, & FacMin, Facmax, FacSafe, FacRej, & ThetaMin, NewtonTol, Qmin, Qmax - LOGICAL :: SaveLU, DirectADJ + LOGICAL :: SaveLU, DirectADJ INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 - + !~~~> Initialize statistics ISTATUS(1:20) = 0 RSTATUS(1:20) = ZERO @@ -351,7 +350,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ITOL = 0 END IF -!~~~> ICNTRL(3) - method selection +!~~~> ICNTRL(3) - method selection SELECT CASE (ICNTRL(3)) CASE (0,1) CALL Sdirk2a @@ -366,7 +365,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & CASE DEFAULT CALL Sdirk2a END SELECT - + !~~~> The maximum number of time steps admitted IF (ICNTRL(4) == 0) THEN Max_no_steps = 200000 @@ -376,7 +375,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4) CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr) END IF - + !~~~> The maximum number of Newton iterations admitted IF(ICNTRL(5) == 0)THEN NewtonMaxit=8 @@ -388,14 +387,14 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & END IF END IF -!~~~> Solve ADJ equations directly or by Newton iterations +!~~~> Solve ADJ equations directly or by Newton iterations DirectADJ = (ICNTRL(7) == 1) - + !~~~> Save or not the forward LU factorization SaveLU = (ICNTRL(8) /= 0) .AND. (.NOT.DirectADJ) !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -406,7 +405,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Upper bound on the step size: (positive value) IF (RCNTRL(2) == ZERO) THEN Hmax = ABS(Tfinal-Tinitial) @@ -416,7 +415,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN Hstart = MAX(Hmin,Roundoff) @@ -426,7 +425,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax IF (RCNTRL(4) == ZERO) THEN FacMin = 0.2_dp @@ -505,19 +504,19 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & END IF END DO END IF - + IF (Ierr < 0) RETURN - -!~~~> Allocate memory buffers + +!~~~> Allocate memory buffers CALL SDIRK_AllocBuffers -!~~~> Call forward integration +!~~~> Call forward integration CALL SDIRK_FwdInt( N, Tinitial, Tfinal, Y, Ierr ) -!~~~> Call adjoint integration +!~~~> Call adjoint integration CALL SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) -!~~~> Free memory buffers +!~~~> Free memory buffers CALL SDIRK_FreeBuffers @@ -533,26 +532,26 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) USE KPP_ROOT_Parameters IMPLICIT NONE -!~~~> Arguments: +!~~~> Arguments: INTEGER :: N KPP_REAL, INTENT(INOUT) :: Y(NVAR) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: - KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & + +!~~~> Local variables: + KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & NewtonRate, SCAL(NVAR), RHS(NVAR), & T, H, Theta, Hratio, NewtonPredictedErr, & Qnewton, Err, Fac, Hnew,Tdirection, & NewtonIncrement, NewtonIncrementOld INTEGER :: j, ISING, istage, NewtonIter, IP(NVAR) LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA + +#ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: FJAC, E -#else +#else KPP_REAL, DIMENSION(LU_NONZERO):: FJAC, E -#endif +#endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -585,14 +584,14 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) IF (ISING /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF - END IF + END IF IF (ISTATUS(Nstp) > Max_no_steps) THEN CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN - END IF + END IF IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN - END IF + END IF stages:DO istage = 1, rkS @@ -601,47 +600,48 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) - -!~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp + +!~~~> Prepare the loop-independent part of the right-hand side IF (istage > 1) THEN - DO j = 1, istage-1 - ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) - ! Zi(:) = sum_j Alpha(i,j)*Zj(:) - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) - END DO + DO j = 1, istage-1 + ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) + + ! Zi(:) = sum_j Alpha(i,j)*Zj(:) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) + END DO END IF !~~~> Initializations for Newton iteration NewtonDone = .FALSE. Fac = 0.5d0 ! Step reduction factor if too many iterations - + NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) ! RHS <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N) - CALL WSCAL(N, H*rkGamma, RHS, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1) - CALL WAXPY (N, ONE, G,1, RHS,1) + RHS(1:N) = RHS(1:N) * (H * rkGamma) + RHS(1:N) = RHS(1:N) - Z(1:N,istage) + RHS(1:N) = RHS(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( 'N', H, N, E, IP, ISING, RHS ) - + !~~~> Check convergence of Newton iterations CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 + NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process + ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) IF (NewtonPredictedErr >= NewtonTol) THEN @@ -656,14 +656,14 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) - + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT NewtonLoop - + END DO NewtonLoop - + IF (.NOT.NewtonDone) THEN !CALL RK_ErrorMsg(-12,T,H,Ierr); H = Fac*H; Reject=.TRUE. @@ -673,7 +673,7 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> End of implified Newton iterations - + END DO stages @@ -681,10 +681,10 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - CALL Set2zero(N,TMP) + TMP(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) - END DO + IF (rkE(i)/=ZERO) TMP(1:N) = TMP(1:N) + rkE(i) * Z(1:N,i) + END DO CALL SDIRK_Solve( 'N', H, N, E, IP, ISING, TMP ) CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err) @@ -708,10 +708,10 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Update time and solution T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO i = 1,rkS - IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) - END DO - + DO i = 1,rkS + IF (rkD(i)/=ZERO) Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) + END DO + !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL) @@ -746,16 +746,16 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) END IF Reject = .TRUE. SkipJac = .TRUE. - SkipLU = .FALSE. + SkipLU = .FALSE. IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1 - + END IF accept - + END DO Tloop ! Successful return Ierr = 1 - + END SUBROUTINE SDIRK_FwdInt @@ -767,34 +767,34 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) USE KPP_ROOT_Parameters IMPLICIT NONE -!~~~> Arguments: +!~~~> Arguments: INTEGER, INTENT(IN) :: N, NADJ KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ) INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: + +!~~~> Local variables: KPP_REAL :: Y(NVAR) KPP_REAL :: Z(NVAR,Smax), U(NVAR,NADJ,Smax), & - TMP(NVAR), G(NVAR), & + TMP(NVAR), G(NVAR), & NewtonRate, SCAL(NVAR), DU(NVAR), & T, H, Theta, NewtonPredictedErr, & NewtonIncrement, NewtonIncrementOld INTEGER :: j, ISING, istage, iadj, NewtonIter, & IP(NVAR), IP_adj(NVAR) LOGICAL :: Reject, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA + +#ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: E, Jac, E_adj -#else +#else KPP_REAL, DIMENSION(LU_NONZERO):: E, Jac, E_adj -#endif +#endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Time loop begins !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Tloop: DO WHILE ( stack_ptr > 0 ) - + !~~~> Recover checkpoints for stage values and vectors CALL SDIRK_Pop( T, H, Y, Z, E, IP ) @@ -806,7 +806,7 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) IF (ISING /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF - END IF + END IF stages:DO istage = rkS, 1, -1 @@ -814,8 +814,8 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) TMP(1:N) = Y(1:N) + Z(1:N,istage) CALL JAC_CHEM(T+rkC(istage)*H,TMP,Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 - - IF (DirectADJ) THEN + + IF (DirectADJ) THEN #ifdef FULL_ALGEBRA E_adj(1:N,1:N) = -Jac(1:N,1:N) DO j=1,N @@ -838,61 +838,61 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) END IF adj: DO iadj = 1, NADJ - + !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(ITOL, AbsTol_adj(1:NVAR,iadj), & RelTol_adj(1:NVAR,iadj), Lambda(1:NVAR,iadj), SCAL) - + !~~~> Prepare the loop-independent part of the right-hand side ! G(:) = H*Jac^T*( B(i)*Lambda + sum_j A(j,i)*Uj(:) ) G(1:N) = rkB(istage)*Lambda(1:N,iadj) IF (istage < rkS) THEN DO j = istage+1, rkS - CALL WAXPY(N,rkA(j,istage),U(1,iadj,j),1,G,1) + G(1:N) = G(1:N) + rkA(j,istage) * U(1:N,iadj,j) END DO END IF -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac),G) ! DZ <- Jac(Y+Z)*Y_tlm -#else - CALL JacTR_SP_Vec ( Jac, G, TMP ) -#endif +#else + CALL JacTR_SP_Vec ( Jac, G, TMP ) +#endif G(1:N) = H*TMP(1:N) -DirADJ:IF (DirectADJ) THEN +DirADJ:IF (DirectADJ) THEN CALL SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, ISING, G ) U(1:N,iadj,istage) = G(1:N) - + ELSE DirADJ !~~~> Initializations for Newton iteration - CALL Set2zero(N,U(1,iadj,istage)) + U(1:N,iadj,istage) = 0.0_dp NewtonDone = .FALSE. - + NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side -#ifdef FULL_ALGEBRA - TMP = MATMUL(TRANSPOSE(Jac),U(1:N,iadj,istage)) -#else - CALL JacTR_SP_Vec ( Jac, U(1:N,iadj,istage), TMP ) -#endif +#ifdef FULL_ALGEBRA + TMP = MATMUL(TRANSPOSE(Jac),U(1:N,iadj,istage)) +#else + CALL JacTR_SP_Vec ( Jac, U(1:N,iadj,istage), TMP ) +#endif DU(1:N) = U(1:N,iadj,istage) - (H*rkGamma)*TMP(1:N) - G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( 'T', H, N, E, IP, ISING, DU ) - + !~~~> Check convergence of Newton iterations - + CALL SDIRK_ErrorNorm(N, DU, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 + NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process + ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) ! Non-convergence of Newton: predicted error too large @@ -904,18 +904,18 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) NewtonIncrementOld = NewtonIncrement ! Update solution U(1:N,iadj,istage) = U(1:N,iadj,istage) - DU(1:N) - + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) ! AbsTol is often inappropriate for adjoints - ! we do at least 4 Newton iterations to ensure convergence ! of all adjoint components IF ((NewtonIter>=4) .AND. NewtonDone) EXIT NewtonLoop - + END DO NewtonLoop - + !~~~> If Newton iterations fail employ the direct solution - IF (.NOT.NewtonDone) THEN + IF (.NOT.NewtonDone) THEN PRINT*,'Problems with Newton Adjoint!!!' #ifdef FULL_ALGEBRA E_adj(1:N,1:N) = -Jac(1:N,1:N) @@ -938,31 +938,30 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) END IF CALL SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, ISING, G ) U(1:N,iadj,istage) = G(1:N) - + END IF !~~~> End of simplified Newton iterations - + END IF DirADJ END DO adj - + END DO stages !~~~> Update adjoint solution ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO istage = 1,rkS + DO istage = 1,rkS DO iadj = 1,NADJ Lambda(1:N,iadj) = Lambda(1:N,iadj) + U(1:N,iadj,istage) - !CALL WAXPY(N,ONE,U(1:N,iadj,istage),1,Lambda(1,iadj),1) - END DO - END DO + END DO + END DO END DO Tloop ! Successful return Ierr = 1 - + END SUBROUTINE SDIRK_DadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -970,23 +969,23 @@ SUBROUTINE SDIRK_AllocBuffers !~~~> Allocate buffer space for checkpointing !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i - + ALLOCATE( chk_H(Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer H'; STOP - END IF + END IF ALLOCATE( chk_T(Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer T'; STOP - END IF + END IF ALLOCATE( chk_Y(NVAR,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer Y'; STOP - END IF + END IF ALLOCATE( chk_Z(NVAR,rkS,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer K'; STOP - END IF + END IF IF (SaveLU) THEN #ifdef FULL_ALGEBRA ALLOCATE( chk_J(NVAR,NVAR,Max_no_steps), STAT=i ) @@ -995,13 +994,13 @@ SUBROUTINE SDIRK_AllocBuffers #endif IF (i/=0) THEN PRINT*,'Failed allocation of buffer J'; STOP - END IF + END IF ALLOCATE( chk_P(NVAR,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer P'; STOP - END IF - END IF - + END IF + END IF + END SUBROUTINE SDIRK_AllocBuffers @@ -1010,34 +1009,34 @@ SUBROUTINE SDIRK_FreeBuffers !~~~> Dallocate buffer space for discrete adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i - + DEALLOCATE( chk_H, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer H'; STOP - END IF + END IF DEALLOCATE( chk_T, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer T'; STOP - END IF + END IF DEALLOCATE( chk_Y, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer Y'; STOP - END IF + END IF DEALLOCATE( chk_Z, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer K'; STOP - END IF + END IF IF (SaveLU) THEN DEALLOCATE( chk_J, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer J'; STOP - END IF + END IF DEALLOCATE( chk_P, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer P'; STOP - END IF - END IF - + END IF + END IF + END SUBROUTINE SDIRK_FreeBuffers @@ -1047,35 +1046,35 @@ SUBROUTINE SDIRK_Push( T, H, Y, Z, E, P ) !~~~> Saves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - KPP_REAL :: T, H, Y(NVAR), Z(NVAR,Smax) + KPP_REAL :: T, H, Y(NVAR), Z(NVAR,Smax) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) -#else +#else KPP_REAL :: E(LU_NONZERO) #endif - + stack_ptr = stack_ptr + 1 IF ( stack_ptr > Max_no_steps ) THEN PRINT*,'Push failed: buffer overflow' STOP - END IF + END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_Z(1:NVAR,1:rkS,stack_ptr) = Z(1:NVAR,1:rkS) - IF (SaveLU) THEN + IF (SaveLU) THEN #ifdef FULL_ALGEBRA chk_J(1:NVAR,1:NVAR,stack_ptr) = E(1:NVAR,1:NVAR) chk_P(1:NVAR,stack_ptr) = P(1:NVAR) -#else +#else chk_J(1:LU_NONZERO,stack_ptr) = E(1:LU_NONZERO) #endif END IF - + END SUBROUTINE SDIRK_Push - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1086,14 +1085,14 @@ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) -#else +#else KPP_REAL :: E(LU_NONZERO) #endif - + IF ( stack_ptr <= 0 ) THEN PRINT*,'Pop failed: empty buffer' STOP - END IF + END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) @@ -1102,13 +1101,13 @@ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) #ifdef FULL_ALGEBRA E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr) P(1:NVAR) = chk_P(1:NVAR,stack_ptr) -#else +#else E(1:LU_NONZERO) = chk_J(1:LU_NONZERO,stack_ptr) #endif END IF stack_ptr = stack_ptr - 1 - + END SUBROUTINE SDIRK_Pop !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1128,14 +1127,14 @@ SUBROUTINE SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL) END DO END IF END SUBROUTINE SDIRK_ErrorScale - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! +! INTEGER :: i, N - KPP_REAL :: Y(N), SCAL(N), Err + KPP_REAL :: Y(N), SCAL(N), Err Err = ZERO DO i=1,N Err = Err+(Y(i)*SCAL(i))**2 @@ -1145,7 +1144,7 @@ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) END SUBROUTINE SDIRK_ErrorNorm -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Handles all error messages @@ -1183,17 +1182,17 @@ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) PRINT *, "T=", T, "and H=", H END SUBROUTINE SDIRK_ErrorMsg - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipJac, SkipLU, E, IP, Reject, ISING ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + IMPLICIT NONE - + KPP_REAL, INTENT(INOUT) :: H KPP_REAL, INTENT(IN) :: T, Y(NVAR) LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject @@ -1210,9 +1209,9 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & ConsecutiveSng = 0 ISING = 1 - + Hloop: DO WHILE (ISING /= 0) - + HGammaInv = ONE/(H*rkGamma) !~~~> Compute the Jacobian @@ -1222,8 +1221,8 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & IF (.NOT. SkipJac) THEN CALL JAC_CHEM( T, Y, FJAC ) ISTATUS(Njac) = ISTATUS(Njac) + 1 - END IF - + END IF + #ifdef FULL_ALGEBRA DO j=1,NVAR DO i=1,NVAR @@ -1253,7 +1252,7 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipLU = .FALSE. Reject = .TRUE. END IF - + END DO Hloop END SUBROUTINE SDIRK_PrepareMatrix @@ -1276,18 +1275,18 @@ SUBROUTINE SDIRK_Solve ( Transp, H, N, E, IP, ISING, RHS ) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) KPP_REAL :: HGammaInv - + HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + RHS(1:N) = RHS(1:N) * HGammaInv SELECT CASE (TRANSP) CASE ('N') -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolve(E, RHS) #endif CASE ('T') -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA CALL DGETRS( 'T', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolveTR(E, RHS, RHS) @@ -1297,7 +1296,7 @@ SUBROUTINE SDIRK_Solve ( Transp, H, N, E, IP, ISING, RHS ) STOP END SELECT ISTATUS(Nsol) = ISTATUS(Nsol) + 1 - + END SUBROUTINE SDIRK_Solve @@ -1385,7 +1384,7 @@ SUBROUTINE Sdirk4a rkAlpha(5,2) = 6.559571569643355712998131800797873d0 rkAlpha(5,3) = -15.90772144271326504260996815012482d0 rkAlpha(5,4) = 25.34908987169226073668861694892683d0 - + !~~~> Coefficients for continuous solution ! rkD(1,1)= 24.74416644927758d0 ! rkD(1,2)= -4.325375951824688d0 @@ -1414,9 +1413,9 @@ SUBROUTINE Sdirk4a ! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) ! END DO ! END DO - + rkELO = 4.0d0 - + END SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1503,7 +1502,7 @@ SUBROUTINE Sdirk4b rkAlpha(5,2) = 9.d0 rkAlpha(5,3) = -56.81818181818181818181818181818182d0 rkAlpha(5,4) = 54.d0 - + END SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1545,7 +1544,7 @@ SUBROUTINE Sdirk2a ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 3.414213562373095048801688724209698d0 - + END SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1587,7 +1586,7 @@ SUBROUTINE Sdirk2b ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = .5857864376269049511983112757903019d0 - + END SUBROUTINE Sdirk2b @@ -1661,16 +1660,16 @@ SUBROUTINE FUN_CHEM( T, Y, P ) KPP_REAL :: T, Told KPP_REAL :: Y(NVAR), P(NVAR) - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() IF ( Do_Update_RCONST ) CALL Update_RCONST() - + CALL Fun( Y, FIX, RCONST, P ) - + TIME = Told - + END SUBROUTINE FUN_CHEM @@ -1691,7 +1690,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #else KPP_REAL :: JV(LU_NONZERO) #endif - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() @@ -1699,11 +1698,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JS) - DO j=1,NVAR - DO j=1,NVAR - JV(i,j) = 0.0d0 - END DO - END DO + JV(1:NVAR,1:NVAR) = 0.0d0 DO i=1,LU_NONZERO JV(LU_IROW(i),LU_ICOL(i)) = JS(i) END DO @@ -1715,5 +1710,3 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) END SUBROUTINE JAC_CHEM END MODULE KPP_ROOT_Integrator - - diff --git a/int/sdirk_tlm.f90 b/int/sdirk_tlm.f90 index 92d161f..1605570 100644 --- a/int/sdirk_tlm.f90 +++ b/int/sdirk_tlm.f90 @@ -20,9 +20,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG USE KPP_ROOT_Jacobian, ONLY : Jac_SP_Vec - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WCOPY, WAXPY, & - WSCAL, WADD + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -396,7 +394,7 @@ SUBROUTINE SdirkTLM(N, NTLM, Tinitial, Tfinal, Y, Y_tlm, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -596,17 +594,17 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp -!~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) +!~~~> Prepare the loop-independent part of the right-hand side IF (istage > 1) THEN DO j = 1, istage-1 ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) IF (StartNewton) THEN - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END IF END DO END IF @@ -618,13 +616,13 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,DZ) ! DZ <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*DZ(1:N) - CALL WSCAL(N, H*rkGamma, DZ, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, DZ, 1) - CALL WAXPY (N, ONE, G,1, DZ,1) + DZ(1:N) = DZ(1:N) * (H*rkGamma) + DZ(1:N) = DZ(1:N) - Z(1:N,istage) + DZ(1:N) = DZ(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, DZ ) @@ -653,7 +651,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+DZ(:) - CALL WAXPY(N,ONE,DZ,1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + DZ(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -693,7 +691,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 - CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z_tlm(1:N,j,itlm) END DO END IF CALL SDIRK_Solve ( H, N, E_tlm, IP_tlm, IER, G ) @@ -712,7 +710,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) NewtonRate = MAX(NewtonRate,Roundoff)**0.8d0 !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z_tlm(1,istage,itlm)) + Z_tlm(1:N,istage,itlm) = 0.0_dp !~~~> Prepare the loop-independent part of the right-hand side #ifdef FULL_ALGEBRA @@ -725,7 +723,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 - CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z_tlm(1:N,j,itlm) END DO END IF @@ -778,7 +776,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) END IF !(TLMNewtonEst) ! Update solution: Z_tlm(:) <-- Z_tlm(:)+DZ(:) - CALL WAXPY(N,ONE,DZ,1,Z_tlm(1,istage,itlm),1) + Z_tlm(1:N,istage,itlm) = Z_tlm(1:N,istage,itlm) + DZ(1:N) ! Check error in Newton iterations IF (TLMNewtonEst) THEN @@ -809,19 +807,20 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - CALL Set2zero(N,Yerr) + Yerr(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,Yerr,1) + IF (rkE(i)/=ZERO) Yerr(1:N) = Yerr(1:N) + rkE(i) * Z(1:N,i) END DO CALL SDIRK_Solve ( H, N, E, IP, IER, Yerr ) CALL SDIRK_ErrorNorm(N, Yerr, SCAL, Err) IF (TLMtruncErr) THEN - CALL Set2zero(NVAR*NTLM,Yerr_tlm) + Yerr_tlm(1:N,1:NTLM) = 0.0_dp DO itlm=1,NTLM DO j=1,rkS - IF (rkE(j) /= ZERO) CALL WAXPY(N,rkE(j),Z_tlm(1,j,itlm),1,Yerr_tlm(1,itlm),1) + IF (rkE(j) /= ZERO) Yerr_tlm(1:N,itlm) = Yerr_tlm(1:N,itlm) & + + rkE(j) * Z_tlm(1:N,j,itlm) END DO CALL SDIRK_Solve (H, N, E, IP, IER, Yerr_tlm(1,itlm)) END DO @@ -846,9 +845,10 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) DO i = 1,rkS IF (rkD(i)/=ZERO) THEN - CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) + Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) DO itlm = 1, NTLM - CALL WAXPY(N,rkD(i),Z_tlm(1,i,itlm),1,Y_tlm(1,itlm),1) + Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) & + + rkD(i) * Z_tlm(1:N,i,itlm) END DO END IF END DO @@ -1081,10 +1081,10 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) KPP_REAL, INTENT(IN) :: E(LU_NONZERO) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) - KPP_REAL :: HGammaInv - HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + ! This code replicates the output of the previous + ! call to WAXPY (@yantosca, 16 Oct 2025) + RHS(1:N) = RHS(1:N) * (ONE / (H*rkGamma)) #ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else diff --git a/int/seulex.f90 b/int/seulex.f90 index 2b0f4fe..348f090 100644 --- a/int/seulex.f90 +++ b/int/seulex.f90 @@ -426,7 +426,7 @@ SUBROUTINE ATMSEULEX( N,Tinitial,Tfinal,Y,H,RelTol,AbsTol, & !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN diff --git a/util/blas.f90 b/util/blas.f90 index cd8f337..06c4fa5 100644 --- a/util/blas.f90 +++ b/util/blas.f90 @@ -8,225 +8,21 @@ ! Virginia Polytechnic Institute and State University !-------------------------------------------------------------- - -!-------------------------------------------------------------- - SUBROUTINE WCOPY(N,X,incX,Y,incY) -!-------------------------------------------------------------- -! copies a vector, x, to a vector, y: y <- x -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1) -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - INTEGER :: i,incX,incY,M,MP1,N - KPP_REAL :: X(N),Y(N) - - IF (N.LE.0) RETURN - - M = MOD(N,8) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = X(i) - END DO - IF( N .LT. 8 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,8 - Y(i) = X(i) - Y(i + 1) = X(i + 1) - Y(i + 2) = X(i + 2) - Y(i + 3) = X(i + 3) - Y(i + 4) = X(i + 4) - Y(i + 5) = X(i + 5) - Y(i + 6) = X(i + 6) - Y(i + 7) = X(i + 7) - END DO - - END SUBROUTINE WCOPY - - -!-------------------------------------------------------------- - SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY) -!-------------------------------------------------------------- -! constant times a vector plus a vector: y <- y + Alpha*x -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1) -!-------------------------------------------------------------- - - INTEGER :: i,incX,incY,M,MP1,N - KPP_REAL :: X(N),Y(N),Alpha - KPP_REAL, PARAMETER :: ZERO = 0.0_dp - - IF (Alpha .EQ. ZERO) RETURN - IF (N .LE. 0) RETURN - - M = MOD(N,4) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = Y(i) + Alpha*X(i) - END DO - IF( N .LT. 4 ) RETURN - END IF - MP1 = M + 1 - DO i = MP1,N,4 - Y(i) = Y(i) + Alpha*X(i) - Y(i + 1) = Y(i + 1) + Alpha*X(i + 1) - Y(i + 2) = Y(i + 2) + Alpha*X(i + 2) - Y(i + 3) = Y(i + 3) + Alpha*X(i + 3) - END DO - - END SUBROUTINE WAXPY - - - -!-------------------------------------------------------------- - SUBROUTINE WSCAL(N,Alpha,X,incX) -!-------------------------------------------------------------- -! constant times a vector: x(1:N) <- Alpha*x(1:N) -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1) -!-------------------------------------------------------------- - - INTEGER :: i,incX,M,MP1,N - KPP_REAL :: X(N),Alpha - KPP_REAL, PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp - - IF (Alpha .EQ. ONE) RETURN - IF (N .LE. 0) RETURN - - M = MOD(N,5) - IF( M .NE. 0 ) THEN - IF (Alpha .EQ. (-ONE)) THEN - DO i = 1,M - X(i) = -X(i) - END DO - ELSEIF (Alpha .EQ. ZERO) THEN - DO i = 1,M - X(i) = ZERO - END DO - ELSE - DO i = 1,M - X(i) = Alpha*X(i) - END DO - END IF - IF( N .LT. 5 ) RETURN - END IF - MP1 = M + 1 - IF (Alpha .EQ. (-ONE)) THEN - DO i = MP1,N,5 - X(i) = -X(i) - X(i + 1) = -X(i + 1) - X(i + 2) = -X(i + 2) - X(i + 3) = -X(i + 3) - X(i + 4) = -X(i + 4) - END DO - ELSEIF (Alpha .EQ. ZERO) THEN - DO i = MP1,N,5 - X(i) = ZERO - X(i + 1) = ZERO - X(i + 2) = ZERO - X(i + 3) = ZERO - X(i + 4) = ZERO - END DO - ELSE - DO i = MP1,N,5 - X(i) = Alpha*X(i) - X(i + 1) = Alpha*X(i + 1) - X(i + 2) = Alpha*X(i + 2) - X(i + 3) = Alpha*X(i + 3) - X(i + 4) = Alpha*X(i + 4) - END DO - END IF - - END SUBROUTINE WSCAL - -!-------------------------------------------------------------- - KPP_REAL FUNCTION WLAMCH( C ) -!-------------------------------------------------------------- -! returns epsilon machine -! after LAPACK -! replace this by the function from the optimized LAPACK implementation: -! CALL SLAMCH('E') or CALL DLAMCH('E') -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - CHARACTER :: C - INTEGER :: i - KPP_REAL, SAVE :: Eps - KPP_REAL :: Suma - KPP_REAL, PARAMETER :: ONE=1.0_dp, HALF=0.5_dp - LOGICAL, SAVE :: First=.TRUE. - -!$OMP THREADPRIVATE( Eps, First ) - - IF (First) THEN - First = .FALSE. - Eps = HALF**(16) - DO i = 17, 80 - Eps = Eps*HALF - CALL WLAMCH_ADD(ONE,Eps,Suma) - IF (Suma.LE.ONE) GOTO 10 - END DO - PRINT*,'ERROR IN WLAMCH. EPS < ',Eps - RETURN -10 Eps = Eps*2 - i = i-1 - END IF - - WLAMCH = Eps - - END FUNCTION WLAMCH - - SUBROUTINE WLAMCH_ADD( A, B, Suma ) -! USE KPP_ROOT_Precision - - KPP_REAL A, B, Suma - Suma = A + B - - END SUBROUTINE WLAMCH_ADD -!-------------------------------------------------------------- - - -!-------------------------------------------------------------- - SUBROUTINE SET2ZERO(N,Y) -!-------------------------------------------------------------- -! copies zeros into the vector y: y <- 0 -! after BLAS -!-------------------------------------------------------------- - - INTEGER :: i,M,MP1,N - KPP_REAL :: Y(N) - KPP_REAL, PARAMETER :: ZERO = 0.0d0 - - IF (N.LE.0) RETURN - - M = MOD(N,8) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = ZERO - END DO - IF( N .LT. 8 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,8 - Y(i) = ZERO - Y(i + 1) = ZERO - Y(i + 2) = ZERO - Y(i + 3) = ZERO - Y(i + 4) = ZERO - Y(i + 5) = ZERO - Y(i + 6) = ZERO - Y(i + 7) = ZERO - END DO - - END SUBROUTINE SET2ZERO - +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +!%%% NOTE: The following BLAS functions have been removed, %%% +!%%% as they now have been replaced by pure F90 code %%% +!%%% in the various integrator modules; %%% +!%%% %%% +!%%% (1) WCOPY %%% +!%%% (2) WAXPY %%% +!%%% (3) WSCAL %%% +!%%% (4) WLAMCH %%% +!%%% (5) WLAMCH_ADD %%% +!%%% (6) SET2ZERO %%% +!%%% (7) WADD %%% +!%%% %%% +!%%% @yantosca, 17 Oct 2025 %%% +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------------------------------------------------- KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) @@ -288,40 +84,6 @@ KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) END FUNCTION WDOT - -!-------------------------------------------------------------- - SUBROUTINE WADD(N,X,Y,Z) -!-------------------------------------------------------------- -! adds two vectors: z <- x + y -! BLAS - like -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - INTEGER :: i, M, MP1, N - KPP_REAL :: X(N),Y(N),Z(N) - - IF (N.LE.0) RETURN - - M = MOD(N,5) - IF( M /= 0 ) THEN - DO i = 1,M - Z(i) = X(i) + Y(i) - END DO - IF( N < 5 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,5 - Z(i) = X(i) + Y(i) - Z(i + 1) = X(i + 1) + Y(i + 1) - Z(i + 2) = X(i + 2) + Y(i + 2) - Z(i + 3) = X(i + 3) + Y(i + 3) - Z(i + 4) = X(i + 4) + Y(i + 4) - END DO - - END SUBROUTINE WADD - - - !-------------------------------------------------------------- SUBROUTINE WGEFA(N,A,Ipvt,info) !-------------------------------------------------------------- @@ -333,7 +95,7 @@ SUBROUTINE WGEFA(N,A,Ipvt,info) INTEGER :: N,Ipvt(N),info KPP_REAL :: A(N,N) KPP_REAL :: t, dmax, da - INTEGER :: j,k,l + INTEGER :: i,j,k,l KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0 info = 0 @@ -362,13 +124,19 @@ SUBROUTINE WGEFA(N,A,Ipvt,info) t = A(l,k); A(l,k) = A(k,k); A(k,k) = t END IF t = -ONE/A(k,k) - CALL WSCAL(n-k,t,A(k+1,k),1) +! CALL WSCAL(n-k,t,A(k+1,k),1) + DO i = k+1, n + A(i,k) = t * A(i,k) + END DO DO j = k+1, n t = A(l,j) IF (l /= k) THEN A(l,j) = A(k,j); A(k,j) = t END IF - CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) + !CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) + DO i = k+1, n + A(i,j) = A(i,j) + t * A(i,k) + END DO END DO END IF @@ -398,7 +166,7 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) CHARACTER :: trans KPP_REAL :: A(N,N),b(N) KPP_REAL :: t - INTEGER :: k,kb,l + INTEGER :: i, k,kb,l SELECT CASE (Trans) @@ -414,7 +182,10 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) b(l) = b(k) b(k) = t END IF - CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) + !CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) + DO i = k+1, n + b(i) = b(i) + t * a(i,k) + END DO END DO END IF ! now solve U*x = y @@ -422,7 +193,10 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) - CALL WAXPY(k-1,t,a(1,k),1,b(1),1) + !CALL WAXPY(k-1,t,a(1,k),1,b(1),1) + DO i = 1, k-1 + b(i) = b(i) + t * a(i,k) + END DO END DO CASE ('t','T') ! Solve transpose(A) * x = b