add diffsl dep

.cargo/config.toml

@@ -0,0 +1,2 @@
+[net]
+git-fetch-with-cli = true

.gitconfig

@@ -0,0 +1,2 @@
+[credential]
+  helper = store

.github/workflows/rust.yml

@@ -11,20 +11,37 @@ env:
 
 jobs:
   build:
-    name: Rust project - latest
+    name: Rust project - ${{ matrix.os }} - ${{ matrix.toolchain }} - ${{ matrix.cargo_args }}
     runs-on: ubuntu-latest
     strategy:
         matrix:
           toolchain:
             - stable
-            - beta
-            # - nightly
+            #- beta
+            #- nightly
+          os:
+            - ubuntu-latest
+            - macos-latest
+            - windows-latest
+          cargo_args:
+            - ''
+            - '--features diffsl-llvm14'
+          exclude:
+            - toolchain: beta
+              cargo_args: '--features diffsl-llvm14'
+            - toolchain: nightly
+              cargo_args: '--features diffsl-llvm14'
 
     steps:
     - uses: actions/checkout@v3
     - name: Set up Rust
       run: rustup update ${{ matrix.toolchain }} && rustup default ${{ matrix.toolchain }}
+    - name: Install LLVM and Clang
+      uses: KyleMayes/install-llvm-action@v1
+      if: matrix.cargo_args == '--features diffsl-llvm14'
+      with:
+        version: "14.0"
     - name: Build
-      run: cargo build --verbose
+      run: cargo build --verbose ${{ matrix.cargo_args }}
     - name: Run tests
-      run: cargo test --verbose
+      run: cargo test --verbose ${{ matrix.cargo_args }}

.gitignore

@@ -18,3 +18,5 @@ Cargo.lock
 
 /target
 /Cargo.lock
+
+.vscode

Cargo.toml

@@ -8,17 +8,45 @@ license = "MIT"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
+[features]
+diffsl = []
+diffsl-llvm4 = ["diffsl4-0", "diffsl"]
+diffsl-llvm5 = ["diffsl5-0", "diffsl"]
+diffsl-llvm6 = ["diffsl6-0", "diffsl"]
+diffsl-llvm7 = ["diffsl7-0", "diffsl"]
+diffsl-llvm8 = ["diffsl8-0", "diffsl"]
+diffsl-llvm9 = ["diffsl9-0", "diffsl"]
+diffsl-llvm10 = ["diffsl10-0", "diffsl"]
+diffsl-llvm11 = ["diffsl11-0", "diffsl"]
+diffsl-llvm12 = ["diffsl12-0", "diffsl"]
+diffsl-llvm13 = ["diffsl13-0", "diffsl"]
+diffsl-llvm14 = ["diffsl14-0", "diffsl"]
+diffsl-llvm15 = ["diffsl15-0", "diffsl"]
+diffsl-llvm16 = ["diffsl16-0", "diffsl"]
+diffsl-llvm17 = ["diffsl17-0", "diffsl"]
+
 [dependencies]
 nalgebra = ">=0.32"
 nalgebra-sparse = ">=0.9"
 anyhow = ">=1.0.77"
 num-traits = "0.2.17"
 ouroboros = "0.18.2"
 serde = { version = "1.0.196", features = ["derive"] }
-iree-rs = { verson = "0.1.1", optional = true }
+diffsl4-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm4-0"], optional = true }
+diffsl5-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm5-0"], optional = true }
+diffsl6-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm6-0"], optional = true }
+diffsl7-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm7-0"], optional = true }
+diffsl8-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm8-0"], optional = true }
+diffsl9-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm9-0"], optional = true }
+diffsl10-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm10-0"], optional = true }
+diffsl11-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm11-0"], optional = true }
+diffsl12-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm12-0"], optional = true }
+diffsl13-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm13-0"], optional = true }
+diffsl14-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm14-0"], optional = true }
+diffsl15-0 = { package = "diffsl", git = "https://github.com/martinjrobins/diffsl.git", version = ">=0.1.1", features = ["llvm15-0"], optional = true }
+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 }
 
