Merge pull request #12 from rustodeans/faer

Adding support for faer

Cargo.toml

@@ -9,6 +9,9 @@ license = "MIT"
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [features]
+default = ["nalgebra", "faer"]
+faer = []
+nalgebra = []
 diffsl = []
 diffsl-llvm4 = ["diffsl4-0", "diffsl"]
 diffsl-llvm5 = ["diffsl5-0", "diffsl"]
@@ -47,6 +50,7 @@ diffsl15-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffs
 diffsl16-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm16-0"], optional = true }
 diffsl17-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm17-0"], optional = true }
 petgraph = "0.6.4"
+faer = "0.18.2"
 
 
 [dev-dependencies]

src/lib.rs

@@ -27,60 +27,24 @@ pub extern crate diffsl8_0 as diffsl;
 #[cfg(feature = "diffsl-llvm9")]
 pub extern crate diffsl9_0 as diffsl;
 
-pub trait Scalar:
-    nalgebra::Scalar
-    + From<f64>
-    + Display
-    + SimdRealField
-    + ComplexField
-    + Copy
-    + ClosedSub
-    + From<f64>
-    + ClosedMul
-    + ClosedDiv
-    + ClosedAdd
-    + Signed
-    + PartialOrd
-    + Pow<Self, Output = Self>
-    + Pow<i32, Output = Self>
-{
-    const EPSILON: Self;
-    const INFINITY: Self;
-    const NAN: Self;
-    fn is_nan(self) -> bool;
-}
-
-type IndexType = usize;
-
-impl Scalar for f64 {
-    const EPSILON: Self = f64::EPSILON;
-    const INFINITY: Self = f64::INFINITY;
-    const NAN: Self = f64::NAN;
-    fn is_nan(self) -> bool {
-        self.is_nan()
-    }
-}
-
 pub mod jacobian;
 pub mod linear_solver;
 pub mod matrix;
 pub mod nonlinear_solver;
 pub mod ode_solver;
 pub mod op;
+pub mod scalar;
 pub mod solver;
 pub mod vector;
 
-use std::fmt::Display;
-
 use linear_solver::{lu::LU, LinearSolver};
 use matrix::{DenseMatrix, Matrix, MatrixViewMut};
-use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField};
 use nonlinear_solver::{newton::NewtonNonlinearSolver, NonLinearSolver};
-use num_traits::{Pow, Signed};
 pub use ode_solver::{
     bdf::Bdf, equations::OdeEquations, OdeSolverMethod, OdeSolverProblem, OdeSolverState,
 };
 use op::NonLinearOp;
+use scalar::{IndexType, Scalar};
 use solver::SolverProblem;
 use vector::{Vector, VectorIndex, VectorRef, VectorView, VectorViewMut};
 

src/matrix/dense_serial.rs

