Feature/float128 a

.github/workflows/build-ci-containers-e4s.yml

@@ -40,7 +40,7 @@ jobs:
     strategy:
       matrix:
         indexsize: [32, 64]
-        precision: ['single', 'double', 'extended']
+        precision: ['single', 'double', 'extended', 'float128']
     runs-on: ubuntu-latest
     needs: build_e4s_base
     steps:

.github/workflows/ubuntu-latest.yml

@@ -26,18 +26,22 @@ jobs:
       max-parallel: 2
       matrix:
         indexsize: [32, 64]
-        precision: ['single', 'double', 'extended']
+        precision: ['single', 'double', 'extended', 'float128']
         buildtype: ['Debug', 'Release', 'RelWithDebInfo']
         tpls: ['ON']
         exclude:
           - buildtype: Debug
             precision: single
           - buildtype: Debug
             precision: extended
+          - buildtype: Debug
+            precision: float128
           - buildtype: Release
             precision: single
           - buildtype: Release
             precision: extended
+          - buildtype: Release
+            precision: float128
           - buildtype: RelWithDebInfo
             precision: double
 

.gitignore

@@ -3,6 +3,7 @@
 /.project
 /builddir
 /build*
+/cmake-build*
 /instdir
 /install*
 /TPLs
@@ -18,9 +19,12 @@
 *.directory
 /.pydevproject
 .vscode
+.vs
+/out
 compile_commands.json
 .venv
 .clangd
+.idea
 
 # custom test environment setup script
 /test/env/env.sh

benchmarks/diffusion_2D/diffusion_2D.cpp

@@ -691,23 +691,20 @@ int UserOutput::open(UserData* udata)
   if (outproc)
   {
     cout << scientific;
-    cout << setprecision(numeric_limits<sunrealtype>::digits10);
+    cout << setprecision(SUN_DIGITS10);
     cout << endl;
     if (error)
     {
-      cout << "          t           ";
-      cout << "          ||u||_rms      ";
-      cout << "          max error      " << endl;
-      cout << " ---------------------";
-      cout << "-------------------------";
-      cout << "-------------------------" << endl;
+      cout << setw(SUN_DIGITS10/2+3) <<" "<<"t"<< setw(SUN_DIGITS10/2+4)  <<" ";
+      cout << setw(SUN_DIGITS10/2+1) <<" "<<"||u||_rms" << setw(SUN_DIGITS10/2+1) <<" ";
+      cout << setw(SUN_DIGITS10/2+1) <<" "<<"max error" << setw(SUN_DIGITS10/2+1) <<" "<< endl;
+      cout << setfill('_') << setw(SUN_DIGITS10*3+28) << setfill(' ')<<endl;
     }
     else
     {
-      cout << "          t           ";
-      cout << "          ||u||_rms      " << endl;
-      cout << " ---------------------";
-      cout << "-------------------------" << endl;
+      cout << setw(SUN_DIGITS10/2+3) <<" "<<"t"<< setw(SUN_DIGITS10/2+4)  <<" ";
+      cout << setw(SUN_DIGITS10/2+1) <<" "<<"||u||_rms" << setw(SUN_DIGITS10/2+1) <<" " << endl;
+      cout << setfill('_') << setw(SUN_DIGITS10*2+18) << " "<<endl;
     }
   }
 
@@ -739,7 +736,7 @@ int UserOutput::open(UserData* udata)
     uoutstream << "# je  " << udata->je << endl;
 
     uoutstream << scientific;
-    uoutstream << setprecision(numeric_limits<sunrealtype>::digits10);
+    uoutstream << setprecision(SUN_DIGITS10);
 
     if (error)
     {
@@ -768,7 +765,7 @@ int UserOutput::open(UserData* udata)
       eoutstream << "# je  " << udata->je << endl;
 
       eoutstream << scientific;
-      eoutstream << setprecision(numeric_limits<sunrealtype>::digits10);
+      eoutstream << setprecision(SUN_DIGITS10);
     }
   }
 