-[dev-dependencies]
-insta = { version = "1.34.0", features = ["yaml"] }
 
-[features]
-iree = ["dep:iree-rs"]
+[dev-dependencies]
+insta = { version = "1.34.0", features = ["yaml"] }

src/lib.rs

@@ -1,3 +1,31 @@
+#[cfg(feature = "diffsl-llvm4")]
+pub extern crate diffsl4_0 as diffsl;
+#[cfg(feature = "diffsl-llvm5")]
+pub extern crate diffsl5_0 as diffsl;
+#[cfg(feature = "diffsl-llvm6")]
+pub extern crate diffsl6_0 as diffsl;
+#[cfg(feature = "diffsl-llvm7")]
+pub extern crate diffsl7_0 as diffsl;
+#[cfg(feature = "diffsl-llvm8")]
+pub extern crate diffsl8_0 as diffsl;
+#[cfg(feature = "diffsl-llvm9")]
+pub extern crate diffsl9_0 as diffsl;
+#[cfg(feature = "diffsl-llvm10")]
+pub extern crate diffsl10_0 as diffsl;
+#[cfg(feature = "diffsl-llvm11")]
+pub extern crate diffsl11_0 as diffsl;
+#[cfg(feature = "diffsl-llvm12")]
+pub extern crate diffsl12_0 as diffsl;
+#[cfg(feature = "diffsl-llvm13")]
+pub extern crate diffsl13_0 as diffsl;
+#[cfg(feature = "diffsl-llvm14")]
+pub extern crate diffsl14_0 as diffsl;
+#[cfg(feature = "diffsl-llvm15")]
+pub extern crate diffsl15_0 as diffsl;
+#[cfg(feature = "diffsl-llvm16")]
+pub extern crate diffsl16_0 as diffsl;
+#[cfg(feature = "diffsl-llvm17")]
+pub extern crate diffsl17_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;
@@ -27,14 +55,13 @@ use nalgebra::{ClosedSub, ClosedMul, ClosedDiv, ClosedAdd, SimdRealField, Comple
 use num_traits::{Signed, Pow};
 use vector::{Vector, VectorView, VectorViewMut, VectorIndex, VectorRef};
 use nonlinear_solver::{NonLinearSolver, newton::NewtonNonlinearSolver};
-use op::{NonLinearOp, LinearOp, ConstantOp};
 use matrix::{DenseMatrix, MatrixViewMut, Matrix};
+use op::{LinearOp, NonLinearOp};
 use solver::SolverProblem;
 use linear_solver::{lu::LU, LinearSolver};
-pub use ode_solver::{OdeSolverProblem, OdeSolverState, bdf::Bdf, OdeSolverMethod};
+pub use ode_solver::{OdeSolverProblem, OdeSolverState, bdf::Bdf, OdeSolverMethod, equations::OdeEquations};
 
 
-#[cfg(test)]
 #[cfg(test)]
 mod tests {
     use std::rc::Rc;

src/linear_solver/gmres.rs

@@ -1,6 +1,6 @@
 use anyhow::Result;
 
-use crate::{LinearSolver, vector::VectorRef, LinearOp, SolverProblem};
+use crate::{LinearSolver, vector::VectorRef, op::LinearOp, SolverProblem};
 
 
 pub struct GMRES<C: LinearOp> 

src/linear_solver/lu.rs

@@ -50,7 +50,7 @@ impl<T: Scalar, C: LinearOp<M = DMatrix<T>, V = DVector<T>, T = T>> LinearSolver
     }
 
     fn set_problem(&mut self, problem: SolverProblem<C>) {
-        self.lu = Some(nalgebra::LU::new(problem.f.jacobian(&problem.p, problem.t)));
+        self.lu = Some(nalgebra::LU::new(problem.f.jacobian(problem.t)));
         self.problem = Some(problem);
     }
 }

src/linear_solver/mod.rs

