Skip to content

Commit 758f078

Browse files
yantoscayantosca
committed
Replace BLAS functions in radau5, runge_kutta*, and sdirk_tlm
int/radau5.f90 int/runge_kutta.f90 int/runge_kutta_adj.f90 int/runge_kutta_tlm.f90 int/sdirk_tlm.f90 - Replace BLAS routines with pure F90 code Signed-off-by: Bob Yantosca <[email protected]>
1 parent 3de658c commit 758f078

File tree

5 files changed

+265
-251
lines changed

5 files changed

+265
-251
lines changed

‎int/radau5.f90‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -726,9 +726,9 @@ SUBROUTINE RAD_Integrator( &
726726
END IF
727727
END IF
728728
DYNOLD=MAX(DYNO,Roundoff)
729-
CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1
730-
CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2
731-
CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3
729+
F1(1:N) = F1(1:N) + Z1(1:N) ! F1 <- F1 + Z1
730+
F2(1:N) = F2(1:N) + Z2(1:N) ! F2 <- F2 + Z2
731+
F3(1:N) = F3(1:N) + Z3(1:N) ! F3 <- F3 + Z3
732732
! Z(1,2,3) = Transf x F(1,2,3)
733733
CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3)
734734
NewtonDone = (FacConv*DYNO <= TolNewton)

‎int/runge_kutta.f90‎

Lines changed: 54 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,9 @@ MODULE KPP_ROOT_Integrator
1919
USE KPP_ROOT_Precision
2020
USE KPP_ROOT_Parameters
2121
USE KPP_ROOT_Global
22-
USE KPP_ROOT_Jacobian, ONLY : LU_DIAG
23-
USE KPP_ROOT_LinearAlgebra
22+
USE KPP_ROOT_Jacobian, ONLY : LU_DIAG
23+
USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, &
24+
KppDecompCmplx, KppSolveCmplx, WLAMCH
2425

2526
IMPLICIT NONE
2627
PUBLIC
@@ -591,9 +592,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
591592

592593
!~~~> Starting values for Newton iteration
593594
IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN
594-
CALL Set2zero(N,Z1)
595-
CALL Set2zero(N,Z2)
596-
CALL Set2zero(N,Z3)
595+
Z1(1:N) = 0.0_dp
596+
Z2(1:N) = 0.0_dp
597+
Z3(1:N) = 0.0_dp
597598
ELSE
598599
! Evaluate quadratic polynomial
599600
CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT)
@@ -640,9 +641,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
640641

641642
NewtonIncrementOld = MAX(NewtonIncrement,Roundoff)
642643
! Update solution
643-
CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1
644-
CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2
645-
CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3
644+
Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1
645+
Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2
646+
Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3
646647

647648
! Check error in Newton iterations
648649
NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
@@ -670,11 +671,11 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
670671

671672
!~~~> Prepare the loop-independent part of the right-hand side
672673
! G = H*rkBgam(0)*F0 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3
673-
CALL Set2Zero(N, G)
674-
IF (rkMethod/=L3A) CALL WAXPY(N,rkBgam(0)*H, F0,1,G,1)
675-
CALL WAXPY(N,rkTheta(1),Z1,1,G,1)
676-
CALL WAXPY(N,rkTheta(2),Z2,1,G,1)
677-
CALL WAXPY(N,rkTheta(3),Z3,1,G,1)
674+
G(1:N) = 0.0_dp
675+
IF (rkMethod/=L3A) G(1:N) = G(1:N) + rkBgam(0)*H * F0(1:N)
676+
G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N)
677+
G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N)
678+
G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N)
678679

679680
!~~~> Initializations for Newton iteration
680681
NewtonDone = .FALSE.
@@ -683,12 +684,12 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
683684
SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit
684685

685686
!~~~> Prepare the loop-dependent part of the right-hand side
686-
CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4
687+
TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4
687688
CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4)
688689
ISTATUS(Nfun) = ISTATUS(Nfun) + 1
689690
! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N)
690-
CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1)
691-
CALL WAXPY (N, rkGamma/H, G,1, DZ4,1)
691+
DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N)
692+
DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N)
692693

693694
!~~~> Solve the linear system
694695
#ifdef FULL_ALGEBRA
@@ -723,8 +724,8 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
723724
END IF
724725
NewtonIncrementOld = NewtonIncrement
725726
! Update solution: Z4 <-- Z4 + DZ4
726-
CALL WAXPY(N,ONE,DZ4,1,Z4,1)
727-
727+
Z4(1:N) = Z4(1:N) + DZ4(1:N)
728+
728729
! Check error in Newton iterations
729730
NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
730731
IF (NewtonDone) EXIT SDNewtonLoop
@@ -742,21 +743,21 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
742743
!~~~> Error estimation
743744
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
744745
IF (SdirkError) THEN
745-
CALL Set2Zero(N, DZ4)
746+
DZ4(1:N) = 0.0_dp
746747
IF (rkMethod==L3A) THEN
747748
DZ4(1:N) = H*rkF(0)*F0(1:N)
748-
IF (rkF(1) /= ZERO) CALL WAXPY(N, rkF(1), Z1, 1, DZ4, 1)
749-
IF (rkF(2) /= ZERO) CALL WAXPY(N, rkF(2), Z2, 1, DZ4, 1)
750-
IF (rkF(3) /= ZERO) CALL WAXPY(N, rkF(3), Z3, 1, DZ4, 1)
749+
IF (rkF(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(1) * Z1(1:N)
750+
IF (rkF(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(2) * Z2(1:N)
751+
IF (rkF(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(3) * Z3(1:N)
751752
TMP = Y + Z4
752753
CALL FUN_CHEM(T+H,TMP,DZ1)
753-
CALL WAXPY(N, H*rkBgam(4), DZ1, 1, DZ4, 1)
754+
DZ4(1:N) = DZ4(1:N) + H*rkBgam(4) * DZ1(1:N)
754755
ELSE
755756
! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4
756-
IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1)
757-
IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1)
758-
IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1)
759-
CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1)
757+
IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N)
758+
IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N)
759+
IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N)
760+
DZ4(1:N) = DZ4(1:N) - Z4(1:N)
760761
END IF
761762
Err = RK_ErrorNorm(N,SCAL,DZ4)
762763
ELSE
@@ -790,9 +791,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
790791
Hold = H
791792
T = T+H
792793
! Update solution: Y <- Y + sum(d_i Z_i)
793-
IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1)
794-
IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1)
795-
IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1)
794+
IF (rkD(1) /= ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N)
795+
IF (rkD(2) /= ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N)
796+
IF (rkD(3) /= ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N)
796797
! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3
797798
IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT)
798799
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)
978979
KPP_REAL :: T, H
979980
KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F0,F,R1,R2,R3,TMP
980981

981-
CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1
982-
CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2
983-
CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3
982+
R1(1:N) = Z1(1:N) ! R1 <- Z1
983+
R2(1:N) = Z2(1:N) ! R2 <- Z2
984+
R3(1:N) = Z3(1:N) ! R3 <- Z3
984985

985986
IF (rkMethod==L3A) THEN
986-
CALL WAXPY(N,-H*rkA(1,0),F0,1,R1,1) ! R1 <- R1 - h*A_10*F0
987-
CALL WAXPY(N,-H*rkA(2,0),F0,1,R2,1) ! R2 <- R2 - h*A_20*F0
988-
CALL WAXPY(N,-H*rkA(3,0),F0,1,R3,1) ! R3 <- R3 - h*A_30*F0
987+
R1(1:N) = R1(1:N) - H*rkA(1,0) * F0(1:N) ! R1 <- R1 - h*A_10*F0
988+
R2(1:N) = R2(1:N) - H*rkA(2,0) * F0(1:N) ! R2 <- R2 - h*A_20*F0
989+
R3(1:N) = R3(1:N) - H*rkA(3,0) * F0(1:N) ! R3 <- R3 - h*A_30*F0
989990
END IF
990991

991-
CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1
992-
CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1)
993-
CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1
994-
CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1
995-
CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1
996-
997-
CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2
998-
CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2)
999-
CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2
1000-
CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2
1001-
CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2
1002-
1003-
CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3
1004-
CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3)
1005-
CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3
1006-
CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3
1007-
CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3
992+
TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1
993+
CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1)
994+
R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1
995+
R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1
996+
R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1
997+
998+
TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2
999+
CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2)
1000+
R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2
1001+
R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2
1002+
R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2
1003+
1004+
TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3
1005+
CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3)
1006+
R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3
1007+
R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3
1008+
R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3
10081009

10091010
END SUBROUTINE RK_PrepareRHS
10101011

0 commit comments

Comments
 (0)