Fixed abstolerance in cvodes

sn248 · Nov 11, 2024 · 19f76ab · 19f76ab
1 parent 9414bed
commit 19f76ab
Show file tree

Hide file tree

Showing 6 changed files with 48 additions and 42 deletions.
diff --git a/NEWS b/NEWS
@@ -1,6 +1,11 @@
+sundialr v0.1.5
+================
+Bug fixes:
+* Fixed the bug in assigning absolute tolerances
+
 sundialr v0.1.4.1
 ================
-Bug fixex:
+Bug fixes:
 * Fixed the linking problem due to multiple defined sybmols
 
 sundialr v0.1.4

diff --git a/R/zzz.R b/R/zzz.R
@@ -1,4 +1,4 @@
-.onAttach <- function(...){
-  packageStartupMessage(
-    "SUNDIALS - Copyright (c) 2002-2020, Lawrence Livermore National Security and Southern Methodist University.\nAll rights reserved.")
-  }
+#.onAttach <- function(...){
+#  packageStartupMessage(
+#    "SUNDIALS - Copyright (c) 2002-2020, Lawrence Livermore National Security and Southern Methodist University.\nAll rights reserved.")
+#  }
diff --git a/inst/include/sundials/sundials_config.h b/inst/include/sundials/sundials_config.h
@@ -44,10 +44,10 @@
  * -----------------------------------------------------------------*/
 
 
-#define SUNDIALS_VERSION "7.0.0"
+#define SUNDIALS_VERSION "7.1.1"
 #define SUNDIALS_VERSION_MAJOR 7
-#define SUNDIALS_VERSION_MINOR 0
-#define SUNDIALS_VERSION_PATCH 0
+#define SUNDIALS_VERSION_MINOR 1
+#define SUNDIALS_VERSION_PATCH 1
 #define SUNDIALS_VERSION_LABEL ""
 #define SUNDIALS_GIT_VERSION ""
 
@@ -57,9 +57,10 @@
  * -----------------------------------------------------------------*/
 
 #define SUNDIALS_C_COMPILER_HAS_BUILTIN_EXPECT
-/* #undef SUNDIALS_C_COMPILER_HAS_ATTRIBUTE_ASSUME */
-#define SUNDIALS_C_COMPILER_HAS_BUILTIN_ASSUME
+#define SUNDIALS_C_COMPILER_HAS_ATTRIBUTE_ASSUME
+/* #undef SUNDIALS_C_COMPILER_HAS_BUILTIN_ASSUME */
 /* #undef SUNDIALS_C_COMPILER_HAS_ASSUME */
+#define SUNDIALS_C_COMPILER_HAS_ATTRIBUTE_UNUSED
 
 /* Define precision of SUNDIALS data type 'sunrealtype'
  * Depending on the precision level, one of the following
@@ -104,8 +105,8 @@
 #define SUNDIALS_LOGGING_LEVEL 2
 
 /* Build metadata */
-#define SUN_C_COMPILER "AppleClang"
-#define SUN_C_COMPILER_VERSION "12.0.0.12000032"
+#define SUN_C_COMPILER "GNU"
+#define SUN_C_COMPILER_VERSION "13.2.0"
 #define SUN_C_COMPILER_FLAGS ""
 
 #define SUN_CXX_COMPILER ""
@@ -118,8 +119,8 @@
 
 #define SUN_BUILD_TYPE "RelWithDebInfo"
 
-#define SUN_JOB_ID "20240610073427"
-#define SUN_JOB_START_TIME "20240610073427"
+#define SUN_JOB_ID "20241109090408"
+#define SUN_JOB_START_TIME "20241109090408"
 
 #define SUN_TPL_LIST ""
 #define SUN_TPL_LIST_SIZE ""
@@ -220,7 +221,7 @@
 /* #undef SUNDIALS_GINKGO_BACKENDS_HIP */
 /* #undef SUNDIALS_GINKGO_BACKENDS_OMP */
 /* #undef SUNDIALS_GINKGO_BACKENDS_REF */
-/* #undef SUNDIALS_GINKGO_BACKENDS_DPCPP */
+/* #undef SUNDIALS_GINKGO_BACKENDS_SYCL */
 
 /* MAGMA backends */
 /* #undef SUNDIALS_MAGMA_BACKENDS_CUDA */

diff --git a/src/cvode.cpp b/src/cvode.cpp
@@ -76,7 +76,13 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC,
    // Relative tolerance
    sunrealtype reltol = reltolerance;
 
+   // Set Sundials context
+   SUNContext sunctx;
+   SUNContext_Create(SUN_COMM_NULL, &sunctx);
+
    // Absolute tolerance
+   // Set the vector absolute tolerance -----------------------------------------
+   // abstol must be same length as IC
    int abstol_len = abstolerance.length();
 
    // absolute tolerance is either length == 1 or equal to length of IC
@@ -85,10 +91,6 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC,
      stop("Absolute tolerance must be a scalar or a vector of same length as IC\n");
    }
 