@@ -60,16 +60,16 @@ pub mod tests {
     fn linear_problem<M: DenseMatrix + 'static>() -> (SolverProblem<impl LinearOp<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 p = Rc::new(M::V::zeros(0));
         let op = Rc::new(LinearClosure::new(
             // f = J * x
             move | x, _p, _t, y | jac.gemv(M::T::one(), x, M::T::zero(), y),
-            2, 2, 0
+            2, 2, p
         ));
-        let p = Rc::new(M::V::zeros(0));
         let t = M::T::zero();
         let rtol = M::T::from(1e-6);
         let atol = Rc::new(M::V::from_vec(vec![1e-6.into(), 1e-6.into()]));
-        let problem = SolverProblem::new(op, p, t, atol, rtol);
+        let problem = SolverProblem::new(op, t, atol, rtol);
         let solns = vec![
             LinearSolveSolution::new(M::V::from_vec(vec![2.0.into(), 4.0.into()]), M::V::from_vec(vec![1.0.into(), 2.0.into()]))
         ];

src/matrix/mod.rs

@@ -1,6 +1,7 @@
 use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign};
 use std::fmt::Debug;
 
+use num_traits::{One, Zero};
 use crate::{IndexType, Scalar, Vector};
 use anyhow::Result;
 
@@ -12,15 +13,11 @@ pub trait MatrixCommon: Sized + Debug
     type V: Vector<T = Self::T>;
     type T: Scalar;
 
-
-
     /// Get the number of columns of the matrix
     fn nrows(&self) -> IndexType;
 
     /// Get the number of rows of the matrix
     fn ncols(&self) -> IndexType;
