vector::custom: prefer a clone version that uses the crate "libc"

This is slightly slower than using the Rust allocator but should be
safer.  The speed in comparable to the C version though.

.gitignore

@@ -1,2 +1,4 @@
 /target
 Cargo.lock
+*.svg
+perf.*

Cargo.toml

@@ -25,6 +25,7 @@ nvecpthreads = ["sundials-sys/nvecpthreads"]
 [dependencies]
 sundials-sys = "0.4.0"
 mpi = { version = "0.7.0", optional = true, default-features = false }
+libc = "0.2"
 
 [dev-dependencies]
 eyre = "0.6.8"

examples/speed.c

@@ -0,0 +1,82 @@
+#include <assert.h>
+#include <nvector/nvector_serial.h>
+#include <sundials/sundials_dense.h>
+#include <sunlinsol/sunlinsol_dense.h>
+#include <sunlinsol/sunlinsol_spgmr.h>
+#include <cvode/cvode.h>
+
+int cvrhs(double t, N_Vector u, N_Vector du, void* user_data) {
+  /* double *u_data = N_VGetArrayPointer(u); */
+  double *du_data = N_VGetArrayPointer(du);
+  du_data[0] = 1.;
+  return 0;
+}
+
+int main() {
+  for (int i = 0; i < 100000; i++) {
+    SUNContext ctx;
+    void *cvode_mem;
+    double t0, t, *y0_data, *u_data;
+    N_Vector y0, u;
+    SUNLinearSolver linsol;
+    int ret;
+
+    if (SUNContext_Create(NULL, &ctx) < 0) {
+      printf("SUNContext_Create");
+      return 1;
+    };
+    cvode_mem = CVodeCreate(CV_ADAMS, ctx);    
+    t0 = 0.;
+    y0 = N_VNew_Serial(1, ctx);
+    if (y0 == NULL) {
+      printf("N_VNew_Serial");
+      return 1;
+    }
+    y0_data = N_VGetArrayPointer(y0);
+    y0_data[0] = 0.;
+    ret = CVodeInit(cvode_mem, cvrhs, t0, y0);
+    if (ret != CV_SUCCESS) {
+      printf("CVodeInit");
+      return 1;
+    }
+    CVodeSStolerances(cvode_mem, 1e-6, 1e-12);
+
+    /* Because the dimension is small, a direct solver is more efficient. */
+    /* SUNMatrix mat = SUNDenseMatrix(1, 1, ctx); */
+    /* linsol = SUNLinSol_Dense(y0, mat, ctx); */
+    /* ret = CVodeSetLinearSolver(cvode_mem, linsol, mat); */
+
+    /* Iterative solver (default in Rust) */
+    linsol = SUNLinSol_SPGMR(y0, SUN_PREC_NONE, 30, ctx);
+    ret = CVodeSetLinearSolver(cvode_mem, linsol, NULL);
+
+    if (ret != CV_SUCCESS) {
+      printf("CVodeSetLinearSolver");
+      return 1;
+    }
+    CVodeSetMaxHnilWarns(cvode_mem, 10);
+
+    u = N_VNew_Serial(1, ctx);
+    if (u == NULL) {
+      printf("N_VNew_Serial: u");
+      return 1;
+    }
+    int status = CVode(cvode_mem, 1., u, &t, CV_NORMAL);
+    assert(status == CV_SUCCESS);
+    u_data = N_VGetArrayPointer(u);
+    /* printf("%f\n", u_data[0]); */
+
+    CVodeFree(&cvode_mem);
+    N_VDestroy(y0);
+    N_VDestroy(u);
+    /* SUNMatDestroy(mat); */
+    SUNLinSolFree(linsol);
+    SUNContext_Free(&ctx);
+  }
+  return 0;
+}
+
+
+/* Local Variables: */
+/* compile-command: "gcc -O3 speed.c -o /tmp/speed -lsundials_generic -lsundials_nvecserial -lsundials_cvode -lsundials_sunmatrixdense -lsundials_sunlinsoldense -lsundials_sunlinsolspgmr" */
+/* End: */

examples/speed.rs

@@ -0,0 +1,15 @@
+use std::time::Instant;
+use sundials::{context, cvode::CVode};
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let now = Instant::now();
+    for _ in 0 .. 100_000 {
+        let mut ode = CVode::adams(0., &[0.],
+            |_t, _u, du| *du = [1.])
+            .build(context!()?)?;
+        let mut u = [f64::NAN];
+        ode.solve(1., &mut u);
+    }
+    println!("Elapsed: {}", now.elapsed().as_secs_f64());
+    Ok(())
+}

src/vector/custom.rs

@@ -175,22 +175,36 @@ where V: NVectorOps + 'a {
         let w = Self::ref_of_nvector(nw);
         // Rust memory cannot be uninitialized, thus clone.
         let v = w.clone();
+        //// Sundials functions — slow.
         // let nv = N_VNewEmpty((*nw).sunctx);
         // if N_VCopyOps(nw, nv) != 0 {
-        //     return ptr::null_mut()
+        //     return std::ptr::null_mut()
         // }
         // (*nv).content = Box::into_raw(Box::new(v)) as *mut c_void;
         // nv
-        // Do not go through the C malloc for this (slow).
-        // FIXME: Problems with free?
+
+        //// Do not go through the C malloc for this
+        // FIXME: This is the fastest but are there problems with free?
+        // let sunctx = (*nw).sunctx;
+        // let ops = (*(*nw).ops).clone();
+        // let v = _generic_N_Vector {
+        //     sunctx,
+        //     ops: Box::into_raw(Box::new(ops)),
+        //     content: Box::into_raw(Box::new(v)) as *mut c_void,
+        // };
+        // Box::into_raw(Box::new(v))
+
+        //// libc version — safe as Sundials uses malloc.
         let sunctx = (*nw).sunctx;
-        let ops = (*(*nw).ops).clone();
-        let v = _generic_N_Vector {
-            sunctx,
-            ops: Box::into_raw(Box::new(ops)),
-            content: Box::into_raw(Box::new(v)) as *mut c_void,
-        };
-        Box::into_raw(Box::new(v))
+        let nv = libc::malloc(std::mem::size_of::<_generic_N_Vector>())
+            as N_Vector;
+        (*nv).sunctx = sunctx;
+        let n = std::mem::size_of::<_generic_N_Vector_Ops>();
+        let ops = libc::malloc(n);
+        libc::memcpy(ops, (*nw).ops as *mut c_void, n);
+        (*nv).ops = ops as N_Vector_Ops;
+        (*nv).content = Box::into_raw(Box::new(v)) as *mut c_void;
+        nv
     }
 
     /// Destroys the N_Vector `nv` and frees memory allocated for its