Remove BLAS functions from int/sdirk*.f90 integrators

yantosca · yantosca · commit b02d0e46ec40 · 2025-10-15T17:37:54.000-04:00
int/sdirk.f90
int/sdirk4.f90
int/sdirk_adj.f90
- Replaced calls to BLAS functions with pure F90 code
  (as determined by AI)

Signed-off-by: Bob Yantosca &lt;yantosca@seas.harvard.edu&gt;
diff --git 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, WLAMCH
 
   IMPLICIT NONE
   PUBLIC
@@ -591,35 +589,35 @@ 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,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)
-               ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
-               IF (StartNewton) THEN
-                  CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
-               END IF
-           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(:)
+             IF (StartNewton) THEN
+                Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j)
+             END IF
+          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
+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
-            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)
+!~~~>  Prepare the loop-dependent part of the right-hand side
+          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)
+          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 +646,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 +673,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 +702,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
@@ -921,7 +919,7 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS )
       KPP_REAL                :: HGammaInv
 
       HGammaInv = ONE/(H*rkGamma)
-      CALL WSCAL(N,HGammaInv,RHS,1)
+      RHS(1:N) = RHS(1:N) * HGammaInv
 #ifdef FULL_ALGEBRA
       CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
 #else
@@ -1069,13 +1067,15 @@ SUBROUTINE Sdirk4a()
 !          rkD(4,3)= -117.6778956422581d0
 !          rkD(4,4)= 134.3962081008550d0
 !          rkD(4,5)= 8.678995715052762d0
-!
-!         DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
-!           CALL Set2zero(N,CONT(1,i))
-!           DO j = 1,rkS
-!             CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
-!           END DO
-!         END DO
+! 
+!          DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
+!             CONT(1:N,i) = 0.0_dp
+!             DO j = 1,rkS
+!                IF (rkD(i,j) /= 0.0_dp) THEN
+!                  CONT(1:N,i) = CONT(1:N,i) + rkD(i,j) * Z(1:N,j)
+!                END IF
+!             END DO
+!          END DO
 
           rkELO = 4.0d0
 
diff --git 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, Set2zero, WLAMCH
 
   IMPLICIT NONE
   PUBLIC
@@ -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