-
-
 }
 
 impl <'a, M> MatrixCommon for &'a M where M: MatrixCommon {
@@ -49,13 +46,11 @@ impl <'a, M> MatrixCommon for &'a mut M where M: MatrixCommon {
 pub trait MatrixOpsByValue<Rhs = Self, Output = Self>: MatrixCommon 
     + Add<Rhs, Output = Output>
     + Sub<Rhs, Output = Output> 
-    + Mul<Rhs, Output = Output> 
 {}
 
 impl <M, Rhs, Output> MatrixOpsByValue<Rhs, Output> for M where M: MatrixCommon 
     + Add<Rhs, Output = Output>
     + Sub<Rhs, Output = Output> 
-    + Mul<Rhs, Output = Output> 
 {}
 
 pub trait MatrixMutOpsByValue<Rhs = Self>: MatrixCommon 
@@ -68,42 +63,10 @@ impl <M, Rhs> MatrixMutOpsByValue<Rhs> for M where M: MatrixCommon
     + SubAssign<Rhs> 
 {}
 
-pub trait MatrixMutOps<Other>: 
-    MatrixMutOpsByValue<Other> 
-    + for<'a> MatrixMutOpsByValue<&'a Other>
-    + MulAssign<Self::T>
-    + DivAssign<Self::T>
-{}
-
-impl <M, View> MatrixMutOps<View> for M 
-where 
-    M: MatrixMutOpsByValue<Self> 
-    + for<'a> MatrixMutOpsByValue<&'a Self> 
-    + MatrixMutOpsByValue<View> 
-    + for<'a> MatrixMutOpsByValue<&'a View>
-    + MulAssign<Self::T>
-    + DivAssign<Self::T>
-{}
-
-
-pub trait MatrixOps<Rhs>: 
-    MatrixOpsByValue<Rhs> 
-    + for<'a> MatrixOpsByValue<&'a Rhs>
-    + Mul<Self::T, Output = Self>
-    + Div<Self::T, Output = Self>
-{}
-
-impl <M, Rhs> MatrixOps<Rhs> for M where 
-    M: MatrixOpsByValue<Rhs> 
-    + for<'a> MatrixOpsByValue<&'a Rhs>
-    + Mul<Self::T, Output = M>
-    + Div<Self::T, Output = M>
-{}
-
 /// A trait allowing for references to implement matrix operations
 pub trait MatrixRef<M: MatrixCommon>:
     MatrixOpsByValue<M, M>
-    + for<'a> MatrixOpsByValue<M, M> 
+    + for<'a> MatrixOpsByValue<&'a M, M> 
     + Mul<M::T, Output = M>
     + Div<M::T, Output = M>
 {}
@@ -118,28 +81,32 @@ impl <RefT, M: MatrixCommon> MatrixRef<M> for RefT where
 
 /// A mutable view of a dense matrix [Matrix]
 pub trait MatrixViewMut<'a>: 
-    MatrixMutOps<Self>
-    + MatrixMutOps<Self::View>
-    where Self: 'a
+    for<'b> MatrixMutOpsByValue<&'b Self>
+    + for<'b> MatrixMutOpsByValue<&'b Self::View>
+    + MulAssign<Self::T>
+    + DivAssign<Self::T>
 {
-    type Owned: DenseMatrix<T = Self::T, V = Self::V, ViewMut<'a> = Self>;
-    type View: MatrixView<'a, Owned = Self::Owned, T = Self::T>;
+    type Owned;
+    type View;
     fn gemm_oo(&mut self, alpha: Self::T, a: &Self::Owned, b: &Self::Owned, beta: Self::T);
     fn gemm_vo(&mut self, alpha: Self::T, a: &Self::View, b: &Self::Owned, beta: Self::T);
 }
 
 /// A view of a dense matrix [Matrix]
 pub trait MatrixView<'a>: 
-    MatrixRef<Self::Owned>
-    + Clone 
-    where Self: 'a
+    for<'b> MatrixOpsByValue<&'b Self::Owned, Self::Owned> 
+    + Mul<Self::T, Output = Self::Owned>
+    + Div<Self::T, Output = Self::Owned>
+    + Clone
 {
-    type Owned: DenseMatrix<T = Self::T, V = Self::V>;
+    type Owned;
 }
 
 /// A base matrix trait (including sparse and dense matrices)
 pub trait Matrix: 
-    MatrixOps<Self>
+    for<'a> MatrixOpsByValue<&'a Self, Self>
+    + Mul<Self::T, Output = Self>
+    + Div<Self::T, Output = Self>
     + Clone 
 {
     /// Extract the diagonal of the matrix as an owned vector
@@ -159,17 +126,17 @@ pub trait Matrix:
 /// A dense column-major matrix. The assumption is that the underlying matrix is stored in column-major order, so functions for taking columns views are efficient
 pub trait DenseMatrix: 
     Matrix
-    + for <'a> MatrixOps<Self::View<'a>>
-    + for <'a> MatrixMutOps<Self::View<'a>>
+    + for<'a, 'b> MatrixOpsByValue<&'b Self::View<'a>, Self>
+    + for<'a, 'b> MatrixMutOpsByValue<&'b Self::View<'a>>
     + Index<(IndexType, IndexType), Output = Self::T> 
     + IndexMut<(IndexType, IndexType), Output = Self::T> 
 {
 
     /// A view of the dense matrix type
-    type View<'a>: MatrixView<'a, Owned = Self, T = Self::T> where Self: 'a;
+    type View<'a>: MatrixView<'a, Owned = Self, T = Self::T, V = Self::V> where Self: 'a;
 
     /// A mutable view of the dense matrix type
-    type ViewMut<'a>: MatrixViewMut<'a, Owned = Self, T = Self::T, View = Self::View<'a>> where Self: 'a;
+    type ViewMut<'a>: MatrixViewMut<'a, Owned = Self, T = Self::T, V = Self::V, View = Self::View<'a>> where Self: 'a;
 
 
     /// Perform a matrix-matrix multiplication `self = alpha * a * b + beta * self`, where `alpha` and `beta` are scalars, and `a` and `b` are matrices
@@ -192,6 +159,14 @@ pub trait DenseMatrix:
     /// Get a mutable vector view of the column `i`
     fn column_mut(&mut self, i: IndexType) -> <Self::V as Vector>::ViewMut<'_>;
 
+    /// mat_mat_mul using gemm, allocating a new matrix
+    fn mat_mul(&self, b: &Self) -> Self {
+        let nrows = self.nrows();
+        let ncols = b.ncols();
+        let mut ret = Self::zeros(nrows, ncols);
+        ret.gemm(Self::T::one(), self, b, Self::T::zero());
+        ret
+    }
 
 }