@@ -794,16 +791,16 @@ int UserOutput::write(sunrealtype t, N_Vector u, UserData* udata)
     }
 
     // Compute rms norm of the state
-    sunrealtype urms = sqrt(N_VDotProd(u, u) / udata->nx / udata->ny);
+    sunrealtype urms = SUNRsqrt(N_VDotProd(u, u) / udata->nx / udata->ny);
 
     // Output current status
     if (outproc)
     {
       if (error)
       {
-        cout << setw(22) << t << setw(25) << urms << setw(25) << max << endl;
+        cout << setw(SUN_DIGITS10+7) << t << setw(SUN_DIGITS10+10) << urms << setw(SUN_DIGITS10+10) << max << endl;
       }
-      else { cout << setw(22) << t << setw(25) << urms << endl; }
+      else { cout << setw(SUN_DIGITS10+7) << t << setw(SUN_DIGITS10+10) << urms << endl; }
     }
 
     // Write solution and error to disk

benchmarks/diffusion_2D/diffusion_2D.hpp

@@ -29,6 +29,8 @@
 #include <limits>
 #include <sstream>
 #include <sundials/sundials_core.hpp>
+#include <sundials/sundials_types.h>
+#include <sundials/sundials_math.h>
 #include <vector>
 
 #include "mpi.h"

benchmarks/diffusion_2D/mpi_serial/CMakeLists.txt

@@ -55,7 +55,7 @@ foreach(test_tuple ${tests})
   target_include_directories(${target} PRIVATE ${benchmark_prefix})
 
   target_link_libraries(${target} PRIVATE sundials_${package}
-                                          sundials_nvecparallel MPI::MPI_CXX)
+                                          sundials_nvecparallel MPI::MPI_CXX ${SUNDIALS_MATH_LIBRARY})
 
   if(BUILD_SUNLINSOL_SUPERLUDIST)
     target_compile_definitions(${target} PRIVATE USE_SUPERLU_DIST)

benchmarks/diffusion_2D/mpi_serial/diffusion.cpp

@@ -66,8 +66,8 @@ int laplacian(sunrealtype t, N_Vector u, N_Vector f, UserData* udata)
     sunrealtype bx = (udata->kx) * TWO * PI * PI;
     sunrealtype by = (udata->ky) * TWO * PI * PI;
 
-    sunrealtype sin_t_cos_t = sin(PI * t) * cos(PI * t);
-    sunrealtype cos_sqr_t   = cos(PI * t) * cos(PI * t);
+    sunrealtype sin_t_cos_t = SUNRsin(PI * t) * SUNRcos(PI * t);
+    sunrealtype cos_sqr_t   = SUNRcos(PI * t) * SUNRcos(PI * t);
 
     for (j = jstart; j < jend; j++)
     {
@@ -76,11 +76,11 @@ int laplacian(sunrealtype t, N_Vector u, N_Vector f, UserData* udata)
         x = (udata->is + i) * udata->dx;
         y = (udata->js + j) * udata->dy;
 
-        sin_sqr_x = sin(PI * x) * sin(PI * x);
-        sin_sqr_y = sin(PI * y) * sin(PI * y);
+        sin_sqr_x = SUNRsin(PI * x) * SUNRsin(PI * x);
+        sin_sqr_y = SUNRsin(PI * y) * SUNRsin(PI * y);
 
-        cos_sqr_x = cos(PI * x) * cos(PI * x);
-        cos_sqr_y = cos(PI * y) * cos(PI * y);
+        cos_sqr_x = SUNRcos(PI * x) * SUNRcos(PI * x);
+        cos_sqr_y = SUNRcos(PI * y) * SUNRcos(PI * y);
 
         farray[IDX(i, j, nx_loc)] =
           -TWO * PI * sin_sqr_x * sin_sqr_y * sin_t_cos_t -

benchmarks/diffusion_2D/mpi_serial/solution.cpp

@@ -29,7 +29,7 @@ int Solution(sunrealtype t, N_Vector u, UserData* udata)
   sunrealtype sin_sqr_x, sin_sqr_y;
 
   // Constants for computing solution
-  sunrealtype cos_sqr_t = cos(PI * t) * cos(PI * t);
+  sunrealtype cos_sqr_t = SUNRcos(PI * t) * SUNRcos(PI * t);
 
   // Iterative over domain interior
   sunindextype istart = (udata->HaveNbrW) ? 0 : 1;
@@ -48,8 +48,8 @@ int Solution(sunrealtype t, N_Vector u, UserData* udata)
       x = (udata->is + i) * udata->dx;
       y = (udata->js + j) * udata->dy;
 
-      sin_sqr_x = sin(PI * x) * sin(PI * x);
-      sin_sqr_y = sin(PI * y) * sin(PI * y);
+      sin_sqr_x = SUNRsin(PI * x) * SUNRsin(PI * x);
+      sin_sqr_y = SUNRsin(PI * y) * SUNRsin(PI * y);
 
       uarray[IDX(i, j, udata->nx_loc)] = sin_sqr_x * sin_sqr_y * cos_sqr_t;
     }
