#9 issues with lifetimes

src/jacobian/mod.rs

@@ -34,7 +34,7 @@ pub fn find_non_zeros_nonlinear<F: NonLinearOp + ?Sized>(op: &F, x: &F::V, t: F:
 /// Find the non-zero entries of the matrix of a linear operator.
 /// This is used as the default `find_non_zeros` function for the `NonLinearOp` and `LinearOp` traits.
 /// Users can override this function with a more efficient and reliable implementation if desired.
-pub fn find_non_zeros_linear<F: LinearOp + ?Sized>(op: &F, x: &F::V, t: F::T) -> Vec<(usize, usize)> {
+pub fn find_non_zeros_linear<F: LinearOp + ?Sized>(op: &F, t: F::T) -> Vec<(usize, usize)> {
     let mut v = F::V::zeros(op.nstates());
     let mut col = F::V::zeros(op.nout());
     let mut triplets = Vec::with_capacity(op.nstates());
@@ -59,7 +59,7 @@ pub struct JacobianColoring<M: Matrix> {
 }
 
 impl<M: Matrix> JacobianColoring<M> {
-    pub fn new<F: Op<M=M>>(op: &F, x: &M::V, t: F::T) -> Self {
+    pub fn new<F: Op<M=M>>(op: &F) -> Self {
         let sparsity = op.sparsity().expect("Jacobian sparsity not defined, cannot use coloring");
         let non_zeros = sparsity.indices();
         let ncols = op.nstates();
@@ -99,7 +99,7 @@ impl<M: Matrix> JacobianColoring<M> {
         t: F::T,
         y: &mut F::M,
     ) -> Vec<(usize, usize, F::T)> {
-        let mut triplets = Vec::with_capacity(op.nstates());
+        let triplets = Vec::with_capacity(op.nstates());
         let mut v = F::V::zeros(op.nstates());
         let mut col = F::V::zeros(op.nout());
         for c in 0..self.dst_indices_per_color.len() {
@@ -117,11 +117,10 @@ impl<M: Matrix> JacobianColoring<M> {
     pub fn matrix_inplace<F: LinearOp<M=M, V=M::V, T=M::T>>(
         &self,
         op: &F,
-        x: &F::V,
         t: F::T,
         y: &mut F::M,
     ) -> Vec<(usize, usize, F::T)> {
-        let mut triplets = Vec::with_capacity(op.nstates());
+        let triplets = Vec::with_capacity(op.nstates());
         let mut v = F::V::zeros(op.nstates());
         let mut col = F::V::zeros(op.nout());
         for c in 0..self.dst_indices_per_color.len() {
@@ -183,8 +182,6 @@ mod tests {
         ];
         for triplets in test_triplets {
             let op = helper_triplets2op(triplets.as_slice(), 2, 2);
-            let x = DVector::from_vec(vec![1.0, 1.0]);
-            let t = 0.0;
             let non_zeros = op.sparsity().unwrap().indices();
             let expect = triplets
                 .iter()
@@ -205,8 +202,6 @@ mod tests {
         let expect = vec![vec![1, 1], vec![1, 2], vec![1, 1], vec![1, 2]];
         for (triplets, expect) in test_triplets.iter().zip(expect) {
             let op = helper_triplets2op(triplets.as_slice(), 2, 2);
-            let x = DVector::from_vec(vec![1.0, 1.0]);
-            let t = 0.0;
 
             let non_zeros = op.sparsity().unwrap().indices();
             let ncols = op.nstates();

src/linear_solver/faer/lu.rs

@@ -1,3 +1,5 @@
+use std::rc::Rc;
+
 use crate::{
     linear_solver::LinearSolver, op::linearise::LinearisedOp, solver::SolverProblem, NonLinearOp,
     Op, Scalar, Matrix, LinearOp,
@@ -31,10 +33,9 @@ where
 
 impl<T: Scalar, C: NonLinearOp<M = Mat<T>, V = Col<T>, T = T>> LinearSolver<C> for LU<T, C> {
     fn set_linearisation(&mut self, x: &C::V, t: C::T) {
+        Rc::<LinearisedOp<C>>::get_mut(&mut self.problem.as_mut().expect("Problem not set").f).unwrap().set_x(x);
         let matrix = self.matrix.as_mut().expect("Matrix not set");
-        let problem = self.problem.as_ref().expect("Problem not set");
-        problem.f.set_x(x);
-        problem.f.matrix_inplace(t, matrix);
+        self.problem.as_ref().unwrap().f.matrix_inplace(t, matrix);
         self.lu = Some(matrix.full_piv_lu());
     }
 

src/linear_solver/mod.rs

@@ -1,7 +1,6 @@
 use crate::{op::Op, solver::SolverProblem};
 use anyhow::Result;
 
-pub mod gmres;
 #[cfg(feature = "nalgebra")]
 pub mod nalgebra;
 
@@ -50,27 +49,28 @@ pub mod tests {
     use std::rc::Rc;
 
     use crate::{
-        linear_solver::FaerLU,
-        linear_solver::NalgebraLU,
-        op::{linear_closure::LinearClosure, LinearOp},
+        linear_solver::{FaerLU, NalgebraLU},
+        op::{closure::Closure, NonLinearOp},
         scalar::scale,
         vector::VectorRef,
         DenseMatrix, LinearSolver, SolverProblem, Vector,
     };
-    use num_traits::One;
+    use num_traits::{One, Zero};
 
     use super::LinearSolveSolution;
 
     fn linear_problem<M: DenseMatrix + 'static>() -> (
-        SolverProblem<impl LinearOp<M = M, V = M::V, T = M::T>>,
+        SolverProblem<impl NonLinearOp<M = M, V = M::V, T = M::T>>,
         Vec<LinearSolveSolution<M::V>>,
     ) {
         let diagonal = M::V::from_vec(vec![2.0.into(), 2.0.into()]);
-        let jac = M::from_diagonal(&diagonal);
+        let jac1 = M::from_diagonal(&diagonal);
+        let jac2 = M::from_diagonal(&diagonal);
         let p = Rc::new(M::V::zeros(0));
-        let op = Rc::new(LinearClosure::new(
-            // f = J * x + beta * y
-            move |x, _p, _t, beta, y| jac.gemv(M::T::one(), x, beta, y),
+        let op = Rc::new(Closure::new(
+            // f = J * x
+            move |x, _p, _t, y| jac1.gemv(M::T::one(), x, M::T::zero(), y),
+            move |_x, _p, _t, v, y| jac2.gemv(M::T::one(), v, M::T::zero(), y),
             2,
             2,
             p,
@@ -90,7 +90,7 @@ pub mod tests {
         problem: SolverProblem<C>,
         solns: Vec<LinearSolveSolution<C::V>>,
     ) where
-        C: LinearOp,
+        C: NonLinearOp,
         for<'a> &'a C::V: VectorRef<C::V>,
     {
         solver.set_problem(&problem);

src/linear_solver/nalgebra/lu.rs

@@ -1,3 +1,5 @@
+use std::rc::Rc;
+
 use anyhow::Result;
 use nalgebra::{DMatrix, DVector, Dyn};
 
@@ -47,11 +49,10 @@ impl<T: Scalar, C: NonLinearOp<M = DMatrix<T>, V = DVector<T>, T = T>> LinearSol
     }
 
     fn set_linearisation(&mut self, x: &<C as Op>::V, t: <C as Op>::T) {
-        let problem = self.problem.as_ref().expect("Problem not set");
+        Rc::<LinearisedOp<C>>::get_mut(&mut self.problem.as_mut().expect("Problem not set").f).unwrap().set_x(x);
         let matrix = self.matrix.as_mut().expect("Matrix not set");
-        problem.f.set_x(x);
-        problem.f.matrix_inplace(t, matrix);
-        self.lu = Some(matrix.lu());
+        self.problem.as_ref().unwrap().f.matrix_inplace(t, matrix);
+        self.lu = Some(matrix.clone().lu());
     }
 
     fn set_problem(&mut self, problem: &SolverProblem<C>) {

src/linear_solver/sundials.rs

@@ -1,3 +1,5 @@
+use std::rc::Rc;
+
 use anyhow::Result;
 use sundials_sys::{
     realtype, SUNLinSolFree, SUNLinSolSetup, SUNLinSolSolve, SUNLinSol_Dense, SUNLinearSolver,
@@ -77,11 +79,10 @@ where
     }
 
     fn set_linearisation(&mut self, x: &Op::V, t: Op::T) {
+        Rc::<LinearisedOp<Op>>::get_mut(&mut self.problem.as_mut().expect("Problem not set").f).unwrap().set_x(x);
         let matrix = self.matrix.as_mut().expect("Matrix not set");
-        let problem = self.problem.as_ref().expect("Problem not set");
         let linear_solver = self.linear_solver.expect("Linear solver not set");
-        problem.f.set_x(x);
-        problem.f.matrix_inplace(t, matrix);
+        self.problem.as_ref().unwrap().f.matrix_inplace(t, matrix);
         sundials_check(unsafe { SUNLinSolSetup(linear_solver, matrix.sundials_matrix()) }).unwrap();
         self.is_setup = true;
     }

src/matrix/dense_faer_serial.rs

@@ -134,7 +134,7 @@ impl<T: Scalar> Matrix for Mat<T> {
     type Sparsity = Dense;
 
     fn set_data_with_indices(
-        &self,
+        &mut self,
         dst_indices: &<Self::Sparsity as MatrixSparsity>::Index,
         src_indices: &<Self::V as Vector>::Index,
         data: &Self::V,
@@ -176,10 +176,10 @@ impl<T: Scalar> Matrix for Mat<T> {
     }
 
     fn scale_add_and_assign(&mut self, x: &Self, beta: Self::T, y: &Self) {
-        zipped!(self, y).for_each(|unzipped!(mut x, y)| x.write(x.read() + beta * y.read()));
+        zipped!(self, x, y).for_each(|unzipped!(mut s, x, y)| s.write(x.read() + beta * y.read()));
     }
 
-    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, sparsity: Option<&Self::Sparsity>) -> Self {
+    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, _sparsity: Option<&Self::Sparsity>) -> Self {
         Self::zeros(nrows, ncols)
     }
 }

src/matrix/dense_nalgebra_serial.rs

@@ -96,7 +96,7 @@ impl<T: Scalar> Matrix for DMatrix<T> {
     type Sparsity = Dense;
 
     fn set_data_with_indices(
-        &self,
+        &mut self,
         dst_indices: &<Self::Sparsity as super::MatrixSparsity>::Index,
         src_indices: &<Self::V as crate::vector::Vector>::Index,
         data: &Self::V,
@@ -141,7 +141,7 @@ impl<T: Scalar> Matrix for DMatrix<T> {
         self.mul_assign(beta);
         self.add_assign(x);
     }
-    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, sparsity: Option<&Self::Sparsity>) -> Self {
+    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, _sparsity: Option<&Self::Sparsity>) -> Self {
         Self::zeros(nrows, ncols)
     }
 }

src/matrix/mod.rs

@@ -227,7 +227,7 @@ pub trait Matrix:
     fn set_column(&mut self, j: IndexType, v: &Self::V);
 
     fn set_data_with_indices(
-        &self,
+        &mut self,
         dst_indices: &<Self::Sparsity as MatrixSparsity>::Index,
         src_indices: &<Self::V as Vector>::Index,
         data: &Self::V,

src/matrix/sparse_serial.rs

@@ -152,12 +152,12 @@ impl<T: Scalar> Matrix for CscMatrix<T> {
     }
 
     fn set_data_with_indices(
-        &self,
+        &mut self,
         dst_indices: &<Self::Sparsity as MatrixSparsity>::Index,
         src_indices: &<Self::V as crate::vector::Vector>::Index,
         data: &Self::V,
     ) {
-        let values = self.values();
+        let values = self.values_mut();
         for (&dst_i, &src_i) in dst_indices.iter().zip(src_indices.iter()) {
             values[dst_i] = data[src_i];
         }
@@ -203,10 +203,15 @@ impl<T: Scalar> Matrix for CscMatrix<T> {
         ret
     }
     fn set_column(&mut self, j: IndexType, v: &Self::V) {
+        // check v is the same length as the column
+        assert_eq!(v.len(), self.nrows());
+
         let mut col = self.col_mut(j);
-        let (row_indices, values) = col.rows_and_values_mut();
-        for (&mut i, mut d) in row_indices.iter_mut().zip(v.iter()) {
-            *d = v[i];
+        let (dst_row_indices, dst_values) = col.rows_and_values_mut();
+
+        // copy across the non-zero values
+        for (&dst_i, dst_v) in dst_row_indices.iter().zip(dst_values.iter_mut()) {
+            *dst_v = v[dst_i];
         }
     }
     fn scale_add_and_assign(&mut self, x: &Self, beta: Self::T, y: &Self) {

src/matrix/sundials.rs

@@ -230,7 +230,7 @@ impl Matrix for SundialsMatrix {
     type Sparsity = Dense;
 
     fn set_data_with_indices(
-        &self,
+        &mut self,
         dst_indices: &<Self::Sparsity as MatrixSparsity>::Index,
         src_indices: &<Self::V as Vector>::Index,
         data: &Self::V,
@@ -311,7 +311,7 @@ impl Matrix for SundialsMatrix {
         y.axpy(alpha, &tmp, beta);
     }
 
-    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, sparsity: Option<&Self::Sparsity>) -> Self {
+    fn new_from_sparsity(nrows: IndexType, ncols: IndexType, _sparsity: Option<&Self::Sparsity>) -> Self {
         Self::new_dense(nrows, ncols)
     }
 }

src/ode_solver/bdf/mod.rs

@@ -413,7 +413,7 @@ where
 
         // loop until step is accepted
         let y_new = loop {
-            let y_new = y_predict.clone();
+            let mut y_new = y_predict.clone();
 
             // solve BDF equation using y0 as starting point
             let solver_result = self.nonlinear_solver.solve_in_place(&mut y_new, t_new);

src/ode_solver/diffsl.rs

@@ -55,12 +55,10 @@ impl DiffSl {
         };
 
         if use_coloring {
-            let t0 = 0.;
-            let mut y0 = V::zeros(nstates);
-            let rhs = DiffSlRhs(&ret);
-            let mass = DiffSlMass(&ret);
-            ret.rhs_coloring = Some(JacobianColoring::new(&rhs, &y0, t0));
-            ret.mass_coloring = Some(JacobianColoring::new(&mass, &y0, t0));
+            let rhs_coloring = JacobianColoring::new(&DiffSlRhs(&ret));
+            let mass_coloring = JacobianColoring::new(&DiffSlMass(&ret));
+            ret.rhs_coloring = Some(rhs_coloring);
+            ret.mass_coloring = Some(mass_coloring);
         }
         Ok(ret)
     }
@@ -102,7 +100,7 @@ impl NonLinearOp for DiffSlRhs<'_> {
     fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V) {
         self.0.compiler.rhs(
             t,
-            y.as_slice(),
+            x.as_slice(),
             self.0.data.borrow_mut().as_mut_slice(),
             y.as_mut_slice(),
         );
@@ -147,15 +145,15 @@ impl LinearOp for DiffSlMass<'_> {
 
     fn matrix_inplace(&self, t: Self::T, y: &mut Self::M) {
         if let Some(coloring) = &self.0.mass_coloring {
-            let dummy = V::zeros(0);
-            coloring.matrix_inplace(self, &dummy, t, y);
+            coloring.matrix_inplace(self, t, y);
         } else {
             LinearOp::matrix_inplace(self, t, y);
         }
     }
 }
 
-impl OdeEquations for DiffSl {
+impl OdeEquations for DiffSl 
+{
     type M = M;
     type T = T;
     type V = V;

src/ode_solver/equations.rs

@@ -118,7 +118,7 @@ where
         let mut mass = LinearClosure::<M, _>::new(mass, nstates, nstates, p.clone());
         if calculate_sparsity {
             rhs.calculate_sparsity(&y0, t0);
-            mass.calculate_sparsity(&y0, t0);
+            mass.calculate_sparsity(t0);
         }
         Self {
             rhs,