From c358ed6e4e983fdf4e63d712aa32683de89ed995 Mon Sep 17 00:00:00 2001
From: Harshada Balasubramanian <harshada.balasubramanian@mpimet.mpg.de>
Date: Mon, 24 Feb 2025 21:30:54 +0000
Subject: [PATCH] Added a new argument to the functions of
`mo_lib_loopindices.cpp` to fix a bug regarding startindex
(icon-libraries/libiconmath!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 | 3 +-
src/support/mo_lib_loopindices.cpp | 100 +++++++++++++++++----------
src/support/mo_lib_loopindices.hpp | 20 +++---
src/support/mo_math_utilities.F90 | 8 +--
src/support/mo_math_utilities.cpp | 16 +++--
src/support/mo_math_utilities.hpp | 8 +--
src/support/support_bindings.cpp | 50 ++++++++++++++
src/support/support_bindings.h | 33 +++++++++
test/c/test_interpolation_vector.cpp | 12 ----
test/c/test_tdma_solver.cpp | 4 +-
10 files changed, 176 insertions(+), 78 deletions(-)
create mode 100644 src/support/support_bindings.cpp
create mode 100644 src/support/support_bindings.h
diff --git a/src/support/CMakeLists.txt b/src/support/CMakeLists.txt
index e78fc16..44f60aa 100644
--- a/src/support/CMakeLists.txt
+++ b/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)
diff --git a/src/support/mo_lib_loopindices.cpp b/src/support/mo_lib_loopindices.cpp
index e6d9d21..30c82bd 100644
--- a/src/support/mo_lib_loopindices.cpp
+++ b/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
diff --git a/src/support/mo_lib_loopindices.hpp b/src/support/mo_lib_loopindices.hpp
index 03eb977..5136c6a 100644
--- a/src/support/mo_lib_loopindices.hpp
+++ b/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
diff --git a/src/support/mo_math_utilities.F90 b/src/support/mo_math_utilities.F90
index bb8e6df..56e3d25 100644
--- a/src/support/mo_math_utilities.F90
+++ b/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
diff --git a/src/support/mo_math_utilities.cpp b/src/support/mo_math_utilities.cpp
index 6a60f2c..5859b9d 100644
--- a/src/support/mo_math_utilities.cpp
+++ b/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);
- }
-}
diff --git a/src/support/mo_math_utilities.hpp b/src/support/mo_math_utilities.hpp
index 4ee5dc9..a3f3ba1 100644
--- a/src/support/mo_math_utilities.hpp
+++ b/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
diff --git a/src/support/support_bindings.cpp b/src/support/support_bindings.cpp
new file mode 100644
index 0000000..664fc1e
--- /dev/null
+++ b/src/support/support_bindings.cpp
@@ -0,0 +1,50 @@
+// 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);
+
+}
diff --git a/src/support/support_bindings.h b/src/support/support_bindings.h
new file mode 100644
index 0000000..df452e4
--- /dev/null
+++ b/src/support/support_bindings.h
@@ -0,0 +1,33 @@
+// 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);
+
+}
diff --git a/test/c/test_interpolation_vector.cpp b/test/c/test_interpolation_vector.cpp
index a33673c..0eb5a8d 100644
--- a/test/c/test_interpolation_vector.cpp
+++ b/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);
diff --git a/test/c/test_tdma_solver.cpp b/test/c/test_tdma_solver.cpp
index 8f120ef..4e09ff3 100644
--- a/test/c/test_tdma_solver.cpp
+++ b/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
--
GitLab