@@ -68,7 +68,7 @@ int SolutionDerivative(sunrealtype t, N_Vector up, UserData* udata)
   sunrealtype sin_sqr_x, sin_sqr_y;
 
   // Constants for computing solution
-  sunrealtype cos_sin_t = -TWO * PI * cos(PI * t) * sin(PI * t);
+  sunrealtype cos_sin_t = -TWO * PI * SUNRcos(PI * t) * SUNRsin(PI * t);
 
   // Iterative over domain interior
   sunindextype istart = (udata->HaveNbrW) ? 0 : 1;
@@ -87,8 +87,8 @@ int SolutionDerivative(sunrealtype t, N_Vector up, UserData* udata)
       x = (udata->is + i) * udata->dx;
       y = (udata->js + j) * udata->dy;
 
-      sin_sqr_x = sin(PI * x) * sin(PI * x);
-      sin_sqr_y = sin(PI * y) * sin(PI * y);
+      sin_sqr_x = SUNRsin(PI * x) * SUNRsin(PI * x);
+      sin_sqr_y = SUNRsin(PI * y) * SUNRsin(PI * y);
 
       uparray[IDX(i, j, udata->nx_loc)] = sin_sqr_x * sin_sqr_y * cos_sin_t;
     }

benchmarks/nvector/openmpdev/test_nvector_performance_openmpdev.c

@@ -36,8 +36,9 @@ static int FinalizeClearCache();
  * --------------------------------------------------------------------*/
 int main(int argc, char* argv[])
 {
-  N_Vector X;          /* test vector       */
-  sunindextype veclen; /* vector length     */
+  SUNContext ctx = NULL; /* SUNDIALS context */
+  N_Vector X     = NULL; /* test vector      */
+  sunindextype veclen;   /* vector length    */
 
   int print_timing; /* output timings    */
   int ntests;       /* number of tests   */
@@ -59,27 +60,27 @@ int main(int argc, char* argv[])
     return (-1);
   }
 
-  veclen = atol(argv[1]);
+  veclen = (sunindextype)atol(argv[1]);
   if (veclen <= 0)
   {
     printf("ERROR: length of vector must be a positive integer \n");
     return (-1);
   }
 
-  nvecs = atol(argv[2]);
+  nvecs = (int)atol(argv[2]);
   if (nvecs < 1) { printf("WARNING: Fused operation tests disabled\n"); }
 
-  nsums = atol(argv[3]);
+  nsums = (int)atol(argv[3]);
   if (nsums < 1) { printf("WARNING: Some fused operation tests disabled\n"); }
 
