mixed precision tridiagonal (#1363)

CMakeLists.txt

@@ -321,6 +321,7 @@ foreach(kind ${kinds})
                                                  field_manager/include
                                                  time_interp/include
                                                  tracer_manager/include
+                                                 tridiagonal/include
                                                  interpolator/include
                                                  coupler/include
                                                  data_override/include)

Makefile.am

@@ -35,8 +35,8 @@ endif
 # Make targets will be run in each subdirectory. Order is significant.
 SUBDIRS = \
   platform \
-  tridiagonal \
   mpp \
+  tridiagonal \
   constants \
   constants4 \
   memutils \

configure.ac

@@ -503,13 +503,13 @@ 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/topography/Makefile
   test_fms/column_diagnostics/Makefile
-
   FMS.pc
   ])
 

test_fms/Makefile.am

@@ -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 topography \
 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

test_fms/tridiagonal/Makefile.am

@@ -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

test_fms/tridiagonal/test_tridiagonal.F90

@@ -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,
+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_kindl/(IN_LEN+1.0_kindl) * 2.0_kindl
+    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) = real(-(i/2) + (mod(i,2)* IN_LEN/2), kindl)
+        ref_array(1,1,end) = real(-(i/2) + (mod(i,2)* IN_LEN/2), kindl)
+      else
+        ref_array(1,1,i) = real(i/2, kindl)
+        ref_array(1,1,end) = real(i/2, kindl)
+      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
+        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

test_fms/tridiagonal/test_tridiagonal.sh

@@ -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

tridiagonal/Makefile.am

@@ -23,14 +23,17 @@
 # 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.
 noinst_LTLIBRARIES = libtridiagonal.la
 
 # The convenience library depends on its source.
-libtridiagonal_la_SOURCES = tridiagonal.F90
+libtridiagonal_la_SOURCES = tridiagonal.F90           \
+                            include/tridiagonal.inc   \
+                            include/tridiagonal_r4.fh \
+                            include/tridiagonal_r8.fh
 
 # Mod file depends on its o file, is built and then installed.
 tridiagonal.lo: tridiagonal_mod.$(FC_MODEXT)