Replace BLAS functions w/ pure F90 code in integrators

.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
 "

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

ci-tests/F90_rkadj/F90_rkadj.kpp

@@ -0,0 +1 @@
+../../examples/rkadj.kpp

ci-tests/F90_sd4/F90_sd4.kpp

@@ -0,0 +1 @@
+../../examples/sd4.kpp

ci-tests/F90_sdtlm/F90_sdtlm.kpp

@@ -0,0 +1 @@
+../../examples/sdtlm.kpp

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

examples/rkadj.kpp

@@ -0,0 +1,4 @@
+#MODEL      small_strato
+#INTEGRATOR runge_kutta_adj
+#LANGUAGE   Fortran90
+#DRIVER     general_adj

examples/sd4.kpp

@@ -0,0 +1,4 @@
+#MODEL      small_strato
+#INTEGRATOR sdirk4
+#LANGUAGE   Fortran90
+#DRIVER     general

examples/sdtlm.kpp

@@ -0,0 +1,4 @@
+#MODEL      small_strato
+#INTEGRATOR sdirk_tlm
+#LANGUAGE   Fortran90
+#DRIVER     general_tlm

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

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

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)

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