-  ntests = atol(argv[4]);
+  ntests = (int)atol(argv[4]);
   if (ntests <= 0)
   {
     printf("ERROR: number of tests must be a positive integer \n");
     return (-1);
   }
 
-  cachesize = atol(argv[5]);
+  cachesize = (int)atol(argv[5]);
   if (cachesize < 0)
   {
     printf("ERROR: cache size (MB) must be a non-negative integer \n");
@@ -98,7 +99,7 @@ int main(int argc, char* argv[])
   printf("  timing on/off         %d  \n", print_timing);
 
   /* Create vectors */
-  X = N_VNew_OpenMPDEV(veclen);
+  X = N_VNew_OpenMPDEV(veclen, ctx);
 
   /* run tests */
   if (print_timing) { printf("\n\n standard operations:\n"); }
@@ -154,6 +155,9 @@ int main(int argc, char* argv[])
 
   FinalizeClearCache();
 
+  flag = SUNContext_Free(&ctx);
+  if (flag) { return flag; }
+
   printf("\nFinished Tests\n");
 
   return (flag);

benchmarks/nvector/test_nvector_performance.c

@@ -2860,6 +2860,6 @@ static void time_stats(N_Vector X, double* times, int num_warmups, int ntests,
     {
       *sdev += (times[i] - *avg) * (times[i] - *avg);
     }
-    *sdev = sqrt(*sdev / (ntests - 1));
+    *sdev = SUNRsqrt(*sdev / (ntests - 1));
   }
 }

bindings/sundials4py/sundials/sundials_core.cpp

@@ -89,6 +89,8 @@ void bind_core(nb::module_& m)
   m.attr("sunrealtype") = np.attr("float64");
 #elif defined(SUNDIALS_EXTENDED_PRECISION)
   m.attr("sunrealtype") = np.attr("longdouble");
+#elif defined(SUNDIALS_FLOAT128_PRECISION)
+  m.attr("sunrealtype") = np.attr("float128");
 #else
 #error Unknown sunrealtype, email sundials-users@llnl.gov
 #endif

cmake/SundialsBuildOptionsPre.cmake

@@ -33,7 +33,7 @@ endif()
 # Option to specify precision (sunrealtype)
 # ---------------------------------------------------------------
 
-set(DOCSTR "single, double, or extended")
+set(DOCSTR "single, double, extended or float128")
 sundials_option(SUNDIALS_PRECISION STRING "${DOCSTR}" "DOUBLE")
 string(TOUPPER ${SUNDIALS_PRECISION} _upper_SUNDIALS_PRECISION)
 set(SUNDIALS_PRECISION
@@ -133,7 +133,7 @@ endif()
 
 if(UNIX)
   sundials_option(SUNDIALS_MATH_LIBRARY PATH
-                  "Which math library (e.g., libm) to link to" "-lm" ADVANCED)
+                  "Which math library (e.g., libm) to link to" "-lm -lquadmath" ADVANCED)
 else()
   sundials_option(SUNDIALS_MATH_LIBRARY PATH
                   "Which math library (e.g., libm) to link to" "" ADVANCED)

cmake/SundialsSetupCompilers.cmake

@@ -94,7 +94,7 @@ if(ENABLE_ALL_WARNINGS)
   # CI tests
 
   # Avoid numerous warnings from printf
-  if(SUNDIALS_PRECISION MATCHES "EXTENDED")
+  if(SUNDIALS_PRECISION MATCHES "EXTENDED|FLOAT128")
     set(WARNING_FLAGS "-Wdouble-promotion ${WARNING_FLAGS}")
   endif()
 

cmake/SundialsSetupCuda.cmake

