step by step

Author: ZocoLini

CMakeLists.txt

@@ -10,7 +10,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY .)
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY .)
 set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY .)
 
-add_subdirectory(modules/cvodes_old)
 add_subdirectory(modules/sacess)
 add_subdirectory(modules/cuqdyn-c)
 add_subdirectory(modules/cli)

deps.cmake

@@ -1,5 +1,17 @@
 include(FetchContent)
 
+# Building v7.30 of CVODES
+include(ExternalProject)
+
+FetchContent_Declare(
+        cvodes
+        URL https://github.com/LLNL/sundials/releases/download/v7.3.0/cvodes-7.3.0.tar.gz
+        DOWNLOAD_EXTRACT_TIMESTAMP TRUE
+)
+
+FetchContent_MakeAvailable(cvodes)
+
+
 # Building GSL dependency
 find_package(GSL 2.7)
 

modules/cli/src/handlers/solve.c

@@ -17,7 +17,7 @@
 #endif
 
 int handle_solve(int argc, char *argv[]);
-void print_matrix(DlsMat mat, FILE *output_file, char *name);
+void print_matrix(SUNMatrix mat, FILE *output_file, char *name);
 void print_vector(N_Vector vec, FILE *output_file, char *name);
 
 Handler create_solve_handler()
@@ -116,10 +116,10 @@ int handle_solve(int argc, char *argv[])
     {
 
         N_Vector params_median = cuqdyn_result->predicted_params_median;
-        DlsMat data_median = cuqdyn_result->predicted_data_median;
+        SUNMatrix data_median = cuqdyn_result->predicted_data_median;
 
-        DlsMat q_low = cuqdyn_result->q_low;
-        DlsMat q_up = cuqdyn_result->q_up;
+        SUNMatrix q_low = cuqdyn_result->q_low;
+        SUNMatrix q_up = cuqdyn_result->q_up;
         N_Vector times = cuqdyn_result->times;
 
         char *output_file_path = malloc(strlen(output_dir) + strlen("/cuqdyn-results.txt") + 1);
@@ -145,7 +145,7 @@ int handle_solve(int argc, char *argv[])
     return 0;
 }
 
