Added a new argument to the functions of `mo_lib_loopindices.cpp` to fix a bug...

Added a new argument to the functions of `mo_lib_loopindices.cpp` to fix a bug regarding startindex (!32)

This made the code to produce bit-identical results for both Fortran and
C++

Co-authored-by: Pradipta Samanta <samanta@dkrz.de>
Merged-by: Pradipta Samanta <samanta@dkrz.de>
Changelog: bugfix

src/support/CMakeLists.txt

@@ -19,7 +19,8 @@ add_library(
   mo_math_types.f90
   mo_math_utilities.cpp
   mo_math_utilities.F90
-  mo_random_number_generators.F90)
+  mo_random_number_generators.F90
+  support_bindings.cpp)
 
 add_library(${PROJECT_NAME}::support ALIAS iconmath-support)
 

src/support/mo_lib_loopindices.cpp

@@ -11,47 +11,75 @@
 
 #include <algorithm> // For std::max
 
-extern "C" {
-    // get_indices_c_lib function
-    void get_indices_c_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out) {
-        if (i_blk == i_startblk) {
-            i_startidx_out = std::max(1, i_startidx_in);
-            i_endidx_out = nproma;
-            if (i_blk == i_endblk) {
-                i_endidx_out = i_endidx_in;
-            }
-        } else if (i_blk == i_endblk) {
-            i_startidx_out = 1;
+// get_indices_c_lib function
+void get_indices_c_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+                        const int i_blk, const int i_startblk, const int i_endblk,
+                        int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
+    
+    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from
+    //(either fortran or c++), the first index should be either 0 or 1.
+    int first_index;
+    if (called_from_cpp)
+        first_index = 0;
+    else
+        first_index = 1;                   
+    
+    if (i_blk == i_startblk) {
+        i_startidx_out = std::max(first_index, i_startidx_in);
+        i_endidx_out = nproma;
+        if (i_blk == i_endblk) {
             i_endidx_out = i_endidx_in;
-        } else {
-            i_startidx_out = 1;
-            i_endidx_out = nproma;
         }
+    } else if (i_blk == i_endblk) {
+        i_startidx_out = first_index;
+        i_endidx_out = i_endidx_in;
+    } else {
+        i_startidx_out = first_index;
+        i_endidx_out = nproma;
     }
+}
 
-    // get_indices_e_lib function
-    void get_indices_e_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out) {
-        i_startidx_out = (i_blk != i_startblk) ? 1 : std::max(1, i_startidx_in);
-        i_endidx_out = (i_blk != i_endblk) ? nproma : i_endidx_in;
-    }
+// get_indices_e_lib function
+void get_indices_e_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+                        const int i_blk, const int i_startblk, const int i_endblk,
+                        int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
+    
+    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from, 
+    //the first index should be either 0 or 1.
+    int first_index;
+    if (called_from_cpp)
+        first_index = 0;
+    else
+        first_index = 1;
+
+    i_startidx_out = (i_blk != i_startblk) ? first_index : std::max(first_index, i_startidx_in);
+    i_endidx_out = (i_blk != i_endblk) ? nproma : i_endidx_in;
+}
+
+// get_indices_v_lib function
+void get_indices_v_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+                        const int i_blk, const int i_startblk, const int i_endblk,
+                        int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true) {
+    
+    //Since code is ported incrementally from Fortran to C++, depending on where the function is called from, 
+    //the first index should be either 0 or 1.
+    int first_index;
+    if (called_from_cpp)
+        first_index = 0;
+    else
+        first_index = 1;
 
-    // get_indices_v_lib function
-    void get_indices_v_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out) {
-        if (i_blk == i_startblk) {
-            i_startidx_out = i_startidx_in;
-            i_endidx_out = nproma;
-            if (i_blk == i_endblk) {
-                i_endidx_out = i_endidx_in;
-            }
-        } else if (i_blk == i_endblk) {
-            i_startidx_out = 1;
+    if (i_blk == i_startblk) {
+        i_startidx_out = i_startidx_in;
+        i_endidx_out = nproma;
+        if (i_blk == i_endblk) {
             i_endidx_out = i_endidx_in;
-        } else {
-            i_startidx_out = 1;
-            i_endidx_out = nproma;
         }
+    } else if (i_blk == i_endblk) {
+        i_startidx_out = first_index;
+        i_endidx_out = i_endidx_in;
+    } else {
+        i_startidx_out = first_index;
+        i_endidx_out = nproma;
     }
-}
+}
\ No newline at end of file

src/support/mo_lib_loopindices.hpp