@@ -57,12 +57,12 @@ if((CMAKE_CXX_COMPILER_ID MATCHES GNU)
 endif()
 
 # ===============================================================
-# Prohibit CUDA interface when using extended precision.
+# Prohibit CUDA interface when using extended or float128 precision.
 # ===============================================================
 
-if(SUNDIALS_PRECISION MATCHES "EXTENDED")
+if(SUNDIALS_PRECISION MATCHES "EXTENDED|FLOAT128")
   message(
-    FATAL_ERROR "CUDA interfaces are incompatible with extended precision.")
+    FATAL_ERROR "CUDA interfaces are incompatible with ${SUNDIALS_PRECISION} precision.")
 endif()
 
 # ===============================================================

cmake/macros/SundialsAddNvectorBenchmark.cmake

@@ -36,7 +36,7 @@ function(sundials_add_nvector_benchmark NAME)
   target_include_directories(${NAME} PRIVATE ${BENCHMARKS_DIR}/nvector)
 
   target_link_libraries(${NAME} PRIVATE ${arg_SUNDIALS_TARGETS}
-                                        ${arg_LINK_LIBRARIES} -lm)
+                                        ${arg_LINK_LIBRARIES} ${SUNDIALS_MATH_LIBRARY})
 
   install(TARGETS ${NAME}
           DESTINATION "${BENCHMARKS_INSTALL_PATH}/${arg_INSTALL_SUBDIR}")

cmake/tpl/SundialsGinkgo.cmake

@@ -100,10 +100,10 @@ endif()
 # -----------------------------------------------------------------------------
 
 if(Ginkgo_FOUND AND (NOT GINKGO_WORKS))
-  if(SUNDIALS_PRECISION MATCHES "extended|EXTENDED")
+  if(SUNDIALS_PRECISION MATCHES "extended|EXTENDED|float128|FLOAT128")
     message(
       FATAL_ERROR
-        "SUNDIALS GINKGO interface is not compatible with extended precision")
+        "SUNDIALS GINKGO interface is not compatible with extended or float128 precision")
   endif()
 
   if(SUNDIALS_GINKGO_BACKENDS MATCHES "CUDA" AND NOT ENABLE_CUDA)

cmake/tpl/SundialsKLU.cmake

@@ -39,8 +39,8 @@ endif()
 # Section 2: Check to make sure options are compatible
 # -----------------------------------------------------------------------------
 
-# KLU does not support single or extended precision
-if(SUNDIALS_PRECISION MATCHES "SINGLE" OR SUNDIALS_PRECISION MATCHES "EXTENDED")
+# KLU does not support single, extended or float128 precision
+if(SUNDIALS_PRECISION MATCHES "SINGLE|EXTENDED|FLOAT128")
   message(
     FATAL_ERROR "KLU is not compatible with ${SUNDIALS_PRECISION} precision")
 endif()

cmake/tpl/SundialsLapack.cmake

@@ -35,8 +35,8 @@ include_guard(GLOBAL)
 # Section 2: Check to make sure options are compatible
 # -----------------------------------------------------------------------------
 
-# LAPACK does not support extended precision
-if(ENABLE_LAPACK AND SUNDIALS_PRECISION MATCHES "EXTENDED")
+# LAPACK does not support extended or float128 precision
+if(ENABLE_LAPACK AND SUNDIALS_PRECISION MATCHES "EXTENDED|FLOAT128")
   message(
     FATAL_ERROR "LAPACK is not compatible with ${SUNDIALS_PRECISION} precision")
 endif()

cmake/tpl/SundialsMAGMA.cmake

@@ -39,10 +39,10 @@ endif()
 # Section 2: Check to make sure options are compatible
 # -----------------------------------------------------------------------------
 
-if(SUNDIALS_PRECISION MATCHES "extended")
+if(SUNDIALS_PRECISION MATCHES "extended|float128")
   message(
     FATAL_ERROR
-      "SUNDIALS MAGMA interface is not compatible with extended precision")
+      "SUNDIALS MAGMA interface is not compatible with ${SUNDIALS_PRECISION} precision")
 endif()
 
 # -----------------------------------------------------------------------------