@@ -31,6 +31,7 @@ impl<'a, T: Scalar> MatrixViewMut<'a> for DMatrixViewMut<'a, T> {
 impl<'a, T: Scalar> MatrixCommon for DMatrixView<'a, T> {
     type V = DVector<T>;
     type T = T;
+
     fn ncols(&self) -> IndexType {
         self.ncols()
     }

src/matrix/mod.rs

@@ -40,6 +40,7 @@ where
 {
     type T = M::T;
     type V = M::V;
+
     fn ncols(&self) -> IndexType {
         M::ncols(self)
     }

src/nonlinear_solver/mod.rs

@@ -3,7 +3,7 @@ use core::panic;
 use num_traits::{One, Pow};
 use std::rc::Rc;
 
-use crate::{op::Op, solver::SolverProblem, IndexType, Scalar, Vector};
+use crate::{op::Op, scalar::scale, solver::SolverProblem, IndexType, Scalar, Vector};
 
 pub struct NonLinearSolveSolution<V> {
     pub x0: V,
@@ -84,7 +84,7 @@ impl<C: Op> Convergence<C> {
         }
     }
     fn reset(&mut self, y: &C::V) {
-        let mut scale = y.abs() * self.rtol;
+        let mut scale = y.abs() * scale(self.rtol);
         scale += self.atol.as_ref();
         self.scale = Some(scale);
         self.iter = 0;
@@ -193,7 +193,7 @@ pub mod tests {
             let x = solver.solve(&soln.x0).unwrap();
             let tol = {
                 let problem = solver.problem().unwrap();
-                soln.x.abs() * problem.rtol + problem.atol.as_ref()
+                soln.x.abs() * scale(problem.rtol) + problem.atol.as_ref()
             };
             x.assert_eq(&soln.x, tol[0]);
         }

src/ode_solver/bdf.rs

@@ -7,7 +7,7 @@ use num_traits::{One, Pow, Zero};
 use serde::Serialize;
 
 use crate::{
-    matrix::MatrixRef, op::ode::BdfCallable, DenseMatrix, IndexType, MatrixViewMut,
+    matrix::MatrixRef, op::ode::BdfCallable, scalar::scale, DenseMatrix, IndexType, MatrixViewMut,
     NewtonNonlinearSolver, NonLinearSolver, Scalar, SolverProblem, Vector, VectorRef, VectorView,
     VectorViewMut, LU,
 };
@@ -245,15 +245,15 @@ where
         for i in 0..self.order {
             let i_t = Eqn::T::from(i as f64);
             time_factor *= (t - (state.t - state.h * i_t)) / (state.h * (Eqn::T::one() + i_t));
-            order_summation += self.diff.column(i + 1) * time_factor;
+            order_summation += self.diff.column(i + 1) * scale(time_factor);
         }
         order_summation
     }
 
     fn problem(&self) -> Option<&OdeSolverProblem<Eqn>> {
         self.ode_problem.as_ref()
     }
-  
+
     fn set_problem(&mut self, state: &mut OdeSolverState<Eqn::M>, problem: &OdeSolverProblem<Eqn>) {
         self.ode_problem = Some(problem.clone());
         let nstates = problem.eqn.nstates();
@@ -289,9 +289,9 @@ where
         }
 
         // update initial step size based on function
-        let mut scale = state.y.abs();
-        scale *= problem.rtol;
-        scale += problem.atol.as_ref();
+        let mut scale_factor = state.y.abs();
+        scale_factor *= scale(problem.rtol);
+        scale_factor += problem.atol.as_ref();
 
         let f0 = problem.eqn.rhs(state.t, &state.y);
         let hf0 = &f0 * state.h;
@@ -303,7 +303,7 @@ where
         self.diff.column_mut(1).copy_from(&hf0);
 
         let mut df = f1 - f0;
-        df.component_div_assign(&scale);
+        df.component_div_assign(&scale_factor);
         let d2 = df.norm();
 
         let one_over_order_plus_one =
@@ -355,7 +355,7 @@ where
                     // test error is within tolerance
                     {
                         let ode_problem = self.ode_problem.as_ref().unwrap();
-                        scale_y = y_new.abs() * ode_problem.rtol;
+                        scale_y = y_new.abs() * scale(ode_problem.rtol);
                         scale_y += ode_problem.atol.as_ref();
                     }
 
@@ -446,14 +446,14 @@ where
             // order k, we need to calculate the optimal step size factors for orders
             // k-1 and k+1. To do this, we note that the error = C_k * D^{k+1} y_n
             let error_m_norm = if order > 1 {
-                let mut error_m = self.diff.column(order) * self.error_const[order - 1];
+                let mut error_m = self.diff.column(order) * scale(self.error_const[order - 1]);
                 error_m.component_div_assign(&scale_y);
                 error_m.norm()
             } else {
                 Eqn::T::INFINITY
             };
             let error_p_norm = if order < Self::MAX_ORDER {
-                let mut error_p = self.diff.column(order + 2) * self.error_const[order + 1];
+                let mut error_p = self.diff.column(order + 2) * scale(self.error_const[order + 1]);
                 error_p.component_div_assign(&scale_y);
                 error_p.norm()
             } else {

src/ode_solver/mod.rs

@@ -363,6 +363,7 @@ mod tests {
     };
     use super::*;
     use crate::nonlinear_solver::newton::NewtonNonlinearSolver;
+    use crate::scalar::scale;
     use tests::bdf::Bdf;
     use tests::test_models::dydt_y2::dydt_y2_problem;
     use tests::test_models::gaussian_decay::gaussian_decay_problem;
@@ -391,7 +392,7 @@ mod tests {
             } else {
                 let tol = {
                     let problem = method.problem().unwrap();
-                    soln.abs() * problem.rtol + problem.atol.as_ref()
+                    soln.abs() * scale(problem.rtol) + problem.atol.as_ref()
                 };
                 soln.assert_eq(&point.state, tol[0]);
             }

src/ode_solver/test_models/dydt_y2.rs

@@ -1,5 +1,6 @@
 use crate::{
     ode_solver::{OdeBuilder, OdeSolverProblem, OdeSolverSolution},
+    scalar::scale,
     DenseMatrix, OdeEquations, Vector,
 };
 use num_traits::One;
@@ -15,7 +16,7 @@ fn rhs<M: DenseMatrix>(x: &M::V, _p: &M::V, _t: M::T, y: &mut M::V) {
 fn rhs_jac<M: DenseMatrix>(x: &M::V, _p: &M::V, _t: M::T, v: &M::V, y: &mut M::V) {
     y.copy_from(v);
     y.component_mul_assign(x);
-    y.mul_assign(M::T::from(2.));
+    y.mul_assign(scale(M::T::from(2.)));
 }
 
 pub fn dydt_y2_problem<M: DenseMatrix + 'static>(
@@ -42,7 +43,7 @@ pub fn dydt_y2_problem<M: DenseMatrix + 'static>(
     for i in 0..=n {
         let t = M::T::from(i as f64 * dt);
         // y = y0 / (1 - y0 * t)
-        let mut denom = y0.clone() * (-t);
+        let mut denom = y0.clone() * (scale(-t));
         denom.add_scalar_mut(M::T::one());
         let mut y = y0.clone();
         y.component_div_assign(&denom);

src/ode_solver/test_models/exponential_decay.rs

@@ -1,5 +1,6 @@
 use crate::{
     ode_solver::{OdeBuilder, OdeSolverProblem, OdeSolverSolution},
+    scalar::scale,
     DenseMatrix, OdeEquations, Vector,
 };
 use nalgebra::ComplexField;
@@ -10,7 +11,7 @@ use std::ops::MulAssign;
 // dy/dt = -ay (p = [a])
 fn exponential_decay<M: DenseMatrix>(x: &M::V, p: &M::V, _t: M::T, y: &mut M::V) {
     y.copy_from(x);
-    y.mul_assign(-p[0]);
+    y.mul_assign(scale(-p[0]));
 }
 
 // Jv = -av
@@ -22,7 +23,7 @@ fn exponential_decay_jacobian<M: DenseMatrix>(
     y: &mut M::V,
 ) {
     y.copy_from(v);
-    y.mul_assign(-p[0]);
+    y.mul_assign(scale(-p[0]));
 }
 
 fn exponential_decay_init<M: DenseMatrix>(_p: &M::V, _t: M::T) -> M::V {
@@ -49,7 +50,7 @@ pub fn exponential_decay_problem<M: DenseMatrix + 'static>(
     for i in 0..10 {
         let t = M::T::from(i as f64 / 10.0);
         let y0: M::V = problem.eqn.init(M::T::zero());
-        let y = y0 * M::T::exp(-p[0] * t);
+        let y = y0 * scale(M::T::exp(-p[0] * t));
         soln.push(y, t);
     }
     (problem, soln)

src/ode_solver/test_models/exponential_decay_with_algebraic.rs

@@ -1,5 +1,6 @@
 use crate::{
     ode_solver::{OdeBuilder, OdeSolverProblem, OdeSolverSolution},
+    scalar::scale,
     DenseMatrix, OdeEquations, Vector,
 };
 use nalgebra::ComplexField;
@@ -18,7 +19,7 @@ fn exponential_decay_with_algebraic<M: DenseMatrix>(
     mut y: &mut M::V,
 ) {
     y.copy_from(x);
-    y.mul_assign(-p[0]);
+    y.mul_assign(scale(-p[0]));
     let nstates = y.len();
     y[nstates - 1] = x[nstates - 1] - x[nstates - 2];
 }
@@ -33,7 +34,7 @@ fn exponential_decay_with_algebraic_jacobian<M: DenseMatrix>(
     mut y: &mut M::V,
 ) {
     y.copy_from(v);
-    y.mul_assign(-p[0]);
+    y.mul_assign(scale(-p[0]));
     let nstates = y.len();
     y[nstates - 1] = v[nstates - 1] - v[nstates - 2];
 }
@@ -75,7 +76,7 @@ pub fn exponential_decay_with_algebraic_problem<M: DenseMatrix + 'static>(
     for i in 0..10 {
         let t = M::T::from(i as f64 / 10.0);
         let y0 = M::V::from_vec(vec![1.0.into(), 1.0.into(), 1.0.into()]);
-        let y: M::V = y0 * M::T::exp(-p[0] * t);
+        let y: M::V = y0 * scale(M::T::exp(-p[0] * t));
         soln.push(y, t);
     }
     (problem, soln)

src/ode_solver/test_models/gaussian_decay.rs

@@ -1,5 +1,6 @@
 use crate::{
     ode_solver::{OdeBuilder, OdeSolverProblem, OdeSolverSolution},
+    scalar::scale,
     DenseMatrix, OdeEquations, Vector,
 };
 use num_traits::Pow;
@@ -10,14 +11,14 @@ use std::ops::MulAssign;
 fn gaussian_decay<M: DenseMatrix>(x: &M::V, p: &M::V, t: M::T, y: &mut M::V) {
     y.copy_from(x);
     y.component_mul_assign(p);
-    y.mul_assign(-t);
+    y.mul_assign(scale(-t));
 }
 
 // Jv = -atv
 fn gaussian_decay_jacobian<M: DenseMatrix>(_x: &M::V, p: &M::V, t: M::T, v: &M::V, y: &mut M::V) {
     y.copy_from(v);
     y.component_mul_assign(p);
-    y.mul_assign(-t);
+    y.mul_assign(scale(-t));
 }
 
 pub fn gaussian_decay_problem<M: DenseMatrix + 'static>(
@@ -41,7 +42,7 @@ pub fn gaussian_decay_problem<M: DenseMatrix + 'static>(
     let mut soln = OdeSolverSolution::default();
     for i in 0..10 {
         let t = M::T::from(i as f64 / 1.0);
-        let px = M::V::from_vec(p.clone()) * t.pow(2) / M::T::from(-2.0);
+        let px = M::V::from_vec(p.clone()) * scale(t.pow(2)) / scale(M::T::from(-2.0));
         let mut y: M::V = problem.eqn.init(M::T::zero());
         y.component_mul_assign(&px.exp());
         soln.push(y, t);

src/op/mod.rs

@@ -1,7 +1,9 @@
 use crate::{
     jacobian::{find_non_zero_entries, find_non_zeros},
+    scalar::scale,
     Matrix, Scalar, Vector,
 };
+
 use num_traits::{One, Zero};
 use std::ops::{AddAssign, MulAssign};
 
@@ -64,9 +66,9 @@ pub trait LinearOp: Op {
     }
     fn gemv(&self, x: &Self::V, t: Self::T, alpha: Self::T, beta: Self::T, y: &mut Self::V) {
         let mut beta_y = y.clone();
-        beta_y.mul_assign(beta);
+        beta_y.mul_assign(scale(beta));
         self.call_inplace(x, t, y);
-        y.mul_assign(alpha);
+        y.mul_assign(scale(alpha));
         y.add_assign(&beta_y);
     }
     fn jacobian_diagonal(&self, t: Self::T) -> Self::V {

src/op/ode.rs

@@ -1,6 +1,7 @@
 use crate::{
     matrix::{DenseMatrix, MatrixRef},
     ode_solver::{equations::OdeEquations, OdeSolverProblem},
+    scalar::scale,
     IndexType, Matrix, Vector, VectorRef,
 };
 use num_traits::{One, Zero};
@@ -92,11 +93,11 @@ impl<Eqn: OdeEquations> BdfCallable<Eqn> {
         y0: &Eqn::V,
     ) {
         // update psi term as defined in second equation on page 9 of [1]
-        let mut new_psi_neg_y0 = diff.column(1) * gamma[1];
+        let mut new_psi_neg_y0 = diff.column(1) * scale(gamma[1]);
         for (i, &gamma_i) in gamma.iter().enumerate().take(order + 1).skip(2) {
-            new_psi_neg_y0 += diff.column(i) * gamma_i
+            new_psi_neg_y0 += diff.column(i) * scale(gamma_i)
         }
-        new_psi_neg_y0 *= alpha[order];
+        new_psi_neg_y0 *= scale(alpha[order]);
 
         // now negate y0
         new_psi_neg_y0.sub_assign(y0);