@@ -8,15 +8,17 @@
 // See LICENSES/ for license information
 // SPDX-License-Identifier: BSD-3-Clause
 // ---------------------------------------------------------------
+#pragma once
 
-extern "C" {
-    // get_indices_c_lib function
-    void get_indices_c_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out);
+void get_indices_c_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+    const int i_blk, const int i_startblk, const int i_endblk,
+    int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true);
 
-    void get_indices_e_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out);
+void get_indices_e_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+                        const int i_blk, const int i_startblk, const int i_endblk,
+                        int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true);
 
-    void get_indices_v_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
-                           int &i_startidx_out, int &i_endidx_out);
-}
+void get_indices_v_lib(const int i_startidx_in, const int i_endidx_in, const int nproma, 
+                            const int i_blk, const int i_startblk, const int i_endblk,
+                            int &i_startidx_out, int &i_endidx_out, const bool called_from_cpp=true);
+    
\ No newline at end of file

src/support/mo_math_utilities.F90

@@ -82,7 +82,7 @@ MODULE mo_math_utilities
 #ifndef __USE_CPP_BINDINGS
   PUBLIC :: tdma_solver_vec
 #else
-  PUBLIC :: tdma_solver_vec_double
+  PUBLIC :: tdma_solver_vec_dp
 #endif
   PUBLIC :: check_orientation
 
@@ -251,14 +251,14 @@ CONTAINS
 
   ! C++ binding for tdma_solver_vec
   INTERFACE
-    SUBROUTINE tdma_solver_vec_double(a, b, c, d, slev, elev, startidx, endidx, nrows, ncols, varout, opt_acc_queue) &
-      BIND(C, NAME="tdma_solver_vec_double")
+    SUBROUTINE tdma_solver_vec_dp(a, b, c, d, slev, elev, startidx, endidx, nrows, ncols, varout, opt_acc_queue) &
+      BIND(C, NAME="tdma_solver_vec_dp")
       IMPORT :: c_double, c_int
       REAL(c_double), INTENT(IN) :: a(*), b(*), c(*), d(*)
       INTEGER(c_int), VALUE :: slev, elev, startidx, endidx, nrows, ncols
       REAL(c_double), INTENT(OUT) :: varout(*)
       INTEGER(c_int), OPTIONAL :: opt_acc_queue
-    END SUBROUTINE tdma_solver_vec_double
+    END SUBROUTINE tdma_solver_vec_dp
   END INTERFACE
 
 CONTAINS

src/support/mo_math_utilities.cpp

@@ -75,11 +75,13 @@ void tdma_solver_vec(const T* a, const T* b, const T* c, const T* d,
   std::cout << "Elapsed time for tdma_solver_vec (Kokkos): " << elapsed_time.count() << " seconds" << std::endl;
 }
 
-extern "C" {
+template
+void tdma_solver_vec<double>(const double* a, const double* b, const double* c, const double* d,
+  int slev, int elev, int startidx, int endidx,
+  int nrows, int ncols, double* varout);
+
+template
+void tdma_solver_vec<float>(const float* a, const float* b, const float* c, const float* d,
+  int slev, int elev, int startidx, int endidx,
+  int nrows, int ncols, float* varout);
 
-  void tdma_solver_vec_double(const double* a, const double* b, const double* c, const double* d,
-      int slev, int elev, int startidx, int endidx,
-      int nrows, int ncols, double* varout) {
-    tdma_solver_vec<double>(a, b, c, d, slev, elev, startidx, endidx, nrows, ncols, varout);
-  }
-}

src/support/mo_math_utilities.hpp

@@ -17,10 +17,4 @@
 template <typename T>
 void tdma_solver_vec(const T* a, const T* b, const T* c, const T* d,
                               int slev, int elev, int startidx, int endidx,
-                              int nrows, int ncols, T* varout);
-
-extern "C" {
-  void tdma_solver_vec_double(const double* a, const double* b, const double* c, const double* d,
-                              int slev, int elev, int startidx, int endidx,
-                              int nrows, int ncols, double* varout);
-}
+                              int nrows, int ncols, T* varout);
\ No newline at end of file

src/support/support_bindings.cpp

