From bdc96237e8cf0661dfacfaf801d7afb6203b3a87 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Sun, 24 Mar 2024 13:32:09 +0000 Subject: [PATCH 01/28] WIP: Vector implementation using faer --- Cargo.toml | 4 ++ src/vector/mod.rs | 17 +++--- src/vector/serial.rs | 6 +-- src/vector/vector_faer.rs | 109 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 126 insertions(+), 10 deletions(-) create mode 100644 src/vector/vector_faer.rs diff --git a/Cargo.toml b/Cargo.toml index e46b9477..0c7a058c 100644 --- a/Cargo.toml +++ b/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 = ["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] diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 0a5bea2f..ff6e774b 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -5,8 +5,10 @@ use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, use crate::{IndexType, Scalar}; mod serial; +mod vector_faer; -pub trait VectorIndex: Sized + Index + Debug + Display { +pub trait VectorIndex: Sized + Index + Debug { + //+ Display fn zeros(len: IndexType) -> Self; fn len(&self) -> IndexType; fn is_empty(&self) -> bool { @@ -14,7 +16,8 @@ pub trait VectorIndex: Sized + Index + Debug + Di } } -pub trait VectorCommon: Sized + Debug + Display { +pub trait VectorCommon: Sized + Debug { + // + Display type T: Scalar; } @@ -72,7 +75,7 @@ pub trait VectorViewMut<'a>: + for<'b> VectorMutOpsByValue<&'b Self::View> + for<'b> VectorMutOpsByValue<&'b Self::Owned> + MulAssign - + DivAssign + // + DivAssign + Index + IndexMut { @@ -89,7 +92,7 @@ pub trait VectorView<'a>: + for<'b> VectorOpsByValue<&'b Self::Owned, Self::Owned> + for<'b> VectorOpsByValue<&'b Self, Self::Owned> + Mul - + Div + // + Div + Index { type Owned; @@ -140,7 +143,7 @@ pub trait Vector: fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T); fn add_scalar_mut(&mut self, scalar: Self::T); fn component_mul_assign(&mut self, other: &Self); - fn component_div_assign(&mut self, other: &Self); + // fn component_div_assign(&mut self, other: &Self); fn filter_indices bool>(&self, f: F) -> Self::Index; fn filter(&self, indices: &Self::Index) -> Self { let mut result = Self::zeros(indices.len()); @@ -164,8 +167,8 @@ pub trait Vector: i, self[i], other[i] ); if self.len() <= 3 { - eprintln!("left: {}", self); - eprintln!("right: {}", other); + eprintln!("left: {:?}", self); + eprintln!("right: {:?}", other); } else if i == 0 { eprintln!( "left: [{}, {}, {}, ...] != [{}, {}, {}, ...]", diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 3c7c6147..9ba7f893 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -92,9 +92,9 @@ impl Vector for DVector { fn axpy(&mut self, alpha: T, x: &Self, beta: T) { self.axpy(alpha, x, beta); } - fn component_div_assign(&mut self, other: &Self) { - self.component_div_assign(other); - } + // fn component_div_assign(&mut self, other: &Self) { + // self.component_div_assign(other); + // } fn component_mul_assign(&mut self, other: &Self) { self.component_mul_assign(other); } diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs new file mode 100644 index 00000000..9903cae9 --- /dev/null +++ b/src/vector/vector_faer.rs @@ -0,0 +1,109 @@ +use faer::{unzipped, zipped}; + +use crate::IndexType; + +use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; +impl Vector for faer::Col { + type View<'a> = faer::ColRef<'a, f64>; + type ViewMut<'a> = faer::ColMut<'a, f64>; + type Index = Vec; + fn len(&self) -> IndexType { + self.nrows() + } + fn norm(&self) -> f64 { + self.norm_l2() + } + fn abs(&self) -> Self { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn as_view(&self) -> Self::View<'_> { + self.as_ref() + } + fn as_view_mut(&mut self) -> Self::ViewMut<'_> { + self.as_mut() + } + fn copy_from(&mut self, other: &Self) { + self.copy_from(other) + } + fn copy_from_view(&mut self, other: &Self::View<'_>) { + self.copy_from(other) + } + fn from_element(nstates: usize, value: Self::T) -> Self { + faer::Col::from_vec(vec![value; nstates]) + } + fn from_vec(vec: Vec) -> Self { + faer::Col::from_vec(vec) + } + fn zeros(nstates: usize) -> Self { + Self::from_element(nstates, 0.0) + } + fn add_scalar_mut(&mut self, scalar: Self::T) { + self = self + scalar + } + fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { + // faer::linalg::matmul::matmul( + // self.as_mut(), + // self.as_ref(), + // x.as_ref(), + // Some(beta), + // alpha, + // faer::Parallelism::None, + // ); + self = self * faer::scale(beta) + x * faer::scale(alpha) + } + fn component_mul_assign(&mut self, other: &Self) { + // faer::linalg::matmul::matmul( + // self.as_mut(), + // self.as_ref(), + // other.as_ref(), + // None, + // 1.0, + // faer::Parallelism::None, + // ); + // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); + self = self * other + } +} + +impl VectorIndex for Vec { + fn zeros(len: IndexType) -> Self { + vec![0; len as usize] + } + fn len(&self) -> IndexType { + self.len() as IndexType + } +} + +impl VectorCommon for faer::Col { + type T = f64; +} +impl<'a> VectorCommon for faer::ColRef<'a, f64> { + type T = f64; +} +impl<'a> VectorCommon for faer::ColMut<'a, f64> { + type T = f64; +} + +impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { + type Owned = faer::Col; + fn abs(&self) -> faer::Col { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn into_owned(self) -> faer::Col { + self.to_owned() + } +} + +impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { + type Owned = faer::Col; + type View = faer::ColRef<'a, f64>; + fn abs(&self) -> faer::Col { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn copy_from(&mut self, other: &Self::Owned) { + self.copy_from(other); + } + fn copy_from_view(&mut self, other: &Self::View) { + self.copy_from(other); + } +} From e4ad22064ddc9721ab604da3b3f45897dd2eea75 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Sun, 24 Mar 2024 23:41:44 +0000 Subject: [PATCH 02/28] WIP: Vector implementation using faer --- src/vector/mod.rs | 2 +- src/vector/serial.rs | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index ff6e774b..01ad0e17 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -143,7 +143,7 @@ pub trait Vector: fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T); fn add_scalar_mut(&mut self, scalar: Self::T); fn component_mul_assign(&mut self, other: &Self); - // fn component_div_assign(&mut self, other: &Self); + fn component_div_assign(&mut self, other: &Self); fn filter_indices bool>(&self, f: F) -> Self::Index; fn filter(&self, indices: &Self::Index) -> Self { let mut result = Self::zeros(indices.len()); diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 9ba7f893..3c7c6147 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -92,9 +92,9 @@ impl Vector for DVector { fn axpy(&mut self, alpha: T, x: &Self, beta: T) { self.axpy(alpha, x, beta); } - // fn component_div_assign(&mut self, other: &Self) { - // self.component_div_assign(other); - // } + fn component_div_assign(&mut self, other: &Self) { + self.component_div_assign(other); + } fn component_mul_assign(&mut self, other: &Self) { self.component_mul_assign(other); } From 40c717e063cd6c7a4c7fda7ca383b441d2c22de2 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Sun, 24 Mar 2024 23:46:38 +0000 Subject: [PATCH 03/28] WIP: my approach is to have a second type inside VectorCommon and MatrixCommon to represent the type they are allowed to operate with, need more testing --- src/vector/mod.rs | 11 +++++++---- src/vector/serial.rs | 3 +++ src/vector/vector_faer.rs | 3 +++ 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 01ad0e17..eb848b98 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -19,6 +19,7 @@ pub trait VectorIndex: Sized + Index + Debug { pub trait VectorCommon: Sized + Debug { // + Display type T: Scalar; + type Op; } impl<'a, V> VectorCommon for &'a V @@ -26,6 +27,7 @@ where V: VectorCommon, { type T = V::T; + type Op = V::Op; } impl<'a, V> VectorCommon for &'a mut V @@ -33,6 +35,7 @@ where V: VectorCommon, { type T = V::T; + type Op = V::Op; } pub trait VectorOpsByValue: @@ -74,7 +77,7 @@ pub trait VectorViewMut<'a>: + VectorMutOpsByValue + for<'b> VectorMutOpsByValue<&'b Self::View> + for<'b> VectorMutOpsByValue<&'b Self::Owned> - + MulAssign + + MulAssign // + DivAssign + Index + IndexMut @@ -91,7 +94,7 @@ pub trait VectorView<'a>: + VectorOpsByValue + for<'b> VectorOpsByValue<&'b Self::Owned, Self::Owned> + for<'b> VectorOpsByValue<&'b Self, Self::Owned> - + Mul + + Mul // + Div + Index { @@ -105,8 +108,8 @@ pub trait Vector: + for<'b> VectorOpsByValue<&'b Self> + for<'a> VectorOpsByValue> + for<'a, 'b> VectorOpsByValue<&'b Self::View<'a>> - + Mul - + Div + + Mul + // + Div + VectorMutOpsByValue + for<'a> VectorMutOpsByValue> + for<'b> VectorMutOpsByValue<&'b Self> diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 3c7c6147..381daa16 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -15,14 +15,17 @@ impl VectorIndex for DVector { impl VectorCommon for DVector { type T = T; + type Op = T; } impl<'a, T: Scalar> VectorCommon for DVectorView<'a, T> { type T = T; + type Op = T; } impl<'a, T: Scalar> VectorCommon for DVectorViewMut<'a, T> { type T = T; + type Op = T; } impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index 9903cae9..d100d806 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -76,12 +76,15 @@ impl VectorIndex for Vec { impl VectorCommon for faer::Col { type T = f64; + type Op = faer::Scale; } impl<'a> VectorCommon for faer::ColRef<'a, f64> { type T = f64; + type Op = faer::Scale; } impl<'a> VectorCommon for faer::ColMut<'a, f64> { type T = f64; + type Op = faer::Scale; } impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { From f1b399cd6109a6918d3170e0e87972cd5928db38 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Mon, 25 Mar 2024 10:51:00 +0000 Subject: [PATCH 04/28] Revert "WIP: Vector implementation using faer" This reverts commit 98c42ac52c8093ce9c2f53eb26d78e7642cf04a7. --- src/vector/mod.rs | 2 +- src/vector/serial.rs | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index eb848b98..cc197837 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -146,7 +146,7 @@ pub trait Vector: fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T); fn add_scalar_mut(&mut self, scalar: Self::T); fn component_mul_assign(&mut self, other: &Self); - fn component_div_assign(&mut self, other: &Self); + // fn component_div_assign(&mut self, other: &Self); fn filter_indices bool>(&self, f: F) -> Self::Index; fn filter(&self, indices: &Self::Index) -> Self { let mut result = Self::zeros(indices.len()); diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 381daa16..f4bec415 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -95,9 +95,9 @@ impl Vector for DVector { fn axpy(&mut self, alpha: T, x: &Self, beta: T) { self.axpy(alpha, x, beta); } - fn component_div_assign(&mut self, other: &Self) { - self.component_div_assign(other); - } + // fn component_div_assign(&mut self, other: &Self) { + // self.component_div_assign(other); + // } fn component_mul_assign(&mut self, other: &Self) { self.component_mul_assign(other); } From e8fc4789a35021452aae7aeaace61294ea6bb621 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Mon, 25 Mar 2024 10:52:36 +0000 Subject: [PATCH 05/28] Revert "WIP: my approach is to have a second type inside VectorCommon and MatrixCommon to represent the type they are allowed to operate with, need more testing" This reverts commit e8c8998f0bedc14b4419e8dd4ee3eff9ffae871a. --- src/vector/mod.rs | 11 ++++------- src/vector/serial.rs | 3 --- src/vector/vector_faer.rs | 3 --- 3 files changed, 4 insertions(+), 13 deletions(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index cc197837..ff6e774b 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -19,7 +19,6 @@ pub trait VectorIndex: Sized + Index + Debug { pub trait VectorCommon: Sized + Debug { // + Display type T: Scalar; - type Op; } impl<'a, V> VectorCommon for &'a V @@ -27,7 +26,6 @@ where V: VectorCommon, { type T = V::T; - type Op = V::Op; } impl<'a, V> VectorCommon for &'a mut V @@ -35,7 +33,6 @@ where V: VectorCommon, { type T = V::T; - type Op = V::Op; } pub trait VectorOpsByValue: @@ -77,7 +74,7 @@ pub trait VectorViewMut<'a>: + VectorMutOpsByValue + for<'b> VectorMutOpsByValue<&'b Self::View> + for<'b> VectorMutOpsByValue<&'b Self::Owned> - + MulAssign + + MulAssign // + DivAssign + Index + IndexMut @@ -94,7 +91,7 @@ pub trait VectorView<'a>: + VectorOpsByValue + for<'b> VectorOpsByValue<&'b Self::Owned, Self::Owned> + for<'b> VectorOpsByValue<&'b Self, Self::Owned> - + Mul + + Mul // + Div + Index { @@ -108,8 +105,8 @@ pub trait Vector: + for<'b> VectorOpsByValue<&'b Self> + for<'a> VectorOpsByValue> + for<'a, 'b> VectorOpsByValue<&'b Self::View<'a>> - + Mul - // + Div + + Mul + + Div + VectorMutOpsByValue + for<'a> VectorMutOpsByValue> + for<'b> VectorMutOpsByValue<&'b Self> diff --git a/src/vector/serial.rs b/src/vector/serial.rs index f4bec415..9ba7f893 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -15,17 +15,14 @@ impl VectorIndex for DVector { impl VectorCommon for DVector { type T = T; - type Op = T; } impl<'a, T: Scalar> VectorCommon for DVectorView<'a, T> { type T = T; - type Op = T; } impl<'a, T: Scalar> VectorCommon for DVectorViewMut<'a, T> { type T = T; - type Op = T; } impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index d100d806..9903cae9 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -76,15 +76,12 @@ impl VectorIndex for Vec { impl VectorCommon for faer::Col { type T = f64; - type Op = faer::Scale; } impl<'a> VectorCommon for faer::ColRef<'a, f64> { type T = f64; - type Op = faer::Scale; } impl<'a> VectorCommon for faer::ColMut<'a, f64> { type T = f64; - type Op = faer::Scale; } impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { From 5b3ba6c66034606c2475a04693c7b4d23c538088 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Mon, 25 Mar 2024 15:10:34 +0000 Subject: [PATCH 06/28] WIP: going back to a point it compiles, keep trying --- src/vector/mod.rs | 2 +- src/vector/serial.rs | 6 +- src/vector/vector_faer.rs | 206 +++++++++++++++++++------------------- 3 files changed, 107 insertions(+), 107 deletions(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index ff6e774b..01ad0e17 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -143,7 +143,7 @@ pub trait Vector: fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T); fn add_scalar_mut(&mut self, scalar: Self::T); fn component_mul_assign(&mut self, other: &Self); - // fn component_div_assign(&mut self, other: &Self); + fn component_div_assign(&mut self, other: &Self); fn filter_indices bool>(&self, f: F) -> Self::Index; fn filter(&self, indices: &Self::Index) -> Self { let mut result = Self::zeros(indices.len()); diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 9ba7f893..3c7c6147 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -92,9 +92,9 @@ impl Vector for DVector { fn axpy(&mut self, alpha: T, x: &Self, beta: T) { self.axpy(alpha, x, beta); } - // fn component_div_assign(&mut self, other: &Self) { - // self.component_div_assign(other); - // } + fn component_div_assign(&mut self, other: &Self) { + self.component_div_assign(other); + } fn component_mul_assign(&mut self, other: &Self) { self.component_mul_assign(other); } diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index 9903cae9..00aa195e 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -1,109 +1,109 @@ -use faer::{unzipped, zipped}; +// use faer::{unzipped, zipped}; -use crate::IndexType; +// use crate::IndexType; -use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; -impl Vector for faer::Col { - type View<'a> = faer::ColRef<'a, f64>; - type ViewMut<'a> = faer::ColMut<'a, f64>; - type Index = Vec; - fn len(&self) -> IndexType { - self.nrows() - } - fn norm(&self) -> f64 { - self.norm_l2() - } - fn abs(&self) -> Self { - zipped!(self).map(|unzipped!(xi)| xi.abs()) - } - fn as_view(&self) -> Self::View<'_> { - self.as_ref() - } - fn as_view_mut(&mut self) -> Self::ViewMut<'_> { - self.as_mut() - } - fn copy_from(&mut self, other: &Self) { - self.copy_from(other) - } - fn copy_from_view(&mut self, other: &Self::View<'_>) { - self.copy_from(other) - } - fn from_element(nstates: usize, value: Self::T) -> Self { - faer::Col::from_vec(vec![value; nstates]) - } - fn from_vec(vec: Vec) -> Self { - faer::Col::from_vec(vec) - } - fn zeros(nstates: usize) -> Self { - Self::from_element(nstates, 0.0) - } - fn add_scalar_mut(&mut self, scalar: Self::T) { - self = self + scalar - } - fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { - // faer::linalg::matmul::matmul( - // self.as_mut(), - // self.as_ref(), - // x.as_ref(), - // Some(beta), - // alpha, - // faer::Parallelism::None, - // ); - self = self * faer::scale(beta) + x * faer::scale(alpha) - } - fn component_mul_assign(&mut self, other: &Self) { - // faer::linalg::matmul::matmul( - // self.as_mut(), - // self.as_ref(), - // other.as_ref(), - // None, - // 1.0, - // faer::Parallelism::None, - // ); - // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); - self = self * other - } -} +// use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; +// impl Vector for faer::Col { +// type View<'a> = faer::ColRef<'a, f64>; +// type ViewMut<'a> = faer::ColMut<'a, f64>; +// type Index = Vec; +// fn len(&self) -> IndexType { +// self.nrows() +// } +// fn norm(&self) -> f64 { +// self.norm_l2() +// } +// fn abs(&self) -> Self { +// zipped!(self).map(|unzipped!(xi)| xi.abs()) +// } +// fn as_view(&self) -> Self::View<'_> { +// self.as_ref() +// } +// fn as_view_mut(&mut self) -> Self::ViewMut<'_> { +// self.as_mut() +// } +// fn copy_from(&mut self, other: &Self) { +// self.copy_from(other) +// } +// fn copy_from_view(&mut self, other: &Self::View<'_>) { +// self.copy_from(other) +// } +// fn from_element(nstates: usize, value: Self::T) -> Self { +// faer::Col::from_vec(vec![value; nstates]) +// } +// fn from_vec(vec: Vec) -> Self { +// faer::Col::from_vec(vec) +// } +// fn zeros(nstates: usize) -> Self { +// Self::from_element(nstates, 0.0) +// } +// fn add_scalar_mut(&mut self, scalar: Self::T) { +// self = self + scalar +// } +// fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { +// // faer::linalg::matmul::matmul( +// // self.as_mut(), +// // self.as_ref(), +// // x.as_ref(), +// // Some(beta), +// // alpha, +// // faer::Parallelism::None, +// // ); +// self = self * faer::scale(beta) + x * faer::scale(alpha) +// } +// fn component_mul_assign(&mut self, other: &Self) { +// // faer::linalg::matmul::matmul( +// // self.as_mut(), +// // self.as_ref(), +// // other.as_ref(), +// // None, +// // 1.0, +// // faer::Parallelism::None, +// // ); +// // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); +// self = self * other +// } +// } -impl VectorIndex for Vec { - fn zeros(len: IndexType) -> Self { - vec![0; len as usize] - } - fn len(&self) -> IndexType { - self.len() as IndexType - } -} +// impl VectorIndex for Vec { +// fn zeros(len: IndexType) -> Self { +// vec![0; len as usize] +// } +// fn len(&self) -> IndexType { +// self.len() as IndexType +// } +// } -impl VectorCommon for faer::Col { - type T = f64; -} -impl<'a> VectorCommon for faer::ColRef<'a, f64> { - type T = f64; -} -impl<'a> VectorCommon for faer::ColMut<'a, f64> { - type T = f64; -} +// impl VectorCommon for faer::Col { +// type T = f64; +// } +// impl<'a> VectorCommon for faer::ColRef<'a, f64> { +// type T = f64; +// } +// impl<'a> VectorCommon for faer::ColMut<'a, f64> { +// type T = f64; +// } -impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { - type Owned = faer::Col; - fn abs(&self) -> faer::Col { - zipped!(self).map(|unzipped!(xi)| xi.abs()) - } - fn into_owned(self) -> faer::Col { - self.to_owned() - } -} +// impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { +// type Owned = faer::Col; +// fn abs(&self) -> faer::Col { +// zipped!(self).map(|unzipped!(xi)| xi.abs()) +// } +// fn into_owned(self) -> faer::Col { +// self.to_owned() +// } +// } -impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { - type Owned = faer::Col; - type View = faer::ColRef<'a, f64>; - fn abs(&self) -> faer::Col { - zipped!(self).map(|unzipped!(xi)| xi.abs()) - } - fn copy_from(&mut self, other: &Self::Owned) { - self.copy_from(other); - } - fn copy_from_view(&mut self, other: &Self::View) { - self.copy_from(other); - } -} +// impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { +// type Owned = faer::Col; +// type View = faer::ColRef<'a, f64>; +// fn abs(&self) -> faer::Col { +// zipped!(self).map(|unzipped!(xi)| xi.abs()) +// } +// fn copy_from(&mut self, other: &Self::Owned) { +// self.copy_from(other); +// } +// fn copy_from_view(&mut self, other: &Self::View) { +// self.copy_from(other); +// } +// } From 0059fdc2db77c588d5f9a803e74b5676dc675f18 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 14:56:26 +0000 Subject: [PATCH 07/28] preparing project structure --- src/lib.rs | 42 +++---------------------------------- src/matrix/dense_serial.rs | 3 +++ src/matrix/mod.rs | 3 +++ src/matrix/sparse_serial.rs | 1 + src/scalar/mod.rs | 38 +++++++++++++++++++++++++++++++++ src/vector/mod.rs | 10 ++++----- src/vector/serial.rs | 3 +++ src/vector/vector_faer.rs | 3 +++ 8 files changed, 59 insertions(+), 44 deletions(-) create mode 100644 src/scalar/mod.rs diff --git a/src/lib.rs b/src/lib.rs index 91d5ff84..eaba57fd 100644 --- a/src/lib.rs +++ b/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 - + Display - + SimdRealField - + ComplexField - + Copy - + ClosedSub - + From - + ClosedMul - + ClosedDiv - + ClosedAdd - + Signed - + PartialOrd - + Pow - + Pow -{ - 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 op::{LinearOp, NonLinearOp}; +use scalar::{IndexType, Scalar}; use solver::SolverProblem; use vector::{Vector, VectorIndex, VectorRef, VectorView, VectorViewMut}; diff --git a/src/matrix/dense_serial.rs b/src/matrix/dense_serial.rs index 3fd2a08d..e265761f 100644 --- a/src/matrix/dense_serial.rs +++ b/src/matrix/dense_serial.rs @@ -8,6 +8,7 @@ use super::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut}; impl<'a, T: Scalar> MatrixCommon for DMatrixViewMut<'a, T> { type V = DVector; type T = T; + type O = T; fn ncols(&self) -> IndexType { self.ncols() @@ -31,6 +32,7 @@ impl<'a, T: Scalar> MatrixViewMut<'a> for DMatrixViewMut<'a, T> { impl<'a, T: Scalar> MatrixCommon for DMatrixView<'a, T> { type V = DVector; type T = T; + type O = T; fn ncols(&self) -> IndexType { self.ncols() } @@ -46,6 +48,7 @@ impl<'a, T: Scalar> MatrixView<'a> for DMatrixView<'a, T> { impl MatrixCommon for DMatrix { type V = DVector; type T = T; + type O = T; fn ncols(&self) -> IndexType { self.ncols() diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs index 95bf2d91..bfc2f2b8 100644 --- a/src/matrix/mod.rs +++ b/src/matrix/mod.rs @@ -11,6 +11,7 @@ mod sparse_serial; pub trait MatrixCommon: Sized + Debug { type V: Vector; type T: Scalar; + type O; /// Get the number of columns of the matrix fn nrows(&self) -> IndexType; @@ -25,6 +26,7 @@ where { type T = M::T; type V = M::V; + type O = M::O; fn ncols(&self) -> IndexType { M::ncols(self) @@ -40,6 +42,7 @@ where { type T = M::T; type V = M::V; + type O = M::O; fn ncols(&self) -> IndexType { M::ncols(self) } diff --git a/src/matrix/sparse_serial.rs b/src/matrix/sparse_serial.rs index ab17c992..2a705c20 100644 --- a/src/matrix/sparse_serial.rs +++ b/src/matrix/sparse_serial.rs @@ -9,6 +9,7 @@ use super::{Matrix, MatrixCommon}; impl MatrixCommon for CscMatrix { type V = DVector; type T = T; + type O = T; fn ncols(&self) -> IndexType { self.ncols() diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs new file mode 100644 index 00000000..01e5a009 --- /dev/null +++ b/src/scalar/mod.rs @@ -0,0 +1,38 @@ +use std::fmt::Display; + +use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField}; +use num_traits::{Pow, Signed}; + +pub trait Scalar: + nalgebra::Scalar + + From + + Display + + SimdRealField + + ComplexField + + Copy + + ClosedSub + + From + + ClosedMul + + ClosedDiv + + ClosedAdd + + Signed + + PartialOrd + + Pow + + Pow +{ + const EPSILON: Self; + const INFINITY: Self; + const NAN: Self; + fn is_nan(self) -> bool; +} + +pub 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() + } +} diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 01ad0e17..863d3a5d 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -1,9 +1,8 @@ +use crate::{IndexType, Scalar}; use num_traits::Zero; use std::fmt::{Debug, Display}; use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; -use crate::{IndexType, Scalar}; - mod serial; mod vector_faer; @@ -19,6 +18,7 @@ pub trait VectorIndex: Sized + Index + Debug { pub trait VectorCommon: Sized + Debug { // + Display type T: Scalar; + type O; } impl<'a, V> VectorCommon for &'a V @@ -26,6 +26,7 @@ where V: VectorCommon, { type T = V::T; + type O = V::O; } impl<'a, V> VectorCommon for &'a mut V @@ -33,6 +34,7 @@ where V: VectorCommon, { type T = V::T; + type O = V::O; } pub trait VectorOpsByValue: @@ -55,7 +57,6 @@ pub trait VectorRef: + for<'a> VectorOpsByValue, V> + for<'a, 'b> VectorOpsByValue<&'a V::View<'b>, V> + Mul - + Div { } @@ -64,8 +65,7 @@ impl VectorRef for RefT where + for<'a> VectorOpsByValue<&'a V, V> + for<'a> VectorOpsByValue, V> + for<'a, 'b> VectorOpsByValue<&'a V::View<'b>, V> - + Mul - + Div + + Mul // + Div { } diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 3c7c6147..8322acc2 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -15,14 +15,17 @@ impl VectorIndex for DVector { impl VectorCommon for DVector { type T = T; + type O = T; } impl<'a, T: Scalar> VectorCommon for DVectorView<'a, T> { type T = T; + type O = T; } impl<'a, T: Scalar> VectorCommon for DVectorViewMut<'a, T> { type T = T; + type O = T; } impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index 00aa195e..06bdc3a2 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -76,12 +76,15 @@ // impl VectorCommon for faer::Col { // type T = f64; +// type O = faer::Scale; // } // impl<'a> VectorCommon for faer::ColRef<'a, f64> { // type T = f64; +// type O = faer::Scale; // } // impl<'a> VectorCommon for faer::ColMut<'a, f64> { // type T = f64; +// type O = faer::Scale; // } // impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { From fcc09b63787a79f25d96d637befef7d66b8c61f4 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 16:12:19 +0000 Subject: [PATCH 08/28] Started implementing Scale and some of its relevant traits --- src/scalar/mod.rs | 81 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 80 insertions(+), 1 deletion(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index 01e5a009..ae6d2d6a 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -1,4 +1,7 @@ -use std::fmt::Display; +use std::{ + fmt::Display, + ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign}, +}; use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField}; use num_traits::{Pow, Signed}; @@ -36,3 +39,79 @@ impl Scalar for f64 { self.is_nan() } } + +#[derive(Copy, Clone, Debug)] +pub struct Scale(pub E); + +impl Scale { + #[inline] + pub fn value(self) -> E { + self.0 + } +} + +#[inline] +pub fn scale(value: E) -> Scale { + Scale(value) +} + +//TODO: Is it possible for us to need RhsE != LhsE? +impl Mul> for Scale { + type Output = Scale; + + #[inline] + fn mul(self, rhs: Scale) -> Self::Output { + Scale(self.0 * rhs.0) + } +} + +impl Add> for Scale { + type Output = Scale; + + #[inline] + fn add(self, rhs: Scale) -> Self::Output { + Scale(self.0 + rhs.0) + } +} + +impl Sub> for Scale { + type Output = Scale; + + #[inline] + fn sub(self, rhs: Scale) -> Self::Output { + Scale(self.0 - rhs.0) + } +} + +impl MulAssign> for Scale { + #[inline] + fn mul_assign(&mut self, rhs: Scale) { + self.0 = self.0 * rhs.0 + } +} + +impl AddAssign> for Scale { + #[inline] + fn add_assign(&mut self, rhs: Scale) { + self.0 = self.0 + rhs.0 + } +} + +impl SubAssign> for Scale { + #[inline] + fn sub_assign(&mut self, rhs: Scale) { + self.0 = self.0 - rhs.0 + } +} + +impl PartialEq for Scale { + #[inline] + fn eq(&self, rhs: &Self) -> bool { + self.0 == rhs.0 + } +} + +#[test] +fn test_scale() { + assert_eq!(scale(2.0) * scale(3.0), scale(6.0)); +} From e3d48ab70dc6198233c438eb257ff8ef93e4e9fc Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 16:57:53 +0000 Subject: [PATCH 09/28] mul trait --- src/scalar/mod.rs | 11 +++++++++++ src/vector/mod.rs | 1 + src/vector/serial.rs | 3 +++ 3 files changed, 15 insertions(+) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index ae6d2d6a..faa13680 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -6,6 +6,8 @@ use std::{ use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField}; use num_traits::{Pow, Signed}; +use crate::vector::{Vector, VectorCommon, VectorRef, VectorView}; + pub trait Scalar: nalgebra::Scalar + From @@ -111,6 +113,15 @@ impl PartialEq for Scale { } } +// I have to implement Mul between Scale and VectorView +impl<'a, V: VectorView<'static>> Mul for Scale { + type Output = V::Owned; + #[inline] + fn mul(self, rhs: V) -> Self::Output { + rhs.scalar_mul(self.0.into()) + } +} + #[test] fn test_scale() { assert_eq!(scale(2.0) * scale(3.0), scale(6.0)); diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 863d3a5d..22bc9ec0 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -98,6 +98,7 @@ pub trait VectorView<'a>: type Owned; fn abs(&self) -> Self::Owned; fn into_owned(self) -> Self::Owned; + fn scalar_mul(&self, rhs: Self::T) -> Self::Owned; } pub trait Vector: diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 8322acc2..5a0500b1 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -36,6 +36,9 @@ impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { fn into_owned(self) -> Self::Owned { self.into_owned() } + fn scalar_mul(&self, rhs: Self::T) -> Self::Owned { + self * rhs + } } impl<'a, T: Scalar> VectorViewMut<'a> for DVectorViewMut<'a, T> { From 6ee9dcc13ce3625dca74c47f4597e788a54a27e9 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 17:01:45 +0000 Subject: [PATCH 10/28] generic over the E Scale --- src/scalar/mod.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index faa13680..6b89a92c 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -114,11 +114,11 @@ impl PartialEq for Scale { } // I have to implement Mul between Scale and VectorView -impl<'a, V: VectorView<'static>> Mul for Scale { +impl<'a, V: VectorView<'static, T = E>, E: Scalar> Mul for Scale { type Output = V::Owned; #[inline] fn mul(self, rhs: V) -> Self::Output { - rhs.scalar_mul(self.0.into()) + rhs.scalar_mul(self.0) } } From 2672cfaef14702c48b8a3d8dc48610629220338a Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 17:50:40 +0000 Subject: [PATCH 11/28] initial version of the abstraction is working --- src/ode_solver/bdf.rs | 10 +++++----- src/op/ode.rs | 5 +++-- src/scalar/mod.rs | 16 ++++++++++++++-- src/vector/mod.rs | 6 ++---- src/vector/serial.rs | 13 +++++++++---- 5 files changed, 33 insertions(+), 17 deletions(-) diff --git a/src/ode_solver/bdf.rs b/src/ode_solver/bdf.rs index 0eb437e3..7eca934f 100644 --- a/src/ode_solver/bdf.rs +++ b/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,7 +245,7 @@ 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 } @@ -253,7 +253,7 @@ where fn problem(&self) -> Option<&OdeSolverProblem> { self.ode_problem.as_ref() } - + fn set_problem(&mut self, state: &mut OdeSolverState, problem: &OdeSolverProblem) { self.ode_problem = Some(problem.clone()); let nstates = problem.eqn.nstates(); @@ -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 { diff --git a/src/op/ode.rs b/src/op/ode.rs index 860ca6ac..06244b12 100644 --- a/src/op/ode.rs +++ b/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,9 +93,9 @@ impl BdfCallable { 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]; diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index 6b89a92c..fad4a269 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -6,7 +6,7 @@ use std::{ use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField}; use num_traits::{Pow, Signed}; -use crate::vector::{Vector, VectorCommon, VectorRef, VectorView}; +use crate::vector::{Vector, VectorCommon, VectorRef, VectorView, VectorViewMut}; pub trait Scalar: nalgebra::Scalar @@ -113,7 +113,6 @@ impl PartialEq for Scale { } } -// I have to implement Mul between Scale and VectorView impl<'a, V: VectorView<'static, T = E>, E: Scalar> Mul for Scale { type Output = V::Owned; #[inline] @@ -122,6 +121,19 @@ impl<'a, V: VectorView<'static, T = E>, E: Scalar> Mul for Scale { } } +//TODO: Not sure why it is now working +// impl<'a, E, V> Mul> for V +// where +// V: VectorView<'static, T = E>, +// E: Scalar, +// { +// type Output = V::Owned; +// #[inline] +// fn mul(self, rhs: Scale) -> Self::Output { +// self.scalar_mul(rhs.0) +// } +// } + #[test] fn test_scale() { assert_eq!(scale(2.0) * scale(3.0), scale(6.0)); diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 22bc9ec0..9202cd7b 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -1,3 +1,4 @@ +use crate::scalar::Scale; use crate::{IndexType, Scalar}; use num_traits::Zero; use std::fmt::{Debug, Display}; @@ -18,7 +19,6 @@ pub trait VectorIndex: Sized + Index + Debug { pub trait VectorCommon: Sized + Debug { // + Display type T: Scalar; - type O; } impl<'a, V> VectorCommon for &'a V @@ -26,7 +26,6 @@ where V: VectorCommon, { type T = V::T; - type O = V::O; } impl<'a, V> VectorCommon for &'a mut V @@ -34,7 +33,6 @@ where V: VectorCommon, { type T = V::T; - type O = V::O; } pub trait VectorOpsByValue: @@ -91,7 +89,7 @@ pub trait VectorView<'a>: + VectorOpsByValue + for<'b> VectorOpsByValue<&'b Self::Owned, Self::Owned> + for<'b> VectorOpsByValue<&'b Self, Self::Owned> - + Mul + + Mul, Output = Self::Owned> // + Div + Index { diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 5a0500b1..a3057de2 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -1,6 +1,8 @@ +use std::ops::Mul; + use nalgebra::{DVector, DVectorView, DVectorViewMut}; -use crate::{IndexType, Scalar}; +use crate::{scalar::Scale, IndexType, Scalar}; use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; @@ -15,17 +17,14 @@ impl VectorIndex for DVector { impl VectorCommon for DVector { type T = T; - type O = T; } impl<'a, T: Scalar> VectorCommon for DVectorView<'a, T> { type T = T; - type O = T; } impl<'a, T: Scalar> VectorCommon for DVectorViewMut<'a, T> { type T = T; - type O = T; } impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { @@ -40,6 +39,12 @@ impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { self * rhs } } +impl<'a, T: Scalar> Mul> for DVectorView<'a, T> { + type Output = DVector; + fn mul(self, rhs: Scale) -> Self::Output { + self * rhs + } +} impl<'a, T: Scalar> VectorViewMut<'a> for DVectorViewMut<'a, T> { type Owned = DVector; From 12abb965b5b8e0576886e5b32afdc2332c57397c Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 18:26:18 +0000 Subject: [PATCH 12/28] initial implementation --- src/nonlinear_solver/mod.rs | 6 ++-- src/ode_solver/bdf.rs | 4 +-- src/ode_solver/mod.rs | 3 +- src/ode_solver/test_models/dydt_y2.rs | 5 ++-- .../test_models/exponential_decay.rs | 7 +++-- .../exponential_decay_with_algebraic.rs | 7 +++-- src/ode_solver/test_models/gaussian_decay.rs | 7 +++-- src/op/mod.rs | 6 ++-- src/op/ode.rs | 2 +- src/vector/mod.rs | 8 +++--- src/vector/serial.rs | 28 +++++++++++++++++-- src/vector/vector_faer.rs | 27 +++++++++++++++--- 12 files changed, 80 insertions(+), 30 deletions(-) diff --git a/src/nonlinear_solver/mod.rs b/src/nonlinear_solver/mod.rs index d05a8db7..7f56cc58 100644 --- a/src/nonlinear_solver/mod.rs +++ b/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 { pub x0: V, @@ -84,7 +84,7 @@ impl Convergence { } } 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]); } diff --git a/src/ode_solver/bdf.rs b/src/ode_solver/bdf.rs index 7eca934f..c41bea24 100644 --- a/src/ode_solver/bdf.rs +++ b/src/ode_solver/bdf.rs @@ -290,7 +290,7 @@ where // update initial step size based on function let mut scale = state.y.abs(); - scale *= problem.rtol; + scale *= crate::scalar::scale(problem.rtol); scale += problem.atol.as_ref(); let f0 = problem.eqn.rhs(state.t, &state.y); @@ -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(); } diff --git a/src/ode_solver/mod.rs b/src/ode_solver/mod.rs index 9fc3450f..2b50a437 100644 --- a/src/ode_solver/mod.rs +++ b/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]); } diff --git a/src/ode_solver/test_models/dydt_y2.rs b/src/ode_solver/test_models/dydt_y2.rs index ee60be3f..bfca49cc 100644 --- a/src/ode_solver/test_models/dydt_y2.rs +++ b/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(x: &M::V, _p: &M::V, _t: M::T, y: &mut M::V) { fn rhs_jac(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( @@ -42,7 +43,7 @@ pub fn dydt_y2_problem( 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); diff --git a/src/ode_solver/test_models/exponential_decay.rs b/src/ode_solver/test_models/exponential_decay.rs index 0634ab0f..a00e3a41 100644 --- a/src/ode_solver/test_models/exponential_decay.rs +++ b/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(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( y: &mut M::V, ) { y.copy_from(v); - y.mul_assign(-p[0]); + y.mul_assign(scale(-p[0])); } fn exponential_decay_init(_p: &M::V, _t: M::T) -> M::V { @@ -49,7 +50,7 @@ pub fn exponential_decay_problem( 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) diff --git a/src/ode_solver/test_models/exponential_decay_with_algebraic.rs b/src/ode_solver/test_models/exponential_decay_with_algebraic.rs index 46a8c840..4301e732 100644 --- a/src/ode_solver/test_models/exponential_decay_with_algebraic.rs +++ b/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( 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( 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( 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) diff --git a/src/ode_solver/test_models/gaussian_decay.rs b/src/ode_solver/test_models/gaussian_decay.rs index 0eaf0f82..a77316c0 100644 --- a/src/ode_solver/test_models/gaussian_decay.rs +++ b/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(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(_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( @@ -41,7 +42,7 @@ pub fn gaussian_decay_problem( 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); diff --git a/src/op/mod.rs b/src/op/mod.rs index f941e8fc..184b078e 100644 --- a/src/op/mod.rs +++ b/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 { diff --git a/src/op/ode.rs b/src/op/ode.rs index 06244b12..788243ec 100644 --- a/src/op/ode.rs +++ b/src/op/ode.rs @@ -97,7 +97,7 @@ impl BdfCallable { for (i, &gamma_i) in gamma.iter().enumerate().take(order + 1).skip(2) { new_psi_neg_y0 += diff.column(i) * scale(gamma_i) } - new_psi_neg_y0 *= alpha[order]; + new_psi_neg_y0 = new_psi_neg_y0 * scale(alpha[order]); // now negate y0 new_psi_neg_y0.sub_assign(y0); diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 9202cd7b..185f70c1 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -104,14 +104,14 @@ pub trait Vector: + for<'b> VectorOpsByValue<&'b Self> + for<'a> VectorOpsByValue> + for<'a, 'b> VectorOpsByValue<&'b Self::View<'a>> - + Mul - + Div + + Mul, Output = Self> + + Div, Output = Self> + VectorMutOpsByValue + for<'a> VectorMutOpsByValue> + for<'b> VectorMutOpsByValue<&'b Self> + for<'a, 'b> VectorMutOpsByValue<&'b Self::View<'a>> - + MulAssign - + DivAssign + + MulAssign> + // + DivAssign + Index + IndexMut + Clone diff --git a/src/vector/serial.rs b/src/vector/serial.rs index a3057de2..8f7eb6b5 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -1,4 +1,4 @@ -use std::ops::Mul; +use std::ops::{Div, Mul, MulAssign}; use nalgebra::{DVector, DVectorView, DVectorViewMut}; @@ -42,7 +42,24 @@ impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { impl<'a, T: Scalar> Mul> for DVectorView<'a, T> { type Output = DVector; fn mul(self, rhs: Scale) -> Self::Output { - self * rhs + self * rhs.value() + } +} +impl<'a, T: Scalar> Mul> for DVectorViewMut<'a, T> { + type Output = DVector; + fn mul(self, rhs: Scale) -> Self::Output { + self * rhs.value() + } +} +impl Mul> for DVector { + type Output = DVector; + fn mul(self, rhs: Scale) -> Self::Output { + self * rhs.value() + } +} +impl MulAssign> for DVector { + fn mul_assign(&mut self, rhs: Scale) { + *self = &*self * rhs.value(); } } @@ -60,6 +77,13 @@ impl<'a, T: Scalar> VectorViewMut<'a> for DVectorViewMut<'a, T> { } } +impl Div> for DVector { + type Output = DVector; + fn div(self, rhs: Scale) -> Self::Output { + self / rhs.value() + } +} + impl Vector for DVector { type View<'a> = DVectorView<'a, T>; type ViewMut<'a> = DVectorViewMut<'a, T>; diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index 06bdc3a2..8ff4c3a0 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -1,8 +1,30 @@ +// use std::ops::{Mul, MulAssign}; + // use faer::{unzipped, zipped}; -// use crate::IndexType; +// use crate::{scalar::Scalar, IndexType}; // use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; +// impl<'a, T: Scalar> Mul> for faer::ColRef<'a, f64> { +// type Output = faer::Col; +// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +// self * faer::scale(rhs.into()) +// } +// } + +// impl<'a, T: Scalar> Mul> for faer::ColMut<'a, f64> { +// type Output = faer::Col; +// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +// self * faer::scale(rhs) +// } +// } + +// impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { +// fn mul_assign(&mut self, rhs: crate::scalar::Scale) { +// self = self * faer::scale(rhs) +// } +// } + // impl Vector for faer::Col { // type View<'a> = faer::ColRef<'a, f64>; // type ViewMut<'a> = faer::ColMut<'a, f64>; @@ -76,15 +98,12 @@ // impl VectorCommon for faer::Col { // type T = f64; -// type O = faer::Scale; // } // impl<'a> VectorCommon for faer::ColRef<'a, f64> { // type T = f64; -// type O = faer::Scale; // } // impl<'a> VectorCommon for faer::ColMut<'a, f64> { // type T = f64; -// type O = faer::Scale; // } // impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { From e00d8d496524c538b4ffc1affb7a5a5ebb4638f6 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 18:47:23 +0000 Subject: [PATCH 13/28] some missing traits, I think the migration is almost done --- src/vector/mod.rs | 2 +- src/vector/serial.rs | 10 ++++++++++ src/vector/vector_faer.rs | 10 ++++++++++ 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 185f70c1..22f891a9 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -72,7 +72,7 @@ pub trait VectorViewMut<'a>: + VectorMutOpsByValue + for<'b> VectorMutOpsByValue<&'b Self::View> + for<'b> VectorMutOpsByValue<&'b Self::Owned> - + MulAssign + + MulAssign> // + DivAssign + Index + IndexMut diff --git a/src/vector/serial.rs b/src/vector/serial.rs index 8f7eb6b5..de3370a2 100644 --- a/src/vector/serial.rs +++ b/src/vector/serial.rs @@ -57,11 +57,21 @@ impl Mul> for DVector { self * rhs.value() } } + impl MulAssign> for DVector { fn mul_assign(&mut self, rhs: Scale) { *self = &*self * rhs.value(); } } +impl<'a, T: Scalar> MulAssign> for DVectorViewMut<'a, T> { + fn mul_assign(&mut self, rhs: Scale) { + //I'm doing this horrendous hack because + //somehow I'm missing which trait to implement so I can multiply + //a &mut DVector by a scalar + let a = &*self * rhs.value(); + self.copy_from(&a); + } +} impl<'a, T: Scalar> VectorViewMut<'a> for DVectorViewMut<'a, T> { type Owned = DVector; diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index 8ff4c3a0..f5835852 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -25,6 +25,13 @@ // } // } +// impl<'a, T: Scalar> Mul> for faer::Col { +// type Output = faer::Col; +// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +// self * faer::scale(rhs) +// } +// } + // impl Vector for faer::Col { // type View<'a> = faer::ColRef<'a, f64>; // type ViewMut<'a> = faer::ColMut<'a, f64>; @@ -114,6 +121,9 @@ // fn into_owned(self) -> faer::Col { // self.to_owned() // } +// fn scalar_mul(&self, rhs: Self::T) -> Self::Owned { +// self * faer::scale(rhs) +// } // } // impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { From 61d2e789f82473f038b1107111422cb2e67d0f61 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 21:25:14 +0000 Subject: [PATCH 14/28] Faer implementation requires to enable the faer feature --- Cargo.toml | 2 +- src/vector/mod.rs | 1 + src/vector/vector_faer.rs | 260 +++++++++++++++++++------------------- 3 files changed, 132 insertions(+), 131 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 0c7a058c..0105f522 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,7 @@ license = "MIT" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [features] -default = ["faer"] +default = ["nalgebra"] faer = [] nalgebra = [] diffsl = [] diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 22f891a9..893774c1 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -5,6 +5,7 @@ use std::fmt::{Debug, Display}; use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; mod serial; +#[cfg(feature = "faer")] mod vector_faer; pub trait VectorIndex: Sized + Index + Debug { diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index f5835852..e654978b 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -1,141 +1,141 @@ -// use std::ops::{Mul, MulAssign}; +use std::ops::{Mul, MulAssign}; -// use faer::{unzipped, zipped}; +use faer::{unzipped, zipped}; -// use crate::{scalar::Scalar, IndexType}; +use crate::{scalar::Scalar, IndexType}; -// use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; -// impl<'a, T: Scalar> Mul> for faer::ColRef<'a, f64> { -// type Output = faer::Col; -// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { -// self * faer::scale(rhs.into()) -// } -// } +use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; +impl<'a, T: Scalar> Mul> for faer::ColRef<'a, f64> { + type Output = faer::Col; + fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { + self * faer::scale(rhs.into()) + } +} -// impl<'a, T: Scalar> Mul> for faer::ColMut<'a, f64> { -// type Output = faer::Col; -// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { -// self * faer::scale(rhs) -// } -// } +impl<'a, T: Scalar> Mul> for faer::ColMut<'a, f64> { + type Output = faer::Col; + fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { + self * faer::scale(rhs) + } +} -// impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { -// fn mul_assign(&mut self, rhs: crate::scalar::Scale) { -// self = self * faer::scale(rhs) -// } -// } +impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { + fn mul_assign(&mut self, rhs: crate::scalar::Scale) { + self = self * faer::scale(rhs) + } +} -// impl<'a, T: Scalar> Mul> for faer::Col { -// type Output = faer::Col; -// fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { -// self * faer::scale(rhs) -// } -// } +impl<'a, T: Scalar> Mul> for faer::Col { + type Output = faer::Col; + fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { + self * faer::scale(rhs) + } +} -// impl Vector for faer::Col { -// type View<'a> = faer::ColRef<'a, f64>; -// type ViewMut<'a> = faer::ColMut<'a, f64>; -// type Index = Vec; -// fn len(&self) -> IndexType { -// self.nrows() -// } -// fn norm(&self) -> f64 { -// self.norm_l2() -// } -// fn abs(&self) -> Self { -// zipped!(self).map(|unzipped!(xi)| xi.abs()) -// } -// fn as_view(&self) -> Self::View<'_> { -// self.as_ref() -// } -// fn as_view_mut(&mut self) -> Self::ViewMut<'_> { -// self.as_mut() -// } -// fn copy_from(&mut self, other: &Self) { -// self.copy_from(other) -// } -// fn copy_from_view(&mut self, other: &Self::View<'_>) { -// self.copy_from(other) -// } -// fn from_element(nstates: usize, value: Self::T) -> Self { -// faer::Col::from_vec(vec![value; nstates]) -// } -// fn from_vec(vec: Vec) -> Self { -// faer::Col::from_vec(vec) -// } -// fn zeros(nstates: usize) -> Self { -// Self::from_element(nstates, 0.0) -// } -// fn add_scalar_mut(&mut self, scalar: Self::T) { -// self = self + scalar -// } -// fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { -// // faer::linalg::matmul::matmul( -// // self.as_mut(), -// // self.as_ref(), -// // x.as_ref(), -// // Some(beta), -// // alpha, -// // faer::Parallelism::None, -// // ); -// self = self * faer::scale(beta) + x * faer::scale(alpha) -// } -// fn component_mul_assign(&mut self, other: &Self) { -// // faer::linalg::matmul::matmul( -// // self.as_mut(), -// // self.as_ref(), -// // other.as_ref(), -// // None, -// // 1.0, -// // faer::Parallelism::None, -// // ); -// // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); -// self = self * other -// } -// } +impl Vector for faer::Col { + type View<'a> = faer::ColRef<'a, f64>; + type ViewMut<'a> = faer::ColMut<'a, f64>; + type Index = Vec; + fn len(&self) -> IndexType { + self.nrows() + } + fn norm(&self) -> f64 { + self.norm_l2() + } + fn abs(&self) -> Self { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn as_view(&self) -> Self::View<'_> { + self.as_ref() + } + fn as_view_mut(&mut self) -> Self::ViewMut<'_> { + self.as_mut() + } + fn copy_from(&mut self, other: &Self) { + self.copy_from(other) + } + fn copy_from_view(&mut self, other: &Self::View<'_>) { + self.copy_from(other) + } + fn from_element(nstates: usize, value: Self::T) -> Self { + faer::Col::from_vec(vec![value; nstates]) + } + fn from_vec(vec: Vec) -> Self { + faer::Col::from_vec(vec) + } + fn zeros(nstates: usize) -> Self { + Self::from_element(nstates, 0.0) + } + fn add_scalar_mut(&mut self, scalar: Self::T) { + self = self + scalar + } + fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { + // faer::linalg::matmul::matmul( + // self.as_mut(), + // self.as_ref(), + // x.as_ref(), + // Some(beta), + // alpha, + // faer::Parallelism::None, + // ); + self = self * faer::scale(beta) + x * faer::scale(alpha) + } + fn component_mul_assign(&mut self, other: &Self) { + // faer::linalg::matmul::matmul( + // self.as_mut(), + // self.as_ref(), + // other.as_ref(), + // None, + // 1.0, + // faer::Parallelism::None, + // ); + // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); + self = self * other + } +} -// impl VectorIndex for Vec { -// fn zeros(len: IndexType) -> Self { -// vec![0; len as usize] -// } -// fn len(&self) -> IndexType { -// self.len() as IndexType -// } -// } +impl VectorIndex for Vec { + fn zeros(len: IndexType) -> Self { + vec![0; len as usize] + } + fn len(&self) -> IndexType { + self.len() as IndexType + } +} -// impl VectorCommon for faer::Col { -// type T = f64; -// } -// impl<'a> VectorCommon for faer::ColRef<'a, f64> { -// type T = f64; -// } -// impl<'a> VectorCommon for faer::ColMut<'a, f64> { -// type T = f64; -// } +impl VectorCommon for faer::Col { + type T = f64; +} +impl<'a> VectorCommon for faer::ColRef<'a, f64> { + type T = f64; +} +impl<'a> VectorCommon for faer::ColMut<'a, f64> { + type T = f64; +} -// impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { -// type Owned = faer::Col; -// fn abs(&self) -> faer::Col { -// zipped!(self).map(|unzipped!(xi)| xi.abs()) -// } -// fn into_owned(self) -> faer::Col { -// self.to_owned() -// } -// fn scalar_mul(&self, rhs: Self::T) -> Self::Owned { -// self * faer::scale(rhs) -// } -// } +impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { + type Owned = faer::Col; + fn abs(&self) -> faer::Col { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn into_owned(self) -> faer::Col { + self.to_owned() + } + fn scalar_mul(&self, rhs: Self::T) -> Self::Owned { + self * faer::scale(rhs) + } +} -// impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { -// type Owned = faer::Col; -// type View = faer::ColRef<'a, f64>; -// fn abs(&self) -> faer::Col { -// zipped!(self).map(|unzipped!(xi)| xi.abs()) -// } -// fn copy_from(&mut self, other: &Self::Owned) { -// self.copy_from(other); -// } -// fn copy_from_view(&mut self, other: &Self::View) { -// self.copy_from(other); -// } -// } +impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { + type Owned = faer::Col; + type View = faer::ColRef<'a, f64>; + fn abs(&self) -> faer::Col { + zipped!(self).map(|unzipped!(xi)| xi.abs()) + } + fn copy_from(&mut self, other: &Self::Owned) { + self.copy_from(other); + } + fn copy_from_view(&mut self, other: &Self::View) { + self.copy_from(other); + } +} From 93446b0acb1fff88ce471a1ec8bdf055378cb026 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 22:35:55 +0000 Subject: [PATCH 15/28] Working on the Vector trait for faer --- src/scalar/mod.rs | 1 + src/vector/vector_faer.rs | 35 +++++++++++++++++++++++++++-------- 2 files changed, 28 insertions(+), 8 deletions(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index fad4a269..f9ce953b 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -21,6 +21,7 @@ pub trait Scalar: + ClosedDiv + ClosedAdd + Signed + + Into + PartialOrd + Pow + Pow diff --git a/src/vector/vector_faer.rs b/src/vector/vector_faer.rs index e654978b..c21cd09e 100644 --- a/src/vector/vector_faer.rs +++ b/src/vector/vector_faer.rs @@ -8,29 +8,34 @@ use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; impl<'a, T: Scalar> Mul> for faer::ColRef<'a, f64> { type Output = faer::Col; fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { - self * faer::scale(rhs.into()) + self * faer::scale(rhs.value().into()) } } impl<'a, T: Scalar> Mul> for faer::ColMut<'a, f64> { type Output = faer::Col; fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { - self * faer::scale(rhs) + self * faer::scale(rhs.value().into()) } } +impl<'a, T: Scalar> Mul> for faer::Col { + type Output = faer::Col; + fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { + self * faer::scale(rhs.value().into()) + } +} impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { fn mul_assign(&mut self, rhs: crate::scalar::Scale) { - self = self * faer::scale(rhs) + *self = (&*self * faer::scale(rhs.value().into())).as_mut() } } - -impl<'a, T: Scalar> Mul> for faer::Col { - type Output = faer::Col; - fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { - self * faer::scale(rhs) +impl MulAssign> for faer::Col { + fn mul_assign(&mut self, rhs: crate::scalar::Scale) { + *self = &*self * faer::scale(rhs.value().into()) } } +// impl Vector for faer::Col { type View<'a> = faer::ColRef<'a, f64>; @@ -80,6 +85,9 @@ impl Vector for faer::Col { // ); self = self * faer::scale(beta) + x * faer::scale(alpha) } + fn exp(&self) -> Self { + zipped!(self).map(|unzipped!(xi)| xi.exp()) + } fn component_mul_assign(&mut self, other: &Self) { // faer::linalg::matmul::matmul( // self.as_mut(), @@ -92,6 +100,17 @@ impl Vector for faer::Col { // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); self = self * other } + fn component_div_assign(&mut self, other: &Self) { + zipped!(self.as_mut(), other.as_ref()).map(|unzipped!(si, oi)| si /= oi) + } + fn filter_indices bool>(&self, f: F) -> Self::Index { + self.iter() + .enumerate() + .filter_map(|(i, x)| if f(*x) { Some(i as IndexType) } else { None }) + .collect() + } + fn gather_from(&mut self, other: &Self, indices: &Self::Index) {} + fn scatter_from(&mut self, other: &Self, indices: &Self::Index) {} } impl VectorIndex for Vec { From 2f6a47c6f8999c1667d8dcb09adaefab7407ade2 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Tue, 26 Mar 2024 22:46:05 +0000 Subject: [PATCH 16/28] restructure --- src/vector/{vector_faer.rs => faer_serial.rs} | 0 src/vector/mod.rs | 5 +++-- src/vector/{serial.rs => nalgebra_serial.rs} | 0 3 files changed, 3 insertions(+), 2 deletions(-) rename src/vector/{vector_faer.rs => faer_serial.rs} (100%) rename src/vector/{serial.rs => nalgebra_serial.rs} (100%) diff --git a/src/vector/vector_faer.rs b/src/vector/faer_serial.rs similarity index 100% rename from src/vector/vector_faer.rs rename to src/vector/faer_serial.rs diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 893774c1..9c6e4d85 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -4,9 +4,10 @@ use num_traits::Zero; use std::fmt::{Debug, Display}; use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; -mod serial; #[cfg(feature = "faer")] -mod vector_faer; +mod faer_serial; +#[cfg(feature = "nalgebra")] +mod nalgebra_serial; pub trait VectorIndex: Sized + Index + Debug { //+ Display diff --git a/src/vector/serial.rs b/src/vector/nalgebra_serial.rs similarity index 100% rename from src/vector/serial.rs rename to src/vector/nalgebra_serial.rs From 8cd614408f2e8850708ef546eb38d13886fc1203 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 12:51:43 +0000 Subject: [PATCH 17/28] almost done with the Vector definition --- src/vector/faer_serial.rs | 36 ++++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index c21cd09e..2522f2dc 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -27,7 +27,8 @@ impl<'a, T: Scalar> Mul> for faer::Col { } impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { fn mul_assign(&mut self, rhs: crate::scalar::Scale) { - *self = (&*self * faer::scale(rhs.value().into())).as_mut() + let mut aux: faer::Col = self.as_ref() * faer::scale(rhs.value().into()); + *self = aux.as_mut(); } } impl MulAssign> for faer::Col { @@ -72,7 +73,7 @@ impl Vector for faer::Col { Self::from_element(nstates, 0.0) } fn add_scalar_mut(&mut self, scalar: Self::T) { - self = self + scalar + zipped!(self.as_mut()).for_each(|unzipped!(mut s)| *s += scalar) } fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { // faer::linalg::matmul::matmul( @@ -83,7 +84,7 @@ impl Vector for faer::Col { // alpha, // faer::Parallelism::None, // ); - self = self * faer::scale(beta) + x * faer::scale(alpha) + *self = &*self * faer::scale(beta) + x * faer::scale(alpha); } fn exp(&self) -> Self { zipped!(self).map(|unzipped!(xi)| xi.exp()) @@ -97,20 +98,31 @@ impl Vector for faer::Col { // 1.0, // faer::Parallelism::None, // ); - // zipped!(&mut self, &other).for_each(|unzipped!(s, o)| *s = *s + *o); - self = self * other + zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s *= *o); + // self = self * other } fn component_div_assign(&mut self, other: &Self) { - zipped!(self.as_mut(), other.as_ref()).map(|unzipped!(si, oi)| si /= oi) + zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s /= *o); } fn filter_indices bool>(&self, f: F) -> Self::Index { - self.iter() - .enumerate() - .filter_map(|(i, x)| if f(*x) { Some(i as IndexType) } else { None }) - .collect() + let mut indices = vec![]; + for i in 0..self.len() { + if f(self[i]) { + indices.push(i as IndexType); + } + } + indices + } + fn gather_from(&mut self, other: &Self, indices: &Self::Index) { + for (i, &index) in indices.iter().enumerate() { + self[i] = other[index]; + } + } + fn scatter_from(&mut self, other: &Self, indices: &Self::Index) { + for (i, &index) in indices.iter().enumerate() { + self[index] = other[i]; + } } - fn gather_from(&mut self, other: &Self, indices: &Self::Index) {} - fn scatter_from(&mut self, other: &Self, indices: &Self::Index) {} } impl VectorIndex for Vec { From 2c7f4a4863350893f7a6edf59fd7f42914812763 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 13:20:42 +0000 Subject: [PATCH 18/28] initial impl of Vector for faer --- src/scalar/mod.rs | 2 +- src/vector/faer_serial.rs | 31 ++++++++++++++++++++++++++----- src/vector/mod.rs | 6 +++--- 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index f9ce953b..ef155b07 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -6,7 +6,7 @@ use std::{ use nalgebra::{ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, ComplexField, SimdRealField}; use num_traits::{Pow, Signed}; -use crate::vector::{Vector, VectorCommon, VectorRef, VectorView, VectorViewMut}; +use crate::vector::VectorView; pub trait Scalar: nalgebra::Scalar diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index 2522f2dc..d67b6cad 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -1,4 +1,4 @@ -use std::ops::{Mul, MulAssign}; +use std::ops::{Div, Mul, MulAssign}; use faer::{unzipped, zipped}; @@ -27,8 +27,8 @@ impl<'a, T: Scalar> Mul> for faer::Col { } impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { fn mul_assign(&mut self, rhs: crate::scalar::Scale) { - let mut aux: faer::Col = self.as_ref() * faer::scale(rhs.value().into()); - *self = aux.as_mut(); + let scale = faer::scale(rhs.value().into()); + *self *= scale; } } impl MulAssign> for faer::Col { @@ -36,7 +36,15 @@ impl MulAssign> for faer::Col { *self = &*self * faer::scale(rhs.value().into()) } } -// + +impl<'a, T: Scalar> Div> for faer::Col { + type Output = faer::Col; + fn div(self, rhs: crate::scalar::Scale) -> Self::Output { + self + // zipped!(self.as_ref()).map(|unzipped!(xi)| xi / rhs.value()) + // self * faer::scale(rhs.value().into()) + } +} impl Vector for faer::Col { type View<'a> = faer::ColRef<'a, f64>; @@ -67,7 +75,7 @@ impl Vector for faer::Col { faer::Col::from_vec(vec![value; nstates]) } fn from_vec(vec: Vec) -> Self { - faer::Col::from_vec(vec) + faer::Col::from_fn(vec.len(), |i| vec[i]) } fn zeros(nstates: usize) -> Self { Self::from_element(nstates, 0.0) @@ -170,3 +178,16 @@ impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { self.copy_from(other); } } + +// tests +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_abs() { + let v = faer::Col::from_vec(vec![1.0, -2.0, 3.0]); + let v_abs = v.abs(); + assert_eq!(v_abs, faer::Col::from_vec(vec![1.0, 2.0, 3.0])); + } +} diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 9c6e4d85..de073f3e 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -1,10 +1,10 @@ use crate::scalar::Scale; use crate::{IndexType, Scalar}; use num_traits::Zero; -use std::fmt::{Debug, Display}; -use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; +use std::fmt::Debug; +use std::ops::{Add, AddAssign, Div, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; -#[cfg(feature = "faer")] +// #[cfg(feature = "faer")] mod faer_serial; #[cfg(feature = "nalgebra")] mod nalgebra_serial; From b30381ef72deb4f52d0fc245c032ae18bd1dad44 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 13:21:23 +0000 Subject: [PATCH 19/28] add faer feature --- Cargo.toml | 2 +- src/vector/mod.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 0105f522..125035b5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,7 @@ license = "MIT" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [features] -default = ["nalgebra"] +default = ["nalgebra", "faer"] faer = [] nalgebra = [] diffsl = [] diff --git a/src/vector/mod.rs b/src/vector/mod.rs index de073f3e..eda65bca 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -4,7 +4,7 @@ use num_traits::Zero; use std::fmt::Debug; use std::ops::{Add, AddAssign, Div, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}; -// #[cfg(feature = "faer")] +#[cfg(feature = "faer")] mod faer_serial; #[cfg(feature = "nalgebra")] mod nalgebra_serial; From fa4304080ebfe6dab9d4db0b7ee1172d4ebd4f6f Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 13:55:16 +0000 Subject: [PATCH 20/28] cleaning Matrix --- src/matrix/dense_serial.rs | 4 +--- src/matrix/mod.rs | 4 +--- src/matrix/sparse_serial.rs | 1 - src/vector/faer_serial.rs | 8 ++++++++ 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/matrix/dense_serial.rs b/src/matrix/dense_serial.rs index e265761f..a91b530e 100644 --- a/src/matrix/dense_serial.rs +++ b/src/matrix/dense_serial.rs @@ -8,7 +8,6 @@ use super::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut}; impl<'a, T: Scalar> MatrixCommon for DMatrixViewMut<'a, T> { type V = DVector; type T = T; - type O = T; fn ncols(&self) -> IndexType { self.ncols() @@ -32,7 +31,7 @@ impl<'a, T: Scalar> MatrixViewMut<'a> for DMatrixViewMut<'a, T> { impl<'a, T: Scalar> MatrixCommon for DMatrixView<'a, T> { type V = DVector; type T = T; - type O = T; + fn ncols(&self) -> IndexType { self.ncols() } @@ -48,7 +47,6 @@ impl<'a, T: Scalar> MatrixView<'a> for DMatrixView<'a, T> { impl MatrixCommon for DMatrix { type V = DVector; type T = T; - type O = T; fn ncols(&self) -> IndexType { self.ncols() diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs index bfc2f2b8..220cfcac 100644 --- a/src/matrix/mod.rs +++ b/src/matrix/mod.rs @@ -11,7 +11,6 @@ mod sparse_serial; pub trait MatrixCommon: Sized + Debug { type V: Vector; type T: Scalar; - type O; /// Get the number of columns of the matrix fn nrows(&self) -> IndexType; @@ -26,7 +25,6 @@ where { type T = M::T; type V = M::V; - type O = M::O; fn ncols(&self) -> IndexType { M::ncols(self) @@ -42,7 +40,7 @@ where { type T = M::T; type V = M::V; - type O = M::O; + fn ncols(&self) -> IndexType { M::ncols(self) } diff --git a/src/matrix/sparse_serial.rs b/src/matrix/sparse_serial.rs index 2a705c20..ab17c992 100644 --- a/src/matrix/sparse_serial.rs +++ b/src/matrix/sparse_serial.rs @@ -9,7 +9,6 @@ use super::{Matrix, MatrixCommon}; impl MatrixCommon for CscMatrix { type V = DVector; type T = T; - type O = T; fn ncols(&self) -> IndexType { self.ncols() diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index d67b6cad..eccd666d 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -190,4 +190,12 @@ mod tests { let v_abs = v.abs(); assert_eq!(v_abs, faer::Col::from_vec(vec![1.0, 2.0, 3.0])); } + + #[test] + fn test_mult() { + let v = faer::Col::from_vec(vec![1.0, -2.0, 3.0]); + let s = crate::scalar::scale(2.0); + let r = faer::Col::from_vec(vec![2.0, -4.0, 6.0]); + assert_eq!(v * s, r); + } } From c5fb671f8855a76200ff8bacbf3465ad1fe2d59e Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 15:06:41 +0000 Subject: [PATCH 21/28] requested changes --- src/ode_solver/bdf.rs | 8 ++-- src/op/ode.rs | 2 +- src/vector/faer_serial.rs | 96 ++++++++++++++++++--------------------- src/vector/mod.rs | 5 -- 4 files changed, 49 insertions(+), 62 deletions(-) diff --git a/src/ode_solver/bdf.rs b/src/ode_solver/bdf.rs index c41bea24..dc4a24b2 100644 --- a/src/ode_solver/bdf.rs +++ b/src/ode_solver/bdf.rs @@ -289,9 +289,9 @@ where } // update initial step size based on function - let mut scale = state.y.abs(); - scale *= crate::scalar::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 = diff --git a/src/op/ode.rs b/src/op/ode.rs index 788243ec..f1a46615 100644 --- a/src/op/ode.rs +++ b/src/op/ode.rs @@ -97,7 +97,7 @@ impl BdfCallable { for (i, &gamma_i) in gamma.iter().enumerate().take(order + 1).skip(2) { new_psi_neg_y0 += diff.column(i) * scale(gamma_i) } - new_psi_neg_y0 = new_psi_neg_y0 * scale(alpha[order]); + new_psi_neg_y0 *= scale(alpha[order]); // now negate y0 new_psi_neg_y0.sub_assign(y0); diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index eccd666d..dd7e6f65 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -1,38 +1,38 @@ use std::ops::{Div, Mul, MulAssign}; -use faer::{unzipped, zipped}; +use faer::{unzipped, zipped, Col, ColMut, ColRef}; -use crate::{scalar::Scalar, IndexType}; +use crate::{scalar::Scalar, scalar::Scale, IndexType}; use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; -impl<'a, T: Scalar> Mul> for faer::ColRef<'a, f64> { - type Output = faer::Col; - fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +impl<'a, T: Scalar> Mul> for ColRef<'a, f64> { + type Output = Col; + fn mul(self, rhs: Scale) -> Self::Output { self * faer::scale(rhs.value().into()) } } -impl<'a, T: Scalar> Mul> for faer::ColMut<'a, f64> { - type Output = faer::Col; - fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +impl<'a, T: Scalar> Mul> for ColMut<'a, f64> { + type Output = Col; + fn mul(self, rhs: Scale) -> Self::Output { self * faer::scale(rhs.value().into()) } } -impl<'a, T: Scalar> Mul> for faer::Col { - type Output = faer::Col; - fn mul(self, rhs: crate::scalar::Scale) -> Self::Output { +impl<'a, T: Scalar> Mul> for Col { + type Output = Col; + fn mul(self, rhs: Scale) -> Self::Output { self * faer::scale(rhs.value().into()) } } -impl<'a, T: Scalar> MulAssign> for faer::ColMut<'a, f64> { - fn mul_assign(&mut self, rhs: crate::scalar::Scale) { +impl<'a, T: Scalar> MulAssign> for ColMut<'a, f64> { + fn mul_assign(&mut self, rhs: Scale) { let scale = faer::scale(rhs.value().into()); *self *= scale; } } -impl MulAssign> for faer::Col { - fn mul_assign(&mut self, rhs: crate::scalar::Scale) { +impl MulAssign> for Col { + fn mul_assign(&mut self, rhs: Scale) { *self = &*self * faer::scale(rhs.value().into()) } } @@ -46,9 +46,9 @@ impl<'a, T: Scalar> Div> for faer::Col { } } -impl Vector for faer::Col { - type View<'a> = faer::ColRef<'a, f64>; - type ViewMut<'a> = faer::ColMut<'a, f64>; +impl Vector for Col { + type View<'a> = ColRef<'a, f64>; + type ViewMut<'a> = ColMut<'a, f64>; type Index = Vec; fn len(&self) -> IndexType { self.nrows() @@ -72,10 +72,10 @@ impl Vector for faer::Col { self.copy_from(other) } fn from_element(nstates: usize, value: Self::T) -> Self { - faer::Col::from_vec(vec![value; nstates]) + Col::from_vec(vec![value; nstates]) } fn from_vec(vec: Vec) -> Self { - faer::Col::from_fn(vec.len(), |i| vec[i]) + Col::from_fn(vec.len(), |i| vec[i]) } fn zeros(nstates: usize) -> Self { Self::from_element(nstates, 0.0) @@ -84,30 +84,13 @@ impl Vector for faer::Col { zipped!(self.as_mut()).for_each(|unzipped!(mut s)| *s += scalar) } fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) { - // faer::linalg::matmul::matmul( - // self.as_mut(), - // self.as_ref(), - // x.as_ref(), - // Some(beta), - // alpha, - // faer::Parallelism::None, - // ); *self = &*self * faer::scale(beta) + x * faer::scale(alpha); } fn exp(&self) -> Self { zipped!(self).map(|unzipped!(xi)| xi.exp()) } fn component_mul_assign(&mut self, other: &Self) { - // faer::linalg::matmul::matmul( - // self.as_mut(), - // self.as_ref(), - // other.as_ref(), - // None, - // 1.0, - // faer::Parallelism::None, - // ); zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s *= *o); - // self = self * other } fn component_div_assign(&mut self, other: &Self) { zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s /= *o); @@ -142,22 +125,22 @@ impl VectorIndex for Vec { } } -impl VectorCommon for faer::Col { +impl VectorCommon for Col { type T = f64; } -impl<'a> VectorCommon for faer::ColRef<'a, f64> { +impl<'a> VectorCommon for ColRef<'a, f64> { type T = f64; } -impl<'a> VectorCommon for faer::ColMut<'a, f64> { +impl<'a> VectorCommon for ColMut<'a, f64> { type T = f64; } -impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { - type Owned = faer::Col; - fn abs(&self) -> faer::Col { +impl<'a> VectorView<'a> for ColRef<'a, f64> { + type Owned = Col; + fn abs(&self) -> Col { zipped!(self).map(|unzipped!(xi)| xi.abs()) } - fn into_owned(self) -> faer::Col { + fn into_owned(self) -> Col { self.to_owned() } fn scalar_mul(&self, rhs: Self::T) -> Self::Owned { @@ -165,10 +148,10 @@ impl<'a> VectorView<'a> for faer::ColRef<'a, f64> { } } -impl<'a> VectorViewMut<'a> for faer::ColMut<'a, f64> { - type Owned = faer::Col; - type View = faer::ColRef<'a, f64>; - fn abs(&self) -> faer::Col { +impl<'a> VectorViewMut<'a> for ColMut<'a, f64> { + type Owned = Col; + type View = ColRef<'a, f64>; + fn abs(&self) -> Col { zipped!(self).map(|unzipped!(xi)| xi.abs()) } fn copy_from(&mut self, other: &Self::Owned) { @@ -186,16 +169,25 @@ mod tests { #[test] fn test_abs() { - let v = faer::Col::from_vec(vec![1.0, -2.0, 3.0]); + let v = Col::from_vec(vec![1.0, -2.0, 3.0]); let v_abs = v.abs(); - assert_eq!(v_abs, faer::Col::from_vec(vec![1.0, 2.0, 3.0])); + assert_eq!(v_abs, Col::from_vec(vec![1.0, 2.0, 3.0])); } #[test] fn test_mult() { - let v = faer::Col::from_vec(vec![1.0, -2.0, 3.0]); + let v = Col::from_vec(vec![1.0, -2.0, 3.0]); let s = crate::scalar::scale(2.0); - let r = faer::Col::from_vec(vec![2.0, -4.0, 6.0]); + let r = Col::from_vec(vec![2.0, -4.0, 6.0]); assert_eq!(v * s, r); } + + #[test] + fn test_mul_assign() { + let mut v = Col::from_vec(vec![1.0, -2.0, 3.0]); + let s = crate::scalar::scale(2.0); + let r = Col::from_vec(vec![2.0, -4.0, 6.0]); + v.mul_assign(s); + assert_eq!(v, r); + } } diff --git a/src/vector/mod.rs b/src/vector/mod.rs index eda65bca..d0c838d4 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -10,7 +10,6 @@ mod faer_serial; mod nalgebra_serial; pub trait VectorIndex: Sized + Index + Debug { - //+ Display fn zeros(len: IndexType) -> Self; fn len(&self) -> IndexType; fn is_empty(&self) -> bool { @@ -19,7 +18,6 @@ pub trait VectorIndex: Sized + Index + Debug { } pub trait VectorCommon: Sized + Debug { - // + Display type T: Scalar; } @@ -75,7 +73,6 @@ pub trait VectorViewMut<'a>: + for<'b> VectorMutOpsByValue<&'b Self::View> + for<'b> VectorMutOpsByValue<&'b Self::Owned> + MulAssign> - // + DivAssign + Index + IndexMut { @@ -92,7 +89,6 @@ pub trait VectorView<'a>: + for<'b> VectorOpsByValue<&'b Self::Owned, Self::Owned> + for<'b> VectorOpsByValue<&'b Self, Self::Owned> + Mul, Output = Self::Owned> - // + Div + Index { type Owned; @@ -113,7 +109,6 @@ pub trait Vector: + for<'b> VectorMutOpsByValue<&'b Self> + for<'a, 'b> VectorMutOpsByValue<&'b Self::View<'a>> + MulAssign> - // + DivAssign + Index + IndexMut + Clone From 844fb31373e292de20092ebff4d6d5ece902c31e Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 15:13:47 +0000 Subject: [PATCH 22/28] requested changes for mul_assign --- src/vector/nalgebra_serial.rs | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/vector/nalgebra_serial.rs b/src/vector/nalgebra_serial.rs index de3370a2..0e59a632 100644 --- a/src/vector/nalgebra_serial.rs +++ b/src/vector/nalgebra_serial.rs @@ -60,16 +60,12 @@ impl Mul> for DVector { impl MulAssign> for DVector { fn mul_assign(&mut self, rhs: Scale) { - *self = &*self * rhs.value(); + *self *= rhs.value(); } } impl<'a, T: Scalar> MulAssign> for DVectorViewMut<'a, T> { fn mul_assign(&mut self, rhs: Scale) { - //I'm doing this horrendous hack because - //somehow I'm missing which trait to implement so I can multiply - //a &mut DVector by a scalar - let a = &*self * rhs.value(); - self.copy_from(&a); + *self *= rhs.value(); } } From 6155771ae348190d79cd7cd33d3416b1da064657 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 15:23:58 +0000 Subject: [PATCH 23/28] new macros to reduce duplicated code on the trait implemented by Scalar --- src/scalar/mod.rs | 69 +++++++++++++++++++---------------------------- 1 file changed, 27 insertions(+), 42 deletions(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index ef155b07..96606c88 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -58,54 +58,39 @@ pub fn scale(value: E) -> Scale { Scale(value) } -//TODO: Is it possible for us to need RhsE != LhsE? -impl Mul> for Scale { - type Output = Scale; - - #[inline] - fn mul(self, rhs: Scale) -> Self::Output { - Scale(self.0 * rhs.0) - } +macro_rules! impl_bin_op { + ($trait:ident, $method:ident, $operator:tt) => { + impl $trait> for Scale { + type Output = Scale; + + #[inline] + fn $method(self, rhs: Scale) -> Self::Output { + Scale(self.0 $operator rhs.0) + } + } + }; } -impl Add> for Scale { - type Output = Scale; - - #[inline] - fn add(self, rhs: Scale) -> Self::Output { - Scale(self.0 + rhs.0) - } -} - -impl Sub> for Scale { - type Output = Scale; - - #[inline] - fn sub(self, rhs: Scale) -> Self::Output { - Scale(self.0 - rhs.0) - } +macro_rules! impl_assign_bin_op { + ($trait:ident, $method:ident, $operator:tt) => { + impl $trait> for Scale { + #[inline] + fn $method(&mut self, rhs: Scale) { + self.0 = self.0 $operator rhs.0 + } + } + }; } -impl MulAssign> for Scale { - #[inline] - fn mul_assign(&mut self, rhs: Scale) { - self.0 = self.0 * rhs.0 - } -} +//TODO: Is it possible for us to need RhsE != LhsE? -impl AddAssign> for Scale { - #[inline] - fn add_assign(&mut self, rhs: Scale) { - self.0 = self.0 + rhs.0 - } -} +impl_bin_op!(Mul, mul, *); +impl_bin_op!(Add, add, +); +impl_bin_op!(Sub, sub, -); -impl SubAssign> for Scale { - #[inline] - fn sub_assign(&mut self, rhs: Scale) { - self.0 = self.0 - rhs.0 - } -} +impl_assign_bin_op!(MulAssign, mul_assign, *); +impl_assign_bin_op!(AddAssign, add_assign, +); +impl_assign_bin_op!(SubAssign, sub_assign, -); impl PartialEq for Scale { #[inline] From c8e4e86ff663e4bfcc054884d28e180f6b9fcfe6 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 18:07:25 +0000 Subject: [PATCH 24/28] impl Div> for Col --- src/vector/faer_serial.rs | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index dd7e6f65..ea47d8c0 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -40,9 +40,7 @@ impl MulAssign> for Col { impl<'a, T: Scalar> Div> for faer::Col { type Output = faer::Col; fn div(self, rhs: crate::scalar::Scale) -> Self::Output { - self - // zipped!(self.as_ref()).map(|unzipped!(xi)| xi / rhs.value()) - // self * faer::scale(rhs.value().into()) + zipped!(self).map(|unzipped!(xi)| *xi / rhs.value().into()) } } From 96a0db678c7185bf30a746bdf890f1bd69d76de2 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Wed, 27 Mar 2024 18:08:09 +0000 Subject: [PATCH 25/28] clean --- src/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index eaba57fd..5b1dde78 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -43,7 +43,7 @@ use nonlinear_solver::{newton::NewtonNonlinearSolver, NonLinearSolver}; pub use ode_solver::{ bdf::Bdf, equations::OdeEquations, OdeSolverMethod, OdeSolverProblem, OdeSolverState, }; -use op::{LinearOp, NonLinearOp}; +use op::NonLinearOp; use scalar::{IndexType, Scalar}; use solver::SolverProblem; use vector::{Vector, VectorIndex, VectorRef, VectorView, VectorViewMut}; From 4401d891dcb37b0c0722a12d3820cdb3d236e699 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Thu, 28 Mar 2024 11:41:10 +0000 Subject: [PATCH 26/28] Changes requested --- src/scalar/mod.rs | 15 ------------ src/vector/faer_serial.rs | 45 ++++++++++++++++++----------------- src/vector/nalgebra_serial.rs | 25 +++++++++---------- 3 files changed, 36 insertions(+), 49 deletions(-) diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index 96606c88..7820d064 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -82,8 +82,6 @@ macro_rules! impl_assign_bin_op { }; } -//TODO: Is it possible for us to need RhsE != LhsE? - impl_bin_op!(Mul, mul, *); impl_bin_op!(Add, add, +); impl_bin_op!(Sub, sub, -); @@ -107,19 +105,6 @@ impl<'a, V: VectorView<'static, T = E>, E: Scalar> Mul for Scale { } } -//TODO: Not sure why it is now working -// impl<'a, E, V> Mul> for V -// where -// V: VectorView<'static, T = E>, -// E: Scalar, -// { -// type Output = V::Owned; -// #[inline] -// fn mul(self, rhs: Scale) -> Self::Output { -// self.scalar_mul(rhs.0) -// } -// } - #[test] fn test_scale() { assert_eq!(scale(2.0) * scale(3.0), scale(6.0)); diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index ea47d8c0..22c689bb 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -5,30 +5,37 @@ use faer::{unzipped, zipped, Col, ColMut, ColRef}; use crate::{scalar::Scalar, scalar::Scale, IndexType}; use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; -impl<'a, T: Scalar> Mul> for ColRef<'a, f64> { - type Output = Col; - fn mul(self, rhs: Scale) -> Self::Output { - self * faer::scale(rhs.value().into()) - } + +macro_rules! impl_op_for_struct { + ($struct:ident, $trait_name:ident, $func_name:ident) => { + impl<'a, T: Scalar> $trait_name> for $struct<'a, f64> { + type Output = Col; + + fn $func_name(self, rhs: Scale) -> Self::Output { + self * faer::scale(rhs.value().into()) + } + } + }; } -impl<'a, T: Scalar> Mul> for ColMut<'a, f64> { +impl_op_for_struct!(ColRef, Mul, mul); +impl_op_for_struct!(ColMut, Mul, mul); + +impl Mul> for Col { type Output = Col; fn mul(self, rhs: Scale) -> Self::Output { self * faer::scale(rhs.value().into()) } } - -impl<'a, T: Scalar> Mul> for Col { - type Output = Col; - fn mul(self, rhs: Scale) -> Self::Output { - self * faer::scale(rhs.value().into()) +impl<'a, T: Scalar> Div> for faer::Col { + type Output = faer::Col; + fn div(self, rhs: Scale) -> Self::Output { + zipped!(self).map(|unzipped!(xi)| *xi / rhs.value().into()) } } impl<'a, T: Scalar> MulAssign> for ColMut<'a, f64> { fn mul_assign(&mut self, rhs: Scale) { - let scale = faer::scale(rhs.value().into()); - *self *= scale; + *self *= faer::scale(rhs.value().into()); } } impl MulAssign> for Col { @@ -37,13 +44,6 @@ impl MulAssign> for Col { } } -impl<'a, T: Scalar> Div> for faer::Col { - type Output = faer::Col; - fn div(self, rhs: crate::scalar::Scale) -> Self::Output { - zipped!(self).map(|unzipped!(xi)| *xi / rhs.value().into()) - } -} - impl Vector for Col { type View<'a> = ColRef<'a, f64>; type ViewMut<'a> = ColMut<'a, f64>; @@ -164,6 +164,7 @@ impl<'a> VectorViewMut<'a> for ColMut<'a, f64> { #[cfg(test)] mod tests { use super::*; + use crate::scalar::scale; #[test] fn test_abs() { @@ -175,7 +176,7 @@ mod tests { #[test] fn test_mult() { let v = Col::from_vec(vec![1.0, -2.0, 3.0]); - let s = crate::scalar::scale(2.0); + let s = scale(2.0); let r = Col::from_vec(vec![2.0, -4.0, 6.0]); assert_eq!(v * s, r); } @@ -183,7 +184,7 @@ mod tests { #[test] fn test_mul_assign() { let mut v = Col::from_vec(vec![1.0, -2.0, 3.0]); - let s = crate::scalar::scale(2.0); + let s = scale(2.0); let r = Col::from_vec(vec![2.0, -4.0, 6.0]); v.mul_assign(s); assert_eq!(v, r); diff --git a/src/vector/nalgebra_serial.rs b/src/vector/nalgebra_serial.rs index 0e59a632..c0fcb00f 100644 --- a/src/vector/nalgebra_serial.rs +++ b/src/vector/nalgebra_serial.rs @@ -39,25 +39,26 @@ impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { self * rhs } } -impl<'a, T: Scalar> Mul> for DVectorView<'a, T> { - type Output = DVector; - fn mul(self, rhs: Scale) -> Self::Output { - self * rhs.value() - } -} -impl<'a, T: Scalar> Mul> for DVectorViewMut<'a, T> { - type Output = DVector; - fn mul(self, rhs: Scale) -> Self::Output { - self * rhs.value() - } +macro_rules! impl_op_for_dvector_struct { + ($struct:ident, $trait_name:ident, $func_name:ident) => { + impl<'a, T: Scalar> $trait_name> for $struct<'a, T> { + type Output = DVector; + fn $func_name(self, rhs: Scale) -> Self::Output { + self * rhs.value() + } + } + }; } + +impl_op_for_dvector_struct!(DVectorView, Mul, mul); +impl_op_for_dvector_struct!(DVectorViewMut, Mul, mul); + impl Mul> for DVector { type Output = DVector; fn mul(self, rhs: Scale) -> Self::Output { self * rhs.value() } } - impl MulAssign> for DVector { fn mul_assign(&mut self, rhs: Scale) { *self *= rhs.value(); From a13a1de99be814f69a002d22cb6c2b5622e12daa Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Thu, 28 Mar 2024 11:51:39 +0000 Subject: [PATCH 27/28] Consistent macro naming --- src/vector/faer_serial.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs index 22c689bb..73a1d29f 100644 --- a/src/vector/faer_serial.rs +++ b/src/vector/faer_serial.rs @@ -6,7 +6,7 @@ use crate::{scalar::Scalar, scalar::Scale, IndexType}; use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; -macro_rules! impl_op_for_struct { +macro_rules! impl_op_for_faer_struct { ($struct:ident, $trait_name:ident, $func_name:ident) => { impl<'a, T: Scalar> $trait_name> for $struct<'a, f64> { type Output = Col; @@ -18,8 +18,8 @@ macro_rules! impl_op_for_struct { }; } -impl_op_for_struct!(ColRef, Mul, mul); -impl_op_for_struct!(ColMut, Mul, mul); +impl_op_for_faer_struct!(ColRef, Mul, mul); +impl_op_for_faer_struct!(ColMut, Mul, mul); impl Mul> for Col { type Output = Col; From aec7b7fff246a9b2e3d4bc98b87fb494a25e85c7 Mon Sep 17 00:00:00 2001 From: "Julian D. Otalvaro" Date: Thu, 28 Mar 2024 11:53:06 +0000 Subject: [PATCH 28/28] Some order --- src/vector/nalgebra_serial.rs | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/vector/nalgebra_serial.rs b/src/vector/nalgebra_serial.rs index c0fcb00f..c2d53e7f 100644 --- a/src/vector/nalgebra_serial.rs +++ b/src/vector/nalgebra_serial.rs @@ -6,6 +6,20 @@ use crate::{scalar::Scale, IndexType, Scalar}; use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut}; +macro_rules! impl_op_for_dvector_struct { + ($struct:ident, $trait_name:ident, $func_name:ident) => { + impl<'a, T: Scalar> $trait_name> for $struct<'a, T> { + type Output = DVector; + fn $func_name(self, rhs: Scale) -> Self::Output { + self * rhs.value() + } + } + }; +} + +impl_op_for_dvector_struct!(DVectorView, Mul, mul); +impl_op_for_dvector_struct!(DVectorViewMut, Mul, mul); + impl VectorIndex for DVector { fn zeros(len: IndexType) -> Self { DVector::from_element(len, 0) @@ -39,19 +53,6 @@ impl<'a, T: Scalar> VectorView<'a> for DVectorView<'a, T> { self * rhs } } -macro_rules! impl_op_for_dvector_struct { - ($struct:ident, $trait_name:ident, $func_name:ident) => { - impl<'a, T: Scalar> $trait_name> for $struct<'a, T> { - type Output = DVector; - fn $func_name(self, rhs: Scale) -> Self::Output { - self * rhs.value() - } - } - }; -} - -impl_op_for_dvector_struct!(DVectorView, Mul, mul); -impl_op_for_dvector_struct!(DVectorViewMut, Mul, mul); impl Mul> for DVector { type Output = DVector;