diff --git a/src/support/CMakeLists.txt b/src/support/CMakeLists.txt
index e78fc1642a7edc83201d189487d6a77c45ad3331..44f60aa117b4a1a62a129b14f469f8bdd0020d65 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 e6d9d2159a32f9e5060e4f85a32569e48696757c..30c82bd2e98521f99b09abf9343ee1a5b52f6185 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 03eb9775ac508a56cc81639bcc3e95cabc2efae8..5136c6abbb9e89cc07825ad9b8e053fd7561da43 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 bb8e6df7c70268c727eafb446480ce0d6cdaa7a3..56e3d250530be9eb668b87c86124d5db9f2d781b 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 6a60f2ce3988c7d759f09296f28d2c1317a61f6c..5859b9d3b9b3875fead0a48559065de0c74d8027 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 4ee5dc910e2f3bb9b2e69568aa3d78b5d3060c77..a3f3ba1a603febc37c5551e2ff8e70b61481e2ff 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 0000000000000000000000000000000000000000..664fc1e7d3be31e981bc959b65b9d44743175a3e
--- /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 0000000000000000000000000000000000000000..df452e4e92c10931497c4e87e3efafd778182e0c
--- /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 a33673c960b33d558a713a4739db42e6ba267cca..0eb5a8dc3ed9340ab34eab49da7c84d1b5d97207 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 8f120ef5735161d980f6a6d02d80fc7f171b7fa0..4e09ff335368b2cbc532fd2ade2aee8a8259415d 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