Skip to content

Deprecating mixed precision inner product subroutines #373

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Matthew Walker (mattatmet) wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Deprecating mixed precision inner product subroutines #373

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b60ea5e
psykal_light removal and replace calls in solver and signature
mattatmet May 19, 2026
68f7cf7
Adding my name to contributors list
mattatmet May 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion CONTRIBUTORS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,4 +22,5 @@
| MetBenjaminWent | Benjamin Went | Met Office | 2026-01-30 |
| jcsmeto | James Cunningham-Smith | Met Office | 2026-02-06 |
| thomasmelvin | Thomas Melvin | Met Office | 2026-01-15 |
| ericaneininger | Erica Neininger | Met Office | 2026-03-02 |
| ericaneininger | Erica Neininger | Met Office | 2026-03-02 |
| mattatmet | Matthew Walker | Met Office | 2026-04-21 |
12 changes: 3 additions & 9 deletions components/science/source/algorithm/solver/sci_r_solver_field_vector_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,6 @@ contains
! compute the norm of a field vector
function norm_field_vector(self) result(normal)

use sci_psykal_light_mod, only: invoke_rdouble_X_innerproduct_X

implicit none
class(r_solver_field_vector_type), intent(in) :: self
real(kind=r_def) :: normal
Expand All @@ -129,7 +127,7 @@ contains

nfctr = size(self%vector)
do fctr = 1, nfctr
call invoke_rdouble_X_innerproduct_X(field_norm, self%vector(fctr))
call invoke(X_innerproduct_X(field_norm, self%vector(fctr)))
normal=normal + field_norm
end do
normal = sqrt(normal)
Expand All @@ -138,22 +136,18 @@ contains
! compute the norm of a field in the vector
function field_norm_field_vector(self, n) result(normal)

use sci_psykal_light_mod, only: invoke_rdouble_X_innerproduct_X

implicit none
class(r_solver_field_vector_type), intent(in) :: self
integer(kind=i_def), intent(in) :: n
real(kind=r_def) :: normal

call invoke_rdouble_X_innerproduct_X(normal, self%vector(n))
call invoke(X_innerproduct_X(normal, self%vector(n)))
normal = sqrt(normal)
end function field_norm_field_vector

! compute the dot of inner product of a field vector
function dot_field_vector(self, x) result(dot_prod)

use sci_psykal_light_mod, only: invoke_rdouble_X_innerproduct_Y

implicit none
class(r_solver_field_vector_type), intent(in) :: self
class(abstract_vector_type), intent(in) :: x
Expand All @@ -167,7 +161,7 @@ contains
nfctr = size(self%vector)
do fctr = 1, nfctr
inner_prod_field = 0.0_r_def
call invoke_rdouble_X_innerproduct_Y( inner_prod_field, self%vector(fctr), x%vector(fctr) )
call invoke(X_innerproduct_Y( inner_prod_field, self%vector(fctr), x%vector(fctr)))
dot_prod = dot_prod + inner_prod_field
end do
class default
Expand Down
148 changes: 0 additions & 148 deletions components/science/source/psy/sci_psykal_light_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,154 +65,6 @@ subroutine invoke_rtran_halo_exchange(field, depth)

end subroutine invoke_rtran_halo_exchange

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Psyclone does not currently have native support for builtins with mixed
! precision, this will be addressed in https://github.com/stfc/PSyclone/issues/1786
! Perform innerproduct of a r_solver precision field in r_double precision
subroutine invoke_rdouble_X_innerproduct_X(field_norm, field)

use scalar_mod, only: scalar_type
use omp_lib, only: omp_get_thread_num
use omp_lib, only: omp_get_max_threads
use mesh_mod, only: mesh_type

implicit none

real(kind=r_def), intent(out) :: field_norm
type(r_solver_field_type), intent(in) :: field

