Skip to content

Conversation

@yantosca
Copy link
Contributor

@yantosca yantosca commented Oct 29, 2025

This PR replaces the BLAS WDOT function from util/blas.F90 with the Fortran90 intrinsic function DOT_PRODUCT, which does the same thing, and probably more efficiently.

We have also updated a couple of calls in BLAS routine WGESL as follows:

Prior code:

      CASE ('t','T')  !  Solve transpose(A) * x = b

!        first solve  trans(U)*y = b
         DO k = 1, n
            t = WDOT(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
         END DO
!        now solve trans(L)*x = y
         IF (n >= 2) THEN
         DO kb = 1, n-1
            k = n - kb
            b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
            l = Ipvt(k)
            IF (l /= k) THEN
               t = b(l); b(l) = b(k); b(k) = t
            END IF
         END DO
         END IF

New code:

      CASE ('t','T')  !  Solve transpose(A) * x = b

!        first solve  trans(U)*y = b
         DO k = 1, n
            t = DOT_PRODUCT( a(1:k-1, k), b(1:k-1) )
            b(k) = (b(k) - t)/a(k,k)
         END DO
!        now solve trans(L)*x = y
         IF (n >= 2) THEN
         DO kb = 1, n-1
            k = n - kb
            b(k) = b(k) + DOT_PRODUCT( a(k+1:n, k), b(k+1:n) )
            l = Ipvt(k)
            IF (l /= k) THEN
               t = b(l); b(l) = b(k); b(k) = t
            END IF
         END DO
         END IF

As a side benefit, this PR finally removes this compiler warning, caused by a Fortran66-style IF statement in WDOT:

gfortran -cpp -O -g  -c F90_feuler_Initialize.F90
F90_feuler_LinearAlgebra.F90:86:46:

   86 |       IF (incX .EQ. incY) IF (incX-1) 5,20,60
      |                                              1
Warning: Fortran 2018 deleted feature: Arithmetic IF statement at (1)
      IF (incX .EQ. incY) IF (incX-1) 5,20,60 

util/blas.f90
- Removed WDOT routine
- Replaced calls to WDOT with calls to F90 intrinsic DOT_PRODUCT

NOTE: This also has the advantage of removing the compiler warning
for the computed IF statement.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
@yantosca yantosca added this to the 3.4.0 milestone Oct 29, 2025
@yantosca yantosca requested a review from RolfSander October 29, 2025 19:14
@yantosca yantosca self-assigned this Oct 29, 2025
@yantosca yantosca added feature New feature or request integrators Related to numerical integrators labels Oct 29, 2025
@yantosca yantosca mentioned this pull request Oct 29, 2025
15 tasks
Copy link
Contributor

@RolfSander RolfSander left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great! You can also add a checkmark to "Update obsolete fortran code,
e.g., arithmetic IF statements" in "Ideas for KPP 4.0.0".

@yantosca yantosca merged commit 232610a into dev Oct 31, 2025
10 checks passed
@yantosca yantosca deleted the feature/replace-wdot-with-dot_product branch October 31, 2025 15:06
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature New feature or request integrators Related to numerical integrators

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants