Skip to content
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.

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

mixed precision tridiagonal #1363

Merged
merged 11 commits into from
Sep 28, 2023
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
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
2 changes: 1 addition & 1 deletion Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ endif
# Make targets will be run in each subdirectory. Order is significant.
SUBDIRS = \
platform \
tridiagonal \
mpp \
tridiagonal \
constants \
constants4 \
memutils \
Expand Down
2 changes: 1 addition & 1 deletion configure.ac
Original file line number Diff line number Diff line change
Expand Up @@ -503,12 +503,12 @@ AC_CONFIG_FILES([
test_fms/coupler/Makefile
test_fms/parser/Makefile
test_fms/string_utils/Makefile
test_fms/tridiagonal/Makefile
test_fms/sat_vapor_pres/Makefile
test_fms/diag_integral/Makefile
test_fms/tracer_manager/Makefile
test_fms/random_numbers/Makefile
test_fms/column_diagnostics/Makefile

FMS.pc
])

Expand Down
2 changes: 1 addition & 1 deletion test_fms/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ ACLOCAL_AMFLAGS = -I m4
SUBDIRS = astronomy coupler diag_manager data_override exchange monin_obukhov drifters \
mosaic interpolator fms mpp mpp_io time_interp time_manager horiz_interp \
field_manager axis_utils affinity fms2_io parser string_utils sat_vapor_pres tracer_manager \
random_numbers diag_integral column_diagnostics
random_numbers diag_integral column_diagnostics tridiagonal


# testing utility scripts to distribute
Expand Down
52 changes: 52 additions & 0 deletions test_fms/tridiagonal/Makefile.am
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
#***********************************************************************
#* GNU Lesser General Public License
#*
#* This file is part of the GFDL Flexible Modeling System (FMS).
#*
#* FMS is free software: you can redistribute it and/or modify it under
#* the terms of the GNU Lesser General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or (at
#* your option) any later version.
#*
#* FMS is distributed in the hope that it will be useful, but WITHOUT
#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
#* for more details.
#*
#* You should have received a copy of the GNU Lesser General Public
#* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
#***********************************************************************

# This is an automake file for the test_fms/tridiagonal directory of the FMS
# package.

# Ryan Mulhall

# Find the .mod directory
AM_CPPFLAGS = -I$(top_srcdir)/include -I$(MODDIR)

# Link to the FMS library.
LDADD = $(top_builddir)/libFMS/libFMS.la

# Build this test program.
check_PROGRAMS = test_tridiagonal_r4 test_tridiagonal_r8

# compiles test file with both kind sizes via macro
test_tridiagonal_r4_SOURCES=test_tridiagonal.F90
test_tridiagonal_r8_SOURCES=test_tridiagonal.F90

test_tridiagonal_r4_CPPFLAGS=-DTRID_REAL_TYPE=tridiag_r4 -DTEST_TRIDIAG_REAL=r4_kind -I$(MODDIR)
test_tridiagonal_r8_CPPFLAGS=-DTRID_REAL_TYPE=tridiag_r8 -DTEST_TRIDIAG_REAL=r8_kind -I$(MODDIR)

# Run the test program.
TESTS = test_tridiagonal.sh

TEST_EXTENSIONS = .sh
SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \
$(abs_top_srcdir)/test_fms/tap-driver.sh

# These files will be included in the distribution.
EXTRA_DIST = test_tridiagonal.sh

# Clean up
CLEANFILES = *.nml *.out* *.dpi *.spi *.dyn *.spl *.o test_tridiagonal4 test_tridiagonal8 test_tridiagonal
173 changes: 173 additions & 0 deletions test_fms/tridiagonal/test_tridiagonal.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
!***********************************************************************
!* GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
#ifndef TEST_TRIDIAG_KIND
#define TEST_TRIDIAG_KIND 8
#endif

!> Tests the tridiagonal module routines (tri_invert and close_tridiagonal)
!! Tests reals with the kind value set above,
Copy link
Contributor

Choose a reason for hiding this comment

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

Perhaps for the future TODO, FMS tridiagonal solver could be compared to the LAPACK version to check that it's implemented correctly

program test_tridiagonal

use tridiagonal_mod
use platform_mod
use mpp_mod
use fms_mod

implicit none

integer, parameter :: IN_LEN = 8 !< length of input arrays
integer, parameter :: kindl = TEST_TRIDIAG_KIND !< kind value for all reals in this test
!! set by TEST_TRIDIAG_KIND cpp macro
real(TEST_TRIDIAG_KIND), allocatable :: d(:,:,:), x(:,:,:), ref_array(:,:,:)
real(TEST_TRIDIAG_KIND), allocatable :: a(:,:,:), b(:,:,:), c(:,:,:)
real(r4_kind), allocatable :: d_r4(:,:,:), x_r4(:,:,:)
real(r4_kind), allocatable :: a_r4(:,:,:), b_r4(:,:,:), c_r4(:,:,:)
real(r8_kind), allocatable :: d_r8(:,:,:), x_r8(:,:,:)
real(r8_kind), allocatable :: a_r8(:,:,:), b_r8(:,:,:), c_r8(:,:,:)
integer :: i, end, ierr, io
real(TEST_TRIDIAG_KIND) :: k
! nml
logical :: do_error_check = .false.
namelist / test_tridiagonal_nml/ do_error_check

call mpp_init

read (input_nml_file, test_tridiagonal_nml, iostat=io)
ierr = check_nml_error (io, 'test_tridiagonal_nml')

! allocate input and output arrays
allocate(d(1,1,IN_LEN))
allocate(a(1,1,IN_LEN))
allocate(b(1,1,IN_LEN))
allocate(c(1,1,IN_LEN))
allocate(x(1,1,IN_LEN))

!! simple test with only 1 coeff
a = 0.0_kindl
b = 1.0_kindl
c = 0.0_kindl
d = 5.0_kindl
call tri_invert(x, d, a, b, c)
if(any(x .ne. 5.0_kindl)) call mpp_error(FATAL, "test_tridiagonal: invalid results for 1 coefficient check")
!! check with stored data arrays
d = -5.0_kindl
call tri_invert(x, d)
if(any(x .ne. -5.0_kindl)) call mpp_error(FATAL, "test_tridiagonal: invalid results for 1 coefficient check")

! test with a,b,c
! 0.5x(n-2) + x(n-1) + 0.5x(n) = 1
!
! x(n) = k * [4, 1, 3, 2, 2, 3, 1, 4]
! k * [8 , 1, 7, 2, 6, .. ] = k *(-n/2 + ((n%2)*arr_length/2))
a = 0.5_kindl
b = 1.0_kindl
c = 0.5_kindl
d = 1.0_kindl
call tri_invert(x, d, a, b, c)
! set up reference answers
k = 1.0/(IN_LEN+1.0) * 2.0
rem1776 marked this conversation as resolved.
Show resolved Hide resolved
allocate(ref_array(1,1,IN_LEN))
do i=1, IN_LEN/2
end=IN_LEN-i+1
if(mod(i, 2) .eq. 1) then
ref_array(1,1,i) = -(i/2) + (mod(i,2)* IN_LEN/2)
Copy link
Contributor

Choose a reason for hiding this comment

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

real()'s needed

Copy link
Contributor Author

Choose a reason for hiding this comment

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

i'll add one to the whole thing but the math here (divisions) needs to done with int's so that the value * k is accurate.

ref_array(1,1,end) = -(i/2) + (mod(i,2)* IN_LEN/2)
else
ref_array(1,1,i) = i/2
ref_array(1,1,end) = i/2
endif
enddo
ref_array = ref_array * k
! check
do i=1, IN_LEN
if(ABS(x(1,1,i) - ref_array(1,1,i)) .gt. 0.1e-12_kindl) then
Copy link
Contributor

Choose a reason for hiding this comment

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

Is there a reason for this mismatch, answers not agreeing?

Copy link
Contributor

Choose a reason for hiding this comment

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

for i=1, is the answer supposed to be 8/9? It seems like both x and ref_array are off then hence the discrepancy

print *, i, x(1,1,i), ref_array(1,1,i)
call mpp_error(FATAL, "test_tridiagonal: failed reference check for tri_invert")
endif
enddo
!! check with stored data arrays
d = -1.0_kindl
ref_array = ref_array * -1.0_kindl
call tri_invert(x, d)
do i=1, IN_LEN
if(ABS(x(1,1,i) - ref_array(1,1,i)) .gt. 0.1e-12_kindl) then
print *, i, x(1,1,i), ref_array(1,1,i)
call mpp_error(FATAL, "test_tridiagonal: failed reference check for tri_invert with saved values")
endif
enddo
call close_tridiagonal()

!! tests for module state across kinds
!! default keeps stored values separate depending on kind
!! store_both_kinds argument can be specified to store both r4 and r8 kinds
if(kindl .eq. r8_kind) then
allocate(a_r4(1,1,IN_LEN), b_r4(1,1,IN_LEN), c_r4(1,1,IN_LEN))
allocate(d_r4(1,1,IN_LEN), x_r4(1,1,IN_LEN))
a_r4 = 0.0_r4_kind; b_r4 = 1.0_r4_kind; c_r4 = 0.0_r4_kind
d_r4 = 5.0_r4_kind; x_r4 = 0.0_r4_kind
a = 0.0_kindl; b = 2.0_kindl; c = 0.0_kindl
d = 5.0_kindl
! default, module variables distinct per kind
call tri_invert(x_r4, d_r4, a_r4, b_r4, c_r4)
! conditionally errors here for calling with unallocated a/b/c for kind
if( do_error_check ) call tri_invert(x, d)
call tri_invert(x, d, a, b, c)
! check both values are correct from prior state
call tri_invert(x_r4, d_r4)
call tri_invert(x, d)
if(any(x_r4 .ne. 5.0_r4_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r4 kind result")
if(any(x .ne. 2.5_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
call close_tridiagonal()
! run with storing for both kinds
call tri_invert(x_r4, d_r4, a_r4, b_r4, c_r4, store_both_kinds=.true.)
call tri_invert(x_r4, d_r4)
call tri_invert(x, d)
if(any(x_r4 .ne. 5.0_r4_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r4 kind result")
if(any(x .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
else
allocate(a_r8(1,1,IN_LEN), b_r8(1,1,IN_LEN), c_r8(1,1,IN_LEN))
allocate(d_r8(1,1,IN_LEN), x_r8(1,1,IN_LEN))
a_r8 = 0.0_r8_kind; b_r8 = 1.0_r8_kind; c_r8 = 0.0_r8_kind
d_r8 = 5.0_r8_kind; x_r8 = 0.0_r8_kind
a = 0.0_kindl; b = 2.0_kindl; c = 0.0_kindl
d = 5.0_kindl
! default, module variables distinct per kind
call tri_invert(x_r8, d_r8, a_r8, b_r8, c_r8)
! conditionally errors here for calling with unallocated a/b/c for kind
if( do_error_check ) call tri_invert(x, d)
call tri_invert(x, d, a, b, c)
! check both values are correct from prior state
call tri_invert(x_r8, d_r8)
call tri_invert(x, d)
if(any(x_r8 .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
if(any(x .ne. 2.5_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
call close_tridiagonal()
! run with storing for both kinds
call tri_invert(x_r8, d_r8, a_r8, b_r8, c_r8, store_both_kinds=.true.)
call tri_invert(x_r8, d_r8)
call tri_invert(x, d)
if(any(x_r8 .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
if(any(x .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
endif

call close_tridiagonal()

call mpp_exit

end program
51 changes: 51 additions & 0 deletions test_fms/tridiagonal/test_tridiagonal.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
#!/bin/sh

#***********************************************************************
#* GNU Lesser General Public License
#*
#* This file is part of the GFDL Flexible Modeling System (FMS).
#*
#* FMS is free software: you can redistribute it and/or modify it under
#* the terms of the GNU Lesser General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or (at
#* your option) any later version.
#*
#* FMS is distributed in the hope that it will be useful, but WITHOUT
#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
#* for more details.
#*
#* You should have received a copy of the GNU Lesser General Public
#* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
#***********************************************************************

# This is part of the GFDL FMS package. This is a shell script to
# execute tests in the test_fms/time_manager directory.

# Ryan Mulhall 9/2023

# Set common test settings.
. ../test-lib.sh

rm -f input.nml && touch input.nml

test_expect_success "test tridiagonal functionality 32 bit reals" '
mpirun -n 1 ./test_tridiagonal_r4
'
test_expect_success "test tridiagonal functionality 64 bit reals" '
mpirun -n 1 ./test_tridiagonal_r8
'
# tries to call without a,b,c args provided or previously set
cat <<_EOF > input.nml
&test_tridiagonal_nml
do_error_check = .true.
/
_EOF
test_expect_failure "error out if passed in incorrect real size (r4_kind)" '
mpirun -n 1 ./test_tridiagonal_r4
'
test_expect_failure "error out if passed in incorrect real size (r8_kind)" '
mpirun -n 1 ./test_tridiagonal_r8
'

test_done
2 changes: 1 addition & 1 deletion tridiagonal/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
# Ed Hartnett 2/22/19

# Include .h and .mod files.
AM_CPPFLAGS = -I$(top_srcdir)/include
AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/tridiagonal/include
AM_FCFLAGS = $(FC_MODINC). $(FC_MODOUT)$(MODDIR)

# Build this uninstalled convenience library.
Expand Down
Loading