type(scalar_type) :: global_sum
integer(kind=i_def) :: df
real(kind=r_double), allocatable, dimension(:) :: l_field_norm
integer(kind=i_def) :: th_idx
integer(kind=i_def) :: loop0_start, loop0_stop
integer(kind=i_def) :: nthreads
type(r_solver_field_proxy_type) :: field_proxy
integer(kind=i_def) :: max_halo_depth_mesh
type(mesh_type), pointer :: mesh => null()
!
! Determine the number of OpenMP threads
!
nthreads = omp_get_max_threads()
!
! Initialise field and/or operator proxies
!
field_proxy = field%get_proxy()
!
! Create a mesh object
!
mesh => field_proxy%vspace%get_mesh()
max_halo_depth_mesh = mesh%get_halo_depth()
!
! Set-up all of the loop bounds
!
loop0_start = 1
loop0_stop = field_proxy%vspace%get_last_dof_owned()
!
! Call kernels and communication routines
!
!
! Zero summation variables
!
field_norm = 0.0_r_def
ALLOCATE (l_field_norm(nthreads))
l_field_norm = 0.0_r_double
!
!$omp parallel default(shared), private(df,th_idx)
th_idx = omp_get_thread_num()+1
!$omp do schedule(static)
DO df=loop0_start,loop0_stop
l_field_norm(th_idx) = l_field_norm(th_idx) + real(field_proxy%data(df),r_double)**2
END DO
!$omp end do
!$omp end parallel
!
! sum the partial results sequentially
!
DO th_idx=1,nthreads
field_norm = field_norm+real(l_field_norm(th_idx),r_def)
END DO
DEALLOCATE (l_field_norm)
global_sum%value = field_norm
field_norm = global_sum%get_sum()
!
end subroutine invoke_rdouble_X_innerproduct_X


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Psyclone does not currently have native support for builtins with mixed
! precision, this will be addressed in https://github.com/stfc/PSyclone/issues/1786
! Perform innerproduct of a r_solver precision field in r_def precision
subroutine invoke_rdouble_X_innerproduct_Y(field_norm, field1, field2)

use scalar_mod, only: scalar_type
use omp_lib, only: omp_get_thread_num
use omp_lib, only: omp_get_max_threads
use mesh_mod, only: mesh_type

implicit none

real(kind=r_def), intent(out) :: field_norm
type(r_solver_field_type), intent(in) :: field1, field2

type(scalar_type) :: global_sum
integer(kind=i_def) :: df
real(kind=r_double), allocatable, dimension(:) :: l_field_norm
integer(kind=i_def) :: th_idx
integer(kind=i_def) :: loop0_start, loop0_stop
integer(kind=i_def) :: nthreads
type(r_solver_field_proxy_type) :: field1_proxy, field2_proxy
integer(kind=i_def) :: max_halo_depth_mesh
type(mesh_type), pointer :: mesh => null()
!
! Determine the number of OpenMP threads
!
nthreads = omp_get_max_threads()
!
! Initialise field and/or operator proxies
!
field1_proxy = field1%get_proxy()
field2_proxy = field2%get_proxy()
!
! Create a mesh object
!
mesh => field1_proxy%vspace%get_mesh()
max_halo_depth_mesh = mesh%get_halo_depth()
!
! Set-up all of the loop bounds
!
loop0_start = 1
loop0_stop = field1_proxy%vspace%get_last_dof_owned()
!
! Call kernels and communication routines
!
!
! Zero summation variables
!
field_norm = 0.0_r_def
ALLOCATE (l_field_norm(nthreads))
l_field_norm = 0.0_r_double
!
!$omp parallel default(shared), private(df,th_idx)
th_idx = omp_get_thread_num()+1
!$omp do schedule(static)
DO df=loop0_start,loop0_stop
l_field_norm(th_idx) = l_field_norm(th_idx) + real(field1_proxy%data(df),r_double)*real(field2_proxy%data(df),r_double)
END DO
!$omp end do
!$omp end parallel
!
! sum the partial results sequentially
!
DO th_idx=1,nthreads
field_norm = field_norm+real(l_field_norm(th_idx),r_def)
END DO
DEALLOCATE (l_field_norm)
global_sum%value = field_norm
field_norm = global_sum%get_sum()
!
end subroutine invoke_rdouble_X_innerproduct_Y


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine invoke_inc_rdefX_plus_rsolverY(X, Y)
Expand Down