-void print_matrix(DlsMat mat, FILE *output_file, char *name)
+void print_matrix(SUNMatrix mat, FILE *output_file, char *name)
 {
     fprintf(output_file, "[%s]\n", name);
     fprintf(output_file, "%ld %ld\n", SM_ROWS_D(mat), SM_COLUMNS_D(mat));

modules/cuqdyn-c/CMakeLists.txt

@@ -29,6 +29,7 @@ target_include_directories(cuqdyn-c
         PUBLIC
         include
         generated
+        ${cvodes_BINARY_DIR}/include
         ${matio_SOURCE_DIR}/src
         ${muparser_SOURCE_DIR}/include
         ${RUST_LIB_DIR}/include
@@ -37,7 +38,9 @@ target_include_directories(cuqdyn-c
 target_link_directories(cuqdyn-c PUBLIC ${RUST_TARGET_DIR})
 
 target_link_libraries(cuqdyn-c
+        PUBLIC
         sacess
         matio
         mexpreval
+        SUNDIALS::core
 )

modules/cuqdyn-c/include/config.h

@@ -1,16 +1,16 @@
 #ifndef CONFIG_H
 #define CONFIG_H
 
-#include <sundials_old/sundials_nvector.h>
+#include <sundials/sundials_nvector.h>
 
 typedef struct
 {
-    realtype rtol;
+    sunsunrealtype rtol;
     N_Vector atol;
 
 } Tolerances;
 
-Tolerances create_tolerances(realtype rtol, N_Vector atol);
+Tolerances create_tolerances(sunsunrealtype rtol, N_Vector atol);
 void destroy_tolerances(Tolerances);
 
 typedef struct

modules/cuqdyn-c/include/cuqdyn.h

@@ -1,32 +1,27 @@
 #ifndef CUQDYN_H
 #define CUQDYN_H
 
-#include <sundials_old/sundials_nvector.h>
+#include <sundials/sundials_nvector.h>
 
 #include "ode_solver.h"
 
-#define SM_ROWS_D(mat) mat->M
-#define SM_COLUMNS_D(mat) mat->N
-#define SUNDenseMatrix(m, n) NewDenseMat(m, n)
-#define SM_ELEMENT_D(mat, i, j) DENSE_ELEM(mat, i, j)
-#define SM_DATA_D(mat) mat->data
-#define SUNMatDestroy(mat) DestroyMat(mat)
+#define NewDenseMatrix(m, n) SUNDenseMatrix(m, n, get_sundials_ctx())
+#define New_Serial(n) N_VNew_Serial(n, get_sundials_ctx())
 
-#define MIth(v, i) NV_Ith_S(v, i - 1) /* i-th vector component i=1..n */
-#define MIJth(A, i, j) SM_ELEMENT_D(A, i - 1, j - 1) /* (i,j)-th matrix component i,j=1..n */
+SUNContext get_sundials_ctx();
 
 /// Result of the cuqdyn algorithm
 typedef struct
 {
-    DlsMat predicted_data_median;
+    SUNMatrix predicted_data_median;
     N_Vector predicted_params_median;
-    DlsMat q_low;
-    DlsMat q_up;
+    SUNMatrix q_low;
+    SUNMatrix q_up;
     N_Vector times;
 } CuqdynResult;
 
-CuqdynResult* create_cuqdyn_result(DlsMat predicted_data_median, N_Vector predicted_params_median,
-    DlsMat q_low, DlsMat q_up, N_Vector times);
+CuqdynResult* create_cuqdyn_result(SUNMatrix predicted_data_median, N_Vector predicted_params_median,
+    SUNMatrix q_low, SUNMatrix q_up, N_Vector times);
 void destroy_cuqdyn_result(CuqdynResult* result);
 CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, const char *output_file, int rank, int nproc);
 
@@ -42,22 +37,22 @@ long array_get_index(LongArray array, long index);
 
 typedef struct
 {
-    DlsMat *data;
+    SUNMatrix *data;
     long len;
 } MatrixArray;
 
 MatrixArray create_matrix_array(long depth);
 void destroy_matrix_array(MatrixArray array);
-DlsMat matrix_array_get_index(MatrixArray array, long index);
+SUNMatrix matrix_array_get_index(MatrixArray array, long index);
 /// Sets the i-th matrix of the array
-void matrix_array_set_index(MatrixArray array, long index, DlsMat matrix);
+void matrix_array_set_index(MatrixArray array, long index, SUNMatrix matrix);
 /// Returns a matrix where the i-th row and j-th column is the median
 /// of all the i-th row and j-th column of the matrices inside the array
-DlsMat matrix_array_get_median(MatrixArray array);
+SUNMatrix matrix_array_get_median(MatrixArray array);
 /// Resturns a vector containing the i-th row and j-th column of each matrix inside the array
 N_Vector matrix_array_depth_vector_at(MatrixArray array, long i, long j);
 
 /// Creates a vector where the i-th element is the median of the i-th column of the matrix
-N_Vector get_matrix_cols_median(DlsMat matrix);
+N_Vector get_matrix_cols_median(SUNMatrix matrix);
 
 #endif // CUQDYN_H

modules/cuqdyn-c/include/data_reader.h

@@ -1,13 +1,13 @@
 #ifndef FILE_READER_H
 #define FILE_READER_H
 
-#include <sundials_old/sundials_direct.h>
-#include <sundials_old/sundials_nvector.h>
+#include <sundials/sundials_nvector.h>
+#include <sunmatrix/sunmatrix_dense.h>
 
 /* This function wich type of data file is being read (.txt or .mat) and
  * calls the appropriate function to read it.
  */
-int read_data_file(const char *data_file, N_Vector *t, DlsMat *y);
+int read_data_file(const char *data_file, N_Vector *t, SUNMatrix *y);
 
 /*
  * The file should be two ints, m and n, and a matrix mxn
@@ -30,7 +30,7 @@ int read_data_file(const char *data_file, N_Vector *t, DlsMat *y);
  *      .                .
  *      ym0 ym1 ym2 ... ymn
  */
-int read_txt_data_file(const char *data_file, N_Vector *t, DlsMat *y);
+int read_txt_data_file(const char *data_file, N_Vector *t, SUNMatrix *y);
 /*
  * The file should be a single matrix mxn
  * All the data contained in the matrix is the observed data and has this form:
@@ -50,6 +50,6 @@ int read_txt_data_file(const char *data_file, N_Vector *t, DlsMat *y);
  *      .                .
  *      ym0 ym1 ym2 ... ymn
  */
-int read_mat_data_file(const char *data_file, N_Vector *t, DlsMat *y);
+int read_mat_data_file(const char *data_file, N_Vector *t, SUNMatrix *y);
 
 #endif //FILE_READER_H

modules/cuqdyn-c/include/ess_solver.h

@@ -1,11 +1,12 @@
 #ifndef ESS_SOLVER_H
 #define ESS_SOLVER_H
+#include <sundials/sundials_matrix.h>
 
 /// Executes the ESS solver for a given sacess config file using the given objective function.
 /// texp, yexp and initial_values are arguments needed to execute certain objective functions, like
 /// the one defined in lotka_volterra.c.
 /// rank and nproc are used to define the rank of the process and the number of processes in a MPI environment.
 N_Vector execute_ess_solver(const char *config_file, const char *output,
-  N_Vector texp, DlsMat yexp, N_Vector initial_values, int rank, int nproc);
+  N_Vector texp, SUNMatrix yexp, N_Vector initial_values, int rank, int nproc);
 
 #endif //ESS_SOLVER_H

modules/cuqdyn-c/include/functions.h

@@ -1,13 +1,13 @@
 #ifndef LOTKA_VOLTERRA_H
 #define LOTKA_VOLTERRA_H
-#include <sundials_old/sundials_direct.h>
-#include <sundials_old/sundials_nvector.h>
+#include <sundials/sundials_direct.h>
+#include <sundials/sundials_nvector.h>
 
 #include "config.h"
 
 void mexpreval_init_wrapper(OdeExpr ode_expr);
 /// Function used to solve the ODE using cvodes
-int ode_model_fun(realtype t, N_Vector y, N_Vector ydot, void *data);
+int ode_model_fun(sunsunrealtype t, N_Vector y, N_Vector ydot, void *data);
 /// Objetive function for the lotka_volterra problem used by the sacess library
 void* ode_model_obj_func(double *x, void *data);
 

modules/cuqdyn-c/include/matlab.h

@@ -8,34 +8,34 @@
 N_Vector copy_vector_remove_indices(N_Vector original, LongArray indices);
 // Return a copy of a matrix removing the indicated rows.
 // The indices should be sorted in ascending order and 1-based.
-DlsMat copy_matrix_remove_rows(DlsMat matrix, LongArray indices);
+SUNMatrix copy_matrix_remove_rows(SUNMatrix matrix, LongArray indices);
 // Return a copy of a matrix removing the indicated cols.
 // // The indices should be sorted in ascending order and 1-based.
-DlsMat copy_matrix_remove_columns(DlsMat matrix, LongArray indices);
+SUNMatrix copy_matrix_remove_columns(SUNMatrix matrix, LongArray indices);
 // Return a copy of a matrix removing the indicated rows and cols.
 // // The indices should be sorted in ascending order and 1-based.
-DlsMat copy_matrix_remove_rows_and_columns(DlsMat matrix, LongArray row_indices, LongArray col_indices);
+SUNMatrix copy_matrix_remove_rows_and_columns(SUNMatrix matrix, LongArray row_indices, LongArray col_indices);
 
 // Set a row of a matrix with the values of a vector from start (inclusive) to end (exclusive)
-void set_matrix_row(DlsMat matrix, N_Vector vec, long row_index, long start, long end);
+void set_matrix_row(SUNMatrix matrix, N_Vector vec, long row_index, long start, long end);
 // Set a column of a matrix with the values of a vector from start (inclusive) to end (exclusive)
-void set_matrix_column(DlsMat matrix, N_Vector vec, long col_index, long start, long end);
+void set_matrix_column(SUNMatrix matrix, N_Vector vec, long col_index, long start, long end);
 
 // Copy a row of a matrix into a vector from start (inclusive) to end (exclusive)
-N_Vector copy_matrix_row(DlsMat matrix, long row_index, long start, long end);
+N_Vector copy_matrix_row(SUNMatrix matrix, long row_index, long start, long end);
 // Copy a column of a matrix into a vector from start (inclusive) to end (exclusive)
-N_Vector copy_matrix_column(DlsMat matrix, long col_index, long start, long end);
+N_Vector copy_matrix_column(SUNMatrix matrix, long col_index, long start, long end);
 
 /*
  * Sums matrix a + b. Note that 'b' can have fewer cols or rows than 'a' but not otherwise
  * Returns NULL if the dimensions are not compatible
  */
-DlsMat sum_two_matrices(DlsMat a, DlsMat b);
+SUNMatrix sum_two_matrices(SUNMatrix a, SUNMatrix b);
 /*
  * Substracts b from a. Note that 'b' can have fewer cols or rows than 'a' but not otherwise
  * Returns NULL if the dimensions are not compatible
  */
-DlsMat subtract_two_matrices(DlsMat a, DlsMat b);
+SUNMatrix subtract_two_matrices(SUNMatrix a, SUNMatrix b);
 
 /*
  * For an n-element vector A, quantile computes quantiles by using a sorting-based algorithm:
@@ -53,5 +53,5 @@ DlsMat subtract_two_matrices(DlsMat a, DlsMat b);
  * For the quantiles corresponding to the probabilities outside that range,
  * quantile assigns the minimum or maximum values of the elements in A.
  */
-realtype quantile(N_Vector, realtype);
+sunsunrealtype quantile(N_Vector, sunsunrealtype);
 #endif //MATLAB_H