-   // Set the vector absolute tolerance -----------------------------------------
-   // abstol must be same length as IC
-   SUNContext sunctx;
-   SUNContext_Create(SUN_COMM_NULL, &sunctx);
    N_Vector abstol = N_VNew_Serial(y_len, sunctx);
    sunrealtype *abstol_ptr = N_VGetArrayPointer(abstol);
    if(abstol_len == 1){
@@ -97,7 +99,7 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC,
        abstol_ptr[i] = abstolerance[0];
      }
    }
-   else if (abstol_len == y_len){
+   if (abstol_len == y_len){
      for (int i = 0; i<y_len; i++){
        abstol_ptr[i] = abstolerance[i];
      }

diff --git a/src/cvodes.cpp b/src/cvodes.cpp
@@ -115,7 +115,6 @@ int rhs_function_sens(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
 // EwtSet Function. Computes the error weights at the current solution.
 int ewt(N_Vector y, N_Vector w, void *user_data)
 {
-  int i;
   sunrealtype yy, ww, rtol;
   NumericVector atol;
 
@@ -124,14 +123,17 @@ int ewt(N_Vector y, N_Vector w, void *user_data)
   rtol = my_rhs_fun->rtol;
   atol = my_rhs_fun->atol;
 
-  for (i=1; i<=3; i++) {
+  sunindextype y_len = N_VGetLength_Serial(y);
+
+  for (int i=1; i<=y_len; i++) {
+
     yy = NV_Ith_S(y,i-1);
     ww = rtol * SUNRabs(yy) + atol[i-1];
     if (ww <= 0.0) return (-1);
     NV_Ith_S(w,i-1) = 1.0/ww;
-  }
 
-  // Rcout << Ith(w,1) << "\n";
+    // Rcpp::Rcout << NV_Ith_S(w,i-1) << "\n";
+  }
 
   return(0);
 }
@@ -170,34 +172,30 @@ NumericMatrix cvodes(NumericVector time_vector, NumericVector IC,
   // Relative Tolerance
   sunrealtype reltol = reltolerance;
 
-  int abstol_len = abstolerance.length();
   // absolute tolerance is either length == 1 or equal to length of IC
-  // If abstol is not equal to 1 and abstol is not equal to IC, then stop
+  // If abstol is not equal to 1 or abstol is not equal to IC, then stop
+  int abstol_len = abstolerance.length();
   if(abstol_len != 1 && abstol_len != y_len){
      stop("Absolute tolerance must be a scalar or a vector of same length as IC \n");
   }
-
-  // Set the vector absolute tolerance -----------------------------------------
-  // abstol must be same length as IC, if it not 1 or not unequal to y_len
-  SUNContext sunctx;
-  SUNContext_Create(SUN_COMM_NULL, &sunctx);
-  N_Vector abstol = N_VNew_Serial(y_len, sunctx);
-  sunrealtype *abstol_ptr = N_VGetArrayPointer(abstol);
+  // Create an abstolerance vector equal to IC length
+  Rcpp::NumericVector abstol(y_len);
   if(abstol_len == 1){
-
-    // if a scalar is provided - use it to make a vector with same values
     for (int i = 0; i<y_len; i++){
-      abstol_ptr[i] = abstolerance[0];
+      abstol[i] = abstolerance[0];
     }
   }
-  else if (abstol_len == y_len){
+  if (abstol_len == y_len){
     for (int i = 0; i<y_len; i++){
-      abstol_ptr[i] = abstolerance[i];
+      abstol[i] = abstolerance[i];
     }
   }
 
-  //----------------------------------------------------------------------------
-  // // Set the initial conditions-------------------------------------------------
+  // Set context for sundials
+  SUNContext sunctx;
+  SUNContext_Create(SUN_COMM_NULL, &sunctx);
+
+  // Set the initial conditions-------------------------------------------------
   N_Vector y0 = N_VNew_Serial(y_len, sunctx);
   sunrealtype *y0_ptr = N_VGetArrayPointer(y0);
   for (int i = 0; i<y_len; i++){
@@ -247,7 +245,7 @@ NumericMatrix cvodes(NumericVector time_vector, NumericVector IC,
   struct rhs_func_sens my_rhs_function = {input_function,
                                           Parameters,
                                           // as<std::vector<double> >(Parameters),
-                                          reltolerance, abstolerance};
+                                          reltol, abstol};
 
   // setting the user_data in rhs function
   flag = CVodeSetUserData(cvode_mem, (void*)&my_rhs_function);

diff --git a/src/cvsolve.cpp b/src/cvsolve.cpp
@@ -95,7 +95,7 @@ NumericMatrix cvsolve(NumericVector time_vector, NumericVector IC,
       abstol_ptr[i] = abstolerance[0];
     }
   }
-  else if (abstol_len == y_len){
+  if (abstol_len == y_len){
     for (int i = 0; i<abstol_len; i++){
       abstol_ptr[i] = abstolerance[i];
     }