From 144dfbd555c17b5f6c866472a9320e975705ec41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 2 Feb 2018 12:26:27 +0100 Subject: [PATCH] Add quadform/cmpy/cdpy. --- CHANGELOG.md | 8 +++ src/core/blas.rs | 119 ++++++++++++++++++++++++++++++++------ src/core/componentwise.rs | 44 ++++++++++++-- tests/core/blas.rs | 17 +++++- 4 files changed, 164 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 664ceeaa..37a276ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,10 @@ documented here. This project adheres to [Semantic Versioning](http://semver.org/). ## [0.14.0] − WIP +### Modified + * `quadform` has been renamed `quadform_tr`. The new `quadform` method takes + the matrix on the right-hand-side instead of the matrix on the + left-hand-side of the quadratic form. ### Added * The `mint` feature that can be enabled in order to allow conversions from and to types of the [mint](https://crates.io/crates/mint) crate. @@ -12,6 +16,10 @@ This project adheres to [Semantic Versioning](http://semver.org/). `::from_element(...)`. * The `.iamin()` methods that returns the index of the vector entry with smallest absolute value. + * Add blas-like operations: `cmpy, cdpy` for componentwise multiplicatons and + division with scalar factors: + - `self <- alpha * self + beta * a * b` + - `self <- alpha * self + beta / a * b` * `UnitQuaternion::scaled_rotation_between_axis` and `UnitQuaternion::rotation_between_axis` that take Unit vectors instead of Vector as arguments. diff --git a/src/core/blas.rs b/src/core/blas.rs index 65dbe7e8..82948c1b 100644 --- a/src/core/blas.rs +++ b/src/core/blas.rs @@ -338,10 +338,10 @@ impl Vector /// If `beta` is zero, `self` is never read. #[inline] pub fn gemv_tr(&mut self, - alpha: N, - a: &Matrix, - x: &Vector, - beta: N) + alpha: N, + a: &Matrix, + x: &Vector, + beta: N) where N: One, SB: Storage, SC: Storage, @@ -481,6 +481,35 @@ impl> Matrix } } } + + /// Computes `self = alpha * a.transpose() * b + beta * self`, where `a, b, self` are matrices. + /// `alpha` and `beta` are scalar. + /// + /// If `beta` is zero, `self` is never read. + #[inline] + pub fn gemm_tr(&mut self, + alpha: N, + a: &Matrix, + b: &Matrix, + beta: N) + where N: One, + SB: Storage, + SC: Storage, + ShapeConstraint: SameNumberOfRows + + SameNumberOfColumns + + AreMultipliable { + let (nrows1, ncols1) = self.shape(); + let (nrows2, ncols2) = a.shape(); + let (nrows3, ncols3) = b.shape(); + + assert_eq!(nrows2, nrows3, "gemm: dimensions mismatch for multiplication."); + assert_eq!((nrows1, ncols1), (ncols2, ncols3), "gemm: dimensions mismatch for addition."); + + for j1 in 0 .. ncols1 { + // FIXME: avoid bound checks. + self.column_mut(j1).gemv_tr(alpha, a, &b.column(j1), beta); + } + } } @@ -520,13 +549,15 @@ impl> Matrix impl> SquareMatrix where N: Scalar + Zero + One + ClosedAdd + ClosedMul { - /// Computes the quadratic form `self = alpha * lrs * mid * lhs.transpose() + beta * self`. - pub fn quadform_with_workspace(&mut self, - work: &mut Vector, - alpha: N, - lhs: &Matrix, - mid: &SquareMatrix, - beta: N) + /// Computes the quadratic form `self = alpha * lhs * mid * lhs.transpose() + beta * self`. + /// + /// This uses the provided workspace `work` to avoid allocations for intermediate results. + pub fn quadform_tr_with_workspace(&mut self, + work: &mut Vector, + alpha: N, + lhs: &Matrix, + mid: &SquareMatrix, + beta: N) where D2: Dim, R3: Dim, C3: Dim, D4: Dim, S2: StorageMut, S3: Storage, @@ -544,18 +575,68 @@ impl> SquareMatrix } } - /// Computes the quadratic form `self = alpha * lrs * mid * lhs.transpose() + beta * self`. - pub fn quadform(&mut self, - alpha: N, - lhs: &Matrix, - mid: &SquareMatrix, - beta: N) - where R3: Dim, C3: Dim, D4: Dim, + /// Computes the quadratic form `self = alpha * lhs * mid * lhs.transpose() + beta * self`. + /// + /// This allocates a workspace vector of dimension D1 for intermediate results. + /// Use `.quadform_tr_with_workspace(...)` instead to avoid allocations. + pub fn quadform_tr(&mut self, + alpha: N, + lhs: &Matrix, + mid: &SquareMatrix, + beta: N) + where R3: Dim, C3: Dim, D4: Dim, S3: Storage, S4: Storage, ShapeConstraint: DimEq + DimEq + DimEq, DefaultAllocator: Allocator { let mut work = unsafe { Vector::new_uninitialized_generic(self.data.shape().0, U1) }; - self.quadform_with_workspace(&mut work, alpha, lhs, mid, beta) + self.quadform_tr_with_workspace(&mut work, alpha, lhs, mid, beta) + } + + /// Computes the quadratic form `self = alpha * rhs.transpose() * mid * rhs + beta * self`. + /// + /// This uses the provided workspace `work` to avoid allocations for intermediate results. + pub fn quadform_with_workspace(&mut self, + work: &mut Vector, + alpha: N, + mid: &SquareMatrix, + rhs: &Matrix, + beta: N) + where D2: Dim, D3: Dim, R4: Dim, C4: Dim, + S2: StorageMut, + S3: Storage, + S4: Storage, + ShapeConstraint: DimEq + + DimEq + + DimEq + + AreMultipliable { + work.gemv(N::one(), mid, &rhs.column(0), N::zero()); + self.column_mut(0).gemv_tr(alpha, &rhs, work, beta); + + for j in 1 .. rhs.ncols() { + work.gemv(N::one(), mid, &rhs.column(j), N::zero()); + self.column_mut(j).gemv_tr(alpha, &rhs, work, beta); + } + } + + /// Computes the quadratic form `self = alpha * rhs.transpose() * mid * rhs + beta * self`. + /// + /// This allocates a workspace vector of dimension D2 for intermediate results. + /// Use `.quadform_with_workspace(...)` instead to avoid allocations. + pub fn quadform(&mut self, + alpha: N, + mid: &SquareMatrix, + rhs: &Matrix, + beta: N) + where D2: Dim, R3: Dim, C3: Dim, + S2: Storage, + S3: Storage, + ShapeConstraint: DimEq + + DimEq + + AreMultipliable, + DefaultAllocator: Allocator { + + let mut work = unsafe { Vector::new_uninitialized_generic(mid.data.shape().0, U1) }; + self.quadform_with_workspace(&mut work, alpha, mid, rhs, beta) } } diff --git a/src/core/componentwise.rs b/src/core/componentwise.rs index 65a4d1a9..56113148 100644 --- a/src/core/componentwise.rs +++ b/src/core/componentwise.rs @@ -1,6 +1,7 @@ // Non-convensional componentwise operators. -use num::Signed; +use std::ops::{Add, Mul}; +use num::{Zero, Signed}; use alga::general::{ClosedMul, ClosedDiv}; @@ -33,7 +34,7 @@ impl> Matrix { } macro_rules! component_binop_impl( - ($($binop: ident, $binop_mut: ident, $binop_assign: ident, $Trait: ident . $op_assign: ident, $desc:expr, $desc_mut:expr);* $(;)*) => {$( + ($($binop: ident, $binop_mut: ident, $binop_assign: ident, $cbpy: ident, $Trait: ident . $op: ident . $op_assign: ident, $desc:expr, $desc_mut:expr);* $(;)*) => {$( impl> Matrix { #[doc = $desc] #[inline] @@ -60,6 +61,41 @@ macro_rules! component_binop_impl( } impl> Matrix { + // componentwise binop plus Y. + #[inline] + pub fn $cbpy(&mut self, alpha: N, a: &Matrix, b: &Matrix, beta: N) + where N: $Trait + Zero + Mul + Add, + R2: Dim, C2: Dim, + R3: Dim, C3: Dim, + SB: Storage, + SC: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + + SameNumberOfRows + SameNumberOfColumns { + assert_eq!(self.shape(), a.shape(), "Componentwise mul/div: mismatched matrix dimensions."); + assert_eq!(self.shape(), b.shape(), "Componentwise mul/div: mismatched matrix dimensions."); + + if beta.is_zero() { + for j in 0 .. self.ncols() { + for i in 0 .. self.nrows() { + unsafe { + let res = alpha * a.get_unchecked(i, j).$op(*b.get_unchecked(i, j)); + *self.get_unchecked_mut(i, j) = res; + } + } + } + } + else { + for j in 0 .. self.ncols() { + for i in 0 .. self.nrows() { + unsafe { + let res = alpha * a.get_unchecked(i, j).$op(*b.get_unchecked(i, j)); + *self.get_unchecked_mut(i, j) = beta * *self.get_unchecked(i, j) + res; + } + } + } + } + } + #[doc = $desc_mut] #[inline] pub fn $binop_assign(&mut self, rhs: &Matrix) @@ -96,9 +132,9 @@ macro_rules! component_binop_impl( ); component_binop_impl!( - component_mul, component_mul_mut, component_mul_assign, ClosedMul.mul_assign, + component_mul, component_mul_mut, component_mul_assign, cmpy, ClosedMul.mul.mul_assign, "Componentwise matrix multiplication.", "Mutable, componentwise matrix multiplication."; - component_div, component_div_mut, component_div_assign, ClosedDiv.div_assign, + component_div, component_div_mut, component_div_assign, cdpy, ClosedDiv.div.div_assign, "Componentwise matrix division.", "Mutable, componentwise matrix division."; // FIXME: add other operators like bitshift, etc. ? ); diff --git a/tests/core/blas.rs b/tests/core/blas.rs index 0db8bc88..e37cb049 100644 --- a/tests/core/blas.rs +++ b/tests/core/blas.rs @@ -74,6 +74,21 @@ quickcheck! { } fn quadform(n: usize, alpha: f64, beta: f64) -> bool { + let n = cmp::max(1, cmp::min(n, 50)); + let rhs = DMatrix::::new_random(6, n); + let mid = DMatrix::::new_random(6, 6); + let mut res = DMatrix::new_random(n, n); + + let expected = &res * beta + rhs.transpose() * &mid * &rhs * alpha; + + res.quadform(alpha, &mid, &rhs, beta); + + println!("{}{}", res, expected); + + relative_eq!(res, expected, epsilon = 1.0e-7) + } + + fn quadform_tr(n: usize, alpha: f64, beta: f64) -> bool { let n = cmp::max(1, cmp::min(n, 50)); let lhs = DMatrix::::new_random(6, n); let mid = DMatrix::::new_random(n, n); @@ -81,7 +96,7 @@ quickcheck! { let expected = &res * beta + &lhs * &mid * lhs.transpose() * alpha; - res.quadform(alpha, &lhs, &mid , beta); + res.quadform_tr(alpha, &lhs, &mid , beta); println!("{}{}", res, expected);