+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include "support_bindings.h"
+#include "mo_lib_loopindices.hpp"
+#include "mo_math_utilities.hpp"
+
+
+// mo_loop_indices.F90
+// C wrappers for C++ functionality
+void get_indices_c_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+    int &i_startidx_out, int &i_endidx_out){
+        get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk,
+                            i_endblk, i_startidx_out, i_endidx_out, false);
+}
+void get_indices_e_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+    int &i_startidx_out, int &i_endidx_out){
+        get_indices_e_lib(i_startidx_in, i_endidx_in, nproma,i_blk, i_startblk, i_endblk,
+                            i_startidx_out, i_endidx_out, false);
+}
+
+void get_indices_v_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+    int &i_startidx_out, int &i_endidx_out){
+        get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk,
+            i_startidx_out,i_endidx_out, false);
+}
+
+void tdma_solver_vec_dp(const double* a, const double* b, const double* c, const double* d,
+    int slev, int elev, int startidx, int endidx,
+    int nrows, int ncols, double* varout){
+    
+    tdma_solver_vec<double>(a, b, c, d, slev, elev, startidx, endidx, nrows, ncols, varout);
+
+}
+
+void tdma_solver_vec_sp(const float* a, const float* b, const float* c, const float* d,
+      int slev, int elev, int startidx, int endidx,
+      int nrows, int ncols, float* varout){
+
+    tdma_solver_vec<float>(a, b, c, d, slev, elev, startidx, endidx, nrows, ncols, varout);
+
+}

src/support/support_bindings.h

+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+#pragma once
+
+extern "C" {
+    // mo_loop_indices.F90
+    void get_indices_c_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+                            int &i_startidx_out, int &i_endidx_out);
+
+    void get_indices_e_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+                            int &i_startidx_out, int &i_endidx_out);
+
+    void get_indices_v_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
+                            int &i_startidx_out, int &i_endidx_out);
+
+    //mo_math_utilities.F90
+    void tdma_solver_vec_dp(const double* a, const double* b, const double* c, const double* d,
+        int slev, int elev, int startidx, int endidx,
+        int nrows, int ncols, double* varout);
+
+    void tdma_solver_vec_sp(const float* a, const float* b, const float* c, const float* d,
+          int slev, int elev, int startidx, int endidx,
+          int nrows, int ncols, float* varout);
+
+}

test/c/test_interpolation_vector.cpp

@@ -75,12 +75,6 @@ TEST(Edges2CellsTest, DPTest) {
 
   // Check that for each computed cell in p_u_out and p_v_out, the value is 6.
   // This is because for each cell, the kernel adds 6 terms of 1*1.
-  p_u_ref[0] = 0.0;
-  p_u_ref[8] = 0.0;
-  p_u_ref[10] = 0.0;
-  p_v_ref[0] = 0.0;
-  p_v_ref[8] = 0.0;
-  p_v_ref[10] = 0.0;
   for (size_t idx = 0; idx < p_u_out.size(); ++idx) {
     EXPECT_NEAR(p_u_out[idx], p_u_ref[idx], 1e-12);
     EXPECT_NEAR(p_v_out[idx], p_v_ref[idx], 1e-12);
@@ -120,12 +114,6 @@ TEST(Edges2CellsTest, SPTest) {
     nproma,
     nlev, nblks_e, nblks_c);
 
-  p_u_ref[0] = 0.0f;
-  p_u_ref[8] = 0.0f;
-  p_u_ref[10] = 0.0f;
-  p_v_ref[0] = 0.0f;
-  p_v_ref[8] = 0.0f;
-  p_v_ref[10] = 0.0f;
   // Verify that every computed output equals 6.
   for (size_t idx = 0; idx < p_u_out.size(); ++idx) {
     EXPECT_NEAR(p_u_out[idx], p_u_ref[idx], 1e-5f);

test/c/test_tdma_solver.cpp

@@ -51,7 +51,7 @@ protected:
 
 TEST_F(TDMASolverTestFixture, FullTest) {
   // Call the solver over the full range:
-  tdma_solver_vec_double(a.data(), b.data(), c.data(), d.data(),
+  tdma_solver_vec<double>(a.data(), b.data(), c.data(), d.data(),
                          0, n, 0, n, n, n, x.data());
 
   // Compute the sum of all elements in the output matrix.
@@ -71,7 +71,7 @@ TEST_F(TDMASolverTestFixture, FullTest) {
 TEST_F(TDMASolverTestFixture, PartialTest) {
   // Call the solver for a partial region:
   // For C++: slev = 1, elev = n-1, startidx = 1, endidx = n-1.
-  tdma_solver_vec_double(a.data(), b.data(), c.data(), d.data(),
+  tdma_solver_vec<double>(a.data(), b.data(), c.data(), d.data(),
                          1, n - 1, 1, n - 1, n, n, x.data());
 
   // Compute the sum over a region