diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs index ad76e245..dc21c12b 100644 --- a/nalgebra-lapack/src/lu.rs +++ b/nalgebra-lapack/src/lu.rs @@ -214,7 +214,7 @@ where } } - /// Solves the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to + /// Solves the linear system `self.adjoint() * x = b`, where `x` is the unknown to /// be determined. pub fn solve_conjugate_transpose( &self, @@ -249,7 +249,7 @@ where self.generic_solve_mut(b'T', b) } - /// Solves in-place the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to + /// Solves in-place the linear system `self.adjoint() * x = b`, where `x` is the unknown to /// be determined. /// /// Returns `false` if no solution was found (the decomposed matrix is singular). diff --git a/src/base/blas.rs b/src/base/blas.rs index 7dce36b6..210db8cd 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -11,7 +11,7 @@ use base::constraint::{ }; use base::dimension::{Dim, Dynamic, U1, U2, U3, U4}; use base::storage::{Storage, StorageMut}; -use base::{DefaultAllocator, Matrix, Scalar, SquareMatrix, Vector}; +use base::{DefaultAllocator, Matrix, Scalar, SquareMatrix, Vector, DVectorSlice}; // FIXME: find a way to avoid code duplication just for complex number support. @@ -32,11 +32,11 @@ impl> Vector { pub fn icamax(&self) -> usize { assert!(!self.is_empty(), "The input vector must not be empty."); - let mut the_max = unsafe { self.vget_unchecked(0).asum() }; + let mut the_max = unsafe { self.vget_unchecked(0).norm1() }; let mut the_i = 0; for i in 1..self.nrows() { - let val = unsafe { self.vget_unchecked(i).asum() }; + let val = unsafe { self.vget_unchecked(i).norm1() }; if val > the_max { the_max = val; @@ -211,12 +211,12 @@ impl> Matrix { pub fn icamax_full(&self) -> (usize, usize) { assert!(!self.is_empty(), "The input matrix must not be empty."); - let mut the_max = unsafe { self.get_unchecked((0, 0)).asum() }; + let mut the_max = unsafe { self.get_unchecked((0, 0)).norm1() }; let mut the_ij = (0, 0); for j in 0..self.ncols() { for i in 0..self.nrows() { - let val = unsafe { self.get_unchecked((i, j)).asum() }; + let val = unsafe { self.get_unchecked((i, j)).norm1() }; if val > the_max { the_max = val; @@ -263,13 +263,11 @@ impl> Matri } } -impl> Matrix { - /// The dot product between two complex or real vectors or matrices (seen as vectors). - /// - /// This is the same as `.dot` except that the conjugate of each component of `self` is taken - /// before performing the products. - #[inline] - pub fn cdot(&self, rhs: &Matrix) -> N +impl> Matrix +where N: Scalar + Zero + ClosedAdd + ClosedMul +{ + #[inline(always)] + fn dotx(&self, rhs: &Matrix, conjugate: impl Fn(N) -> N) -> N where SB: Storage, ShapeConstraint: DimEq + DimEq, @@ -283,27 +281,27 @@ impl> Matrix { // because the `for` loop below won't be very efficient on those. if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { unsafe { - let a = self.get_unchecked((0, 0)).conjugate() * *rhs.get_unchecked((0, 0)); - let b = self.get_unchecked((1, 0)).conjugate() * *rhs.get_unchecked((1, 0)); + let a = conjugate(*self.get_unchecked((0, 0))) * *rhs.get_unchecked((0, 0)); + let b = conjugate(*self.get_unchecked((1, 0))) * *rhs.get_unchecked((1, 0)); return a + b; } } if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { unsafe { - let a = self.get_unchecked((0, 0)).conjugate() * *rhs.get_unchecked((0, 0)); - let b = self.get_unchecked((1, 0)).conjugate() * *rhs.get_unchecked((1, 0)); - let c = self.get_unchecked((2, 0)).conjugate() * *rhs.get_unchecked((2, 0)); + let a = conjugate(*self.get_unchecked((0, 0))) * *rhs.get_unchecked((0, 0)); + let b = conjugate(*self.get_unchecked((1, 0))) * *rhs.get_unchecked((1, 0)); + let c = conjugate(*self.get_unchecked((2, 0))) * *rhs.get_unchecked((2, 0)); return a + b + c; } } if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { unsafe { - let mut a = self.get_unchecked((0, 0)).conjugate() * *rhs.get_unchecked((0, 0)); - let mut b = self.get_unchecked((1, 0)).conjugate() * *rhs.get_unchecked((1, 0)); - let c = self.get_unchecked((2, 0)).conjugate() * *rhs.get_unchecked((2, 0)); - let d = self.get_unchecked((3, 0)).conjugate() * *rhs.get_unchecked((3, 0)); + let mut a = conjugate(*self.get_unchecked((0, 0))) * *rhs.get_unchecked((0, 0)); + let mut b = conjugate(*self.get_unchecked((1, 0))) * *rhs.get_unchecked((1, 0)); + let c = conjugate(*self.get_unchecked((2, 0))) * *rhs.get_unchecked((2, 0)); + let d = conjugate(*self.get_unchecked((3, 0))) * *rhs.get_unchecked((3, 0)); a += c; b += d; @@ -343,14 +341,14 @@ impl> Matrix { acc7 = N::zero(); while self.nrows() - i >= 8 { - acc0 += unsafe { self.get_unchecked((i + 0, j)).conjugate() * *rhs.get_unchecked((i + 0, j)) }; - acc1 += unsafe { self.get_unchecked((i + 1, j)).conjugate() * *rhs.get_unchecked((i + 1, j)) }; - acc2 += unsafe { self.get_unchecked((i + 2, j)).conjugate() * *rhs.get_unchecked((i + 2, j)) }; - acc3 += unsafe { self.get_unchecked((i + 3, j)).conjugate() * *rhs.get_unchecked((i + 3, j)) }; - acc4 += unsafe { self.get_unchecked((i + 4, j)).conjugate() * *rhs.get_unchecked((i + 4, j)) }; - acc5 += unsafe { self.get_unchecked((i + 5, j)).conjugate() * *rhs.get_unchecked((i + 5, j)) }; - acc6 += unsafe { self.get_unchecked((i + 6, j)).conjugate() * *rhs.get_unchecked((i + 6, j)) }; - acc7 += unsafe { self.get_unchecked((i + 7, j)).conjugate() * *rhs.get_unchecked((i + 7, j)) }; + acc0 += unsafe { conjugate(*self.get_unchecked((i + 0, j))) * *rhs.get_unchecked((i + 0, j)) }; + acc1 += unsafe { conjugate(*self.get_unchecked((i + 1, j))) * *rhs.get_unchecked((i + 1, j)) }; + acc2 += unsafe { conjugate(*self.get_unchecked((i + 2, j))) * *rhs.get_unchecked((i + 2, j)) }; + acc3 += unsafe { conjugate(*self.get_unchecked((i + 3, j))) * *rhs.get_unchecked((i + 3, j)) }; + acc4 += unsafe { conjugate(*self.get_unchecked((i + 4, j))) * *rhs.get_unchecked((i + 4, j)) }; + acc5 += unsafe { conjugate(*self.get_unchecked((i + 5, j))) * *rhs.get_unchecked((i + 5, j)) }; + acc6 += unsafe { conjugate(*self.get_unchecked((i + 6, j))) * *rhs.get_unchecked((i + 6, j)) }; + acc7 += unsafe { conjugate(*self.get_unchecked((i + 7, j))) * *rhs.get_unchecked((i + 7, j)) }; i += 8; } @@ -360,17 +358,14 @@ impl> Matrix { res += acc3 + acc7; for k in i..self.nrows() { - res += unsafe { self.get_unchecked((k, j)).conjugate() * *rhs.get_unchecked((k, j)) } + res += unsafe { conjugate(*self.get_unchecked((k, j))) * *rhs.get_unchecked((k, j)) } } } res } -} -impl> Matrix -where N: Scalar + Zero + ClosedAdd + ClosedMul -{ + /// The dot product between two vectors or matrices (seen as vectors). /// /// Note that this is **not** the matrix multiplication as in, e.g., numpy. For matrix @@ -396,97 +391,36 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul SB: Storage, ShapeConstraint: DimEq + DimEq, { - assert!( - self.nrows() == rhs.nrows(), - "Dot product dimensions mismatch." - ); + self.dotx(rhs, |e| e) + } - // So we do some special cases for common fixed-size vectors of dimension lower than 8 - // because the `for` loop below won't be very efficient on those. - if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { - unsafe { - let a = *self.get_unchecked((0, 0)) * *rhs.get_unchecked((0, 0)); - let b = *self.get_unchecked((1, 0)) * *rhs.get_unchecked((1, 0)); - - return a + b; - } - } - if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { - unsafe { - let a = *self.get_unchecked((0, 0)) * *rhs.get_unchecked((0, 0)); - let b = *self.get_unchecked((1, 0)) * *rhs.get_unchecked((1, 0)); - let c = *self.get_unchecked((2, 0)) * *rhs.get_unchecked((2, 0)); - - return a + b + c; - } - } - if (R::is::() || R2::is::()) && (C::is::() || C2::is::()) { - unsafe { - let mut a = *self.get_unchecked((0, 0)) * *rhs.get_unchecked((0, 0)); - let mut b = *self.get_unchecked((1, 0)) * *rhs.get_unchecked((1, 0)); - let c = *self.get_unchecked((2, 0)) * *rhs.get_unchecked((2, 0)); - let d = *self.get_unchecked((3, 0)) * *rhs.get_unchecked((3, 0)); - - a += c; - b += d; - - return a + b; - } - } - - // All this is inspired from the "unrolled version" discussed in: - // http://blog.theincredibleholk.org/blog/2012/12/10/optimizing-dot-product/ - // - // And this comment from bluss: - // https://users.rust-lang.org/t/how-to-zip-two-slices-efficiently/2048/12 - let mut res = N::zero(); - - // We have to define them outside of the loop (and not inside at first assignment) - // otherwise vectorization won't kick in for some reason. - let mut acc0; - let mut acc1; - let mut acc2; - let mut acc3; - let mut acc4; - let mut acc5; - let mut acc6; - let mut acc7; - - for j in 0..self.ncols() { - let mut i = 0; - - acc0 = N::zero(); - acc1 = N::zero(); - acc2 = N::zero(); - acc3 = N::zero(); - acc4 = N::zero(); - acc5 = N::zero(); - acc6 = N::zero(); - acc7 = N::zero(); - - while self.nrows() - i >= 8 { - acc0 += unsafe { *self.get_unchecked((i + 0, j)) * *rhs.get_unchecked((i + 0, j)) }; - acc1 += unsafe { *self.get_unchecked((i + 1, j)) * *rhs.get_unchecked((i + 1, j)) }; - acc2 += unsafe { *self.get_unchecked((i + 2, j)) * *rhs.get_unchecked((i + 2, j)) }; - acc3 += unsafe { *self.get_unchecked((i + 3, j)) * *rhs.get_unchecked((i + 3, j)) }; - acc4 += unsafe { *self.get_unchecked((i + 4, j)) * *rhs.get_unchecked((i + 4, j)) }; - acc5 += unsafe { *self.get_unchecked((i + 5, j)) * *rhs.get_unchecked((i + 5, j)) }; - acc6 += unsafe { *self.get_unchecked((i + 6, j)) * *rhs.get_unchecked((i + 6, j)) }; - acc7 += unsafe { *self.get_unchecked((i + 7, j)) * *rhs.get_unchecked((i + 7, j)) }; - i += 8; - } - - res += acc0 + acc4; - res += acc1 + acc5; - res += acc2 + acc6; - res += acc3 + acc7; - - for k in i..self.nrows() { - res += unsafe { *self.get_unchecked((k, j)) * *rhs.get_unchecked((k, j)) } - } - } - - res + /// The dot product between two vectors or matrices (seen as vectors). + /// + /// Note that this is **not** the matrix multiplication as in, e.g., numpy. For matrix + /// multiplication, use one of: `.gemm`, `.mul_to`, `.mul`, the `*` operator. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Vector3, Matrix2x3}; + /// let vec1 = Vector3::new(1.0, 2.0, 3.0); + /// let vec2 = Vector3::new(0.1, 0.2, 0.3); + /// assert_eq!(vec1.dot(&vec2), 1.4); + /// + /// let mat1 = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// let mat2 = Matrix2x3::new(0.1, 0.2, 0.3, + /// 0.4, 0.5, 0.6); + /// assert_eq!(mat1.dot(&mat2), 9.1); + /// ``` + #[inline] + pub fn dotc(&self, rhs: &Matrix) -> N + where + N: Complex, + SB: Storage, + ShapeConstraint: DimEq + DimEq, + { + self.dotx(rhs, Complex::conjugate) } /// The dot product between the transpose of `self` and `rhs`. @@ -643,40 +577,15 @@ where } } - /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a - /// vector, and `alpha, beta` two scalars. - /// - /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part - /// (including the diagonal) is actually read. - /// - /// # Examples: - /// - /// ``` - /// # use nalgebra::{Matrix2, Vector2}; - /// let mat = Matrix2::new(1.0, 2.0, - /// 2.0, 4.0); - /// let mut vec1 = Vector2::new(1.0, 2.0); - /// let vec2 = Vector2::new(0.1, 0.2); - /// vec1.gemv_symm(10.0, &mat, &vec2, 5.0); - /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); - /// - /// - /// // The matrix upper-triangular elements can be garbage because it is never - /// // read by this method. Therefore, it is not necessary for the caller to - /// // fill the matrix struct upper-triangle. - /// let mat = Matrix2::new(1.0, 9999999.9999999, - /// 2.0, 4.0); - /// let mut vec1 = Vector2::new(1.0, 2.0); - /// vec1.gemv_symm(10.0, &mat, &vec2, 5.0); - /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); - /// ``` - #[inline] - pub fn gemv_symm( + + #[inline(always)] + fn xgemv( &mut self, alpha: N, a: &SquareMatrix, x: &Vector, beta: N, + dotc: impl Fn(&DVectorSlice, &DVectorSlice) -> N, ) where N: One, SB: Storage, @@ -687,83 +596,6 @@ where let dim2 = a.nrows(); let dim3 = x.nrows(); - assert!( - a.is_square(), - "Symmetric gemv: the input matrix must be square." - ); - assert!( - dim2 == dim3 && dim1 == dim2, - "Symmetric gemv: dimensions mismatch." - ); - - if dim2 == 0 { - return; - } - - // FIXME: avoid bound checks. - let col2 = a.column(0); - let val = unsafe { *x.vget_unchecked(0) }; - self.axpy(alpha * val, &col2, beta); - self[0] += alpha * x.rows_range(1..).dot(&a.slice_range(1.., 0)); - - for j in 1..dim2 { - let col2 = a.column(j); - let dot = x.rows_range(j..).dot(&col2.rows_range(j..)); - - let val; - unsafe { - val = *x.vget_unchecked(j); - *self.vget_unchecked_mut(j) += alpha * dot; - } - self.rows_range_mut(j + 1..) - .axpy(alpha * val, &col2.rows_range(j + 1..), N::one()); - } - } - - /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a - /// vector, and `alpha, beta` two scalars. - /// - /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part - /// (including the diagonal) is actually read. - /// - /// # Examples: - /// - /// ``` - /// # use nalgebra::{Matrix2, Vector2}; - /// let mat = Matrix2::new(1.0, 2.0, - /// 2.0, 4.0); - /// let mut vec1 = Vector2::new(1.0, 2.0); - /// let vec2 = Vector2::new(0.1, 0.2); - /// vec1.gemv_symm(10.0, &mat, &vec2, 5.0); - /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); - /// - /// - /// // The matrix upper-triangular elements can be garbage because it is never - /// // read by this method. Therefore, it is not necessary for the caller to - /// // fill the matrix struct upper-triangle. - /// let mat = Matrix2::new(1.0, 9999999.9999999, - /// 2.0, 4.0); - /// let mut vec1 = Vector2::new(1.0, 2.0); - /// vec1.gemv_symm(10.0, &mat, &vec2, 5.0); - /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); - /// ``` - #[inline] - pub fn cgemv_symm( - &mut self, - alpha: N, - a: &SquareMatrix, - x: &Vector, - beta: N, - ) where - N: Complex, - SB: Storage, - SC: Storage, - ShapeConstraint: DimEq + AreMultipliable, - { - let dim1 = self.nrows(); - let dim2 = a.nrows(); - let dim3 = x.nrows(); - assert!( a.is_square(), "Symmetric cgemv: the input matrix must be square." @@ -781,11 +613,11 @@ where let col2 = a.column(0); let val = unsafe { *x.vget_unchecked(0) }; self.axpy(alpha * val, &col2, beta); - self[0] += alpha * a.slice_range(1.., 0).cdot(&x.rows_range(1..)); + self[0] += alpha * dotc(&a.slice_range(1.., 0), &x.rows_range(1..)); for j in 1..dim2 { let col2 = a.column(j); - let dot = col2.rows_range(j..).cdot(&x.rows_range(j..)); + let dot = dotc(&col2.rows_range(j..), &x.rows_range(j..)); let val; unsafe { @@ -797,6 +629,111 @@ where } } + /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a + /// vector, and `alpha, beta` two scalars. DEPRECATED: use `sygemv` instead. + #[inline] + #[deprecated(note = "This is renamed `sygemv` to match the original BLAS terminology.")] + pub fn gemv_symm( + &mut self, + alpha: N, + a: &SquareMatrix, + x: &Vector, + beta: N, + ) where + N: One, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + self.sygemv(alpha, a, x, beta) + } + + /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a + /// vector, and `alpha, beta` two scalars. + /// + /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part + /// (including the diagonal) is actually read. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Matrix2, Vector2}; + /// let mat = Matrix2::new(1.0, 2.0, + /// 2.0, 4.0); + /// let mut vec1 = Vector2::new(1.0, 2.0); + /// let vec2 = Vector2::new(0.1, 0.2); + /// vec1.sygemv(10.0, &mat, &vec2, 5.0); + /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); + /// + /// + /// // The matrix upper-triangular elements can be garbage because it is never + /// // read by this method. Therefore, it is not necessary for the caller to + /// // fill the matrix struct upper-triangle. + /// let mat = Matrix2::new(1.0, 9999999.9999999, + /// 2.0, 4.0); + /// let mut vec1 = Vector2::new(1.0, 2.0); + /// vec1.sygemv(10.0, &mat, &vec2, 5.0); + /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); + /// ``` + #[inline] + pub fn sygemv( + &mut self, + alpha: N, + a: &SquareMatrix, + x: &Vector, + beta: N, + ) where + N: One, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + self.xgemv(alpha, a, x, beta, |a, b| a.dot(b)) + } + + /// Computes `self = alpha * a * x + beta * self`, where `a` is an **hermitian** matrix, `x` a + /// vector, and `alpha, beta` two scalars. + /// + /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part + /// (including the diagonal) is actually read. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Matrix2, Vector2}; + /// let mat = Matrix2::new(1.0, 2.0, + /// 2.0, 4.0); + /// let mut vec1 = Vector2::new(1.0, 2.0); + /// let vec2 = Vector2::new(0.1, 0.2); + /// vec1.sygemv(10.0, &mat, &vec2, 5.0); + /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); + /// + /// + /// // The matrix upper-triangular elements can be garbage because it is never + /// // read by this method. Therefore, it is not necessary for the caller to + /// // fill the matrix struct upper-triangle. + /// let mat = Matrix2::new(1.0, 9999999.9999999, + /// 2.0, 4.0); + /// let mut vec1 = Vector2::new(1.0, 2.0); + /// vec1.sygemv(10.0, &mat, &vec2, 5.0); + /// assert_eq!(vec1, Vector2::new(10.0, 20.0)); + /// ``` + #[inline] + pub fn hegemv( + &mut self, + alpha: N, + a: &SquareMatrix, + x: &Vector, + beta: N, + ) where + N: Complex, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + self.xgemv(alpha, a, x, beta, |a, b| a.dotc(b)) + } + /// Computes `self = alpha * a.transpose() * x + beta * self`, where `a` is a matrix, `x` a vector, and /// `alpha, beta` two scalars. /// @@ -855,33 +792,17 @@ where } } -// FIXME: duplicate code impl> Matrix - where N: Complex + Zero + ClosedAdd + ClosedMul +where N: Scalar + Zero + ClosedAdd + ClosedMul { - /// Computes `self = alpha * x * y.transpose() + beta * self`. - /// - /// If `beta` is zero, `self` is never read. - /// - /// # Examples: - /// - /// ``` - /// # use nalgebra::{Matrix2x3, Vector2, Vector3}; - /// let mut mat = Matrix2x3::repeat(4.0); - /// let vec1 = Vector2::new(1.0, 2.0); - /// let vec2 = Vector3::new(0.1, 0.2, 0.3); - /// let expected = vec1 * vec2.transpose() * 10.0 + mat * 5.0; - /// - /// mat.ger(10.0, &vec1, &vec2, 5.0); - /// assert_eq!(mat, expected); - /// ``` - #[inline] - pub fn gerc( + #[inline(always)] + fn gerx( &mut self, alpha: N, x: &Vector, y: &Vector, beta: N, + conjugate: impl Fn(N) -> N, ) where N: One, SB: Storage, @@ -899,15 +820,11 @@ impl> Matrix for j in 0..ncols1 { // FIXME: avoid bound checks. - let val = unsafe { y.vget_unchecked(j).conjugate() }; + let val = unsafe { conjugate(*y.vget_unchecked(j)) }; self.column_mut(j).axpy(alpha * val, x, beta); } } -} -impl> Matrix -where N: Scalar + Zero + ClosedAdd + ClosedMul -{ /// Computes `self = alpha * x * y.transpose() + beta * self`. /// /// If `beta` is zero, `self` is never read. @@ -937,20 +854,39 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul SC: Storage, ShapeConstraint: DimEq + DimEq, { - let (nrows1, ncols1) = self.shape(); - let dim2 = x.nrows(); - let dim3 = y.nrows(); + self.gerx(alpha, x, y, beta, |e| e) + } - assert!( - nrows1 == dim2 && ncols1 == dim3, - "ger: dimensions mismatch." - ); - - for j in 0..ncols1 { - // FIXME: avoid bound checks. - let val = unsafe { *y.vget_unchecked(j) }; - self.column_mut(j).axpy(alpha * val, x, beta); - } + /// Computes `self = alpha * x * y.transpose() + beta * self`. + /// + /// If `beta` is zero, `self` is never read. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Matrix2x3, Vector2, Vector3}; + /// let mut mat = Matrix2x3::repeat(4.0); + /// let vec1 = Vector2::new(1.0, 2.0); + /// let vec2 = Vector3::new(0.1, 0.2, 0.3); + /// let expected = vec1 * vec2.transpose() * 10.0 + mat * 5.0; + /// + /// mat.ger(10.0, &vec1, &vec2, 5.0); + /// assert_eq!(mat, expected); + /// ``` + #[inline] + pub fn gerc( + &mut self, + alpha: N, + x: &Vector, + y: &Vector, + beta: N, + ) where + N: Complex, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + DimEq, + { + self.gerx(alpha, x, y, beta, Complex::conjugate) } /// Computes `self = alpha * a * b + beta * self`, where `a, b, self` are matrices. @@ -1146,6 +1082,42 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul impl> Matrix where N: Scalar + Zero + ClosedAdd + ClosedMul { + #[inline(always)] + fn sygerx( + &mut self, + alpha: N, + x: &Vector, + y: &Vector, + beta: N, + conjugate: impl Fn(N) -> N, + ) where + N: One, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + DimEq, + { + let dim1 = self.nrows(); + let dim2 = x.nrows(); + let dim3 = y.nrows(); + + assert!( + self.is_square(), + "Symmetric ger: the input matrix must be square." + ); + assert!(dim1 == dim2 && dim1 == dim3, "ger: dimensions mismatch."); + + for j in 0..dim1 { + let val = unsafe { conjugate(*y.vget_unchecked(j)) }; + let subdim = Dynamic::new(dim1 - j); + // FIXME: avoid bound checks. + self.generic_slice_mut((j, j), (subdim, U1)).axpy( + alpha * val, + &x.rows_range(j..), + beta, + ); + } + } + /// Computes `self = alpha * x * y.transpose() + beta * self`, where `self` is a **symmetric** /// matrix. /// @@ -1166,6 +1138,7 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul /// assert_eq!(mat.lower_triangle(), expected.lower_triangle()); /// assert_eq!(mat.m12, 99999.99999); // This was untouched. #[inline] + #[deprecated(note = "This is renamed `syger` to match the original BLAS terminology.")] pub fn ger_symm( &mut self, alpha: N, @@ -1178,26 +1151,77 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul SC: Storage, ShapeConstraint: DimEq + DimEq, { - let dim1 = self.nrows(); - let dim2 = x.nrows(); - let dim3 = y.nrows(); + self.syger(alpha, x, y, beta) + } - assert!( - self.is_square(), - "Symmetric ger: the input matrix must be square." - ); - assert!(dim1 == dim2 && dim1 == dim3, "ger: dimensions mismatch."); + /// Computes `self = alpha * x * y.transpose() + beta * self`, where `self` is a **symmetric** + /// matrix. + /// + /// If `beta` is zero, `self` is never read. The result is symmetric. Only the lower-triangular + /// (including the diagonal) part of `self` is read/written. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Matrix2, Vector2}; + /// let mut mat = Matrix2::identity(); + /// let vec1 = Vector2::new(1.0, 2.0); + /// let vec2 = Vector2::new(0.1, 0.2); + /// let expected = vec1 * vec2.transpose() * 10.0 + mat * 5.0; + /// mat.m12 = 99999.99999; // This component is on the upper-triangular part and will not be read/written. + /// + /// mat.ger_symm(10.0, &vec1, &vec2, 5.0); + /// assert_eq!(mat.lower_triangle(), expected.lower_triangle()); + /// assert_eq!(mat.m12, 99999.99999); // This was untouched. + #[inline] + pub fn syger( + &mut self, + alpha: N, + x: &Vector, + y: &Vector, + beta: N, + ) where + N: One, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + DimEq, + { + self.sygerx(alpha, x, y, beta, |e| e) + } - for j in 0..dim1 { - let val = unsafe { *y.vget_unchecked(j) }; - let subdim = Dynamic::new(dim1 - j); - // FIXME: avoid bound checks. - self.generic_slice_mut((j, j), (subdim, U1)).axpy( - alpha * val, - &x.rows_range(j..), - beta, - ); - } + /// Computes `self = alpha * x * y.transpose() + beta * self`, where `self` is a **symmetric** + /// matrix. + /// + /// If `beta` is zero, `self` is never read. The result is symmetric. Only the lower-triangular + /// (including the diagonal) part of `self` is read/written. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::{Matrix2, Vector2}; + /// let mut mat = Matrix2::identity(); + /// let vec1 = Vector2::new(1.0, 2.0); + /// let vec2 = Vector2::new(0.1, 0.2); + /// let expected = vec1 * vec2.transpose() * 10.0 + mat * 5.0; + /// mat.m12 = 99999.99999; // This component is on the upper-triangular part and will not be read/written. + /// + /// mat.ger_symm(10.0, &vec1, &vec2, 5.0); + /// assert_eq!(mat.lower_triangle(), expected.lower_triangle()); + /// assert_eq!(mat.m12, 99999.99999); // This was untouched. + #[inline] + pub fn hegerc( + &mut self, + alpha: N, + x: &Vector, + y: &Vector, + beta: N, + ) where + N: Complex, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + DimEq, + { + self.sygerx(alpha, x, y, beta, Complex::conjugate) } } diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 8f6c071a..6a89e945 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -914,9 +914,9 @@ impl> Matrix { } impl> Matrix { - /// Takes the conjugate and transposes `self` and store the result into `out`. + /// Takes the adjoint (aka. conjugate-transpose) of `self` and store the result into `out`. #[inline] - pub fn conjugate_transpose_to(&self, out: &mut Matrix) + pub fn adjoint_to(&self, out: &mut Matrix) where R2: Dim, C2: Dim, @@ -939,20 +939,41 @@ impl> Matrix { } } - /// The conjugate transposition of `self`. + /// The adjoint (aka. conjugate-transpose) of `self`. #[inline] - pub fn conjugate_transpose(&self) -> MatrixMN + pub fn adjoint(&self) -> MatrixMN where DefaultAllocator: Allocator { let (nrows, ncols) = self.data.shape(); unsafe { let mut res: MatrixMN<_, C, R> = Matrix::new_uninitialized_generic(ncols, nrows); - self.conjugate_transpose_to(&mut res); + self.adjoint_to(&mut res); res } } + /// Takes the conjugate and transposes `self` and store the result into `out`. + #[deprecated(note = "Renamed `self.adjoint_to(out)`.")] + #[inline] + pub fn conjugate_transpose_to(&self, out: &mut Matrix) + where + R2: Dim, + C2: Dim, + SB: StorageMut, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + self.adjoint_to(out) + } + + /// The conjugate transposition of `self`. + #[deprecated(note = "Renamed `self.adjoint()`.")] + #[inline] + pub fn conjugate_transpose(&self) -> MatrixMN + where DefaultAllocator: Allocator { + self.adjoint() + } + /// The conjugate of `self`. #[inline] pub fn conjugate(&self) -> MatrixMN @@ -1088,13 +1109,13 @@ impl> SquareMatrix { tr } - /// The hermitian part of `self`, i.e., `0.5 * (self + self.conjugate_transpose())`. + /// The hermitian part of `self`, i.e., `0.5 * (self + self.adjoint())`. #[inline] pub fn hermitian_part(&self) -> MatrixMN where DefaultAllocator: Allocator { assert!(self.is_square(), "Cannot compute the hermitian part of a non-square matrix."); - let mut tr = self.conjugate_transpose(); + let mut tr = self.adjoint(); tr += self; tr *= ::convert::<_, N>(0.5); tr @@ -1522,7 +1543,7 @@ impl> Matrix { SB: Storage, ShapeConstraint: DimEq + DimEq, { - let prod = self.cdot(other); + let prod = self.dotc(other); let n1 = self.norm(); let n2 = other.norm(); @@ -1593,7 +1614,7 @@ impl> Unit> { where DefaultAllocator: Allocator, { - let c_hang = self.cdot(rhs).real(); + let c_hang = self.dotc(rhs).real(); // self == other if c_hang.abs() >= N::Real::one() { diff --git a/src/base/matrix_alga.rs b/src/base/matrix_alga.rs index a2e1dd54..4d540f46 100644 --- a/src/base/matrix_alga.rs +++ b/src/base/matrix_alga.rs @@ -192,7 +192,7 @@ where DefaultAllocator: Allocator #[inline] fn inner_product(&self, other: &Self) -> N { - self.cdot(other) + self.dotc(other) } } diff --git a/src/base/norm.rs b/src/base/norm.rs index b3056067..2a50f4ec 100644 --- a/src/base/norm.rs +++ b/src/base/norm.rs @@ -33,7 +33,7 @@ impl Norm for EuclideanNorm { #[inline] fn norm(&self, m: &Matrix) -> N::Real where R: Dim, C: Dim, S: Storage { - m.cdot(m).real().sqrt() + m.dotc(m).real().sqrt() } #[inline] @@ -101,7 +101,7 @@ impl> Matrix { for i in 0..self.ncols() { let col = self.column(i); - res += col.cdot(&col).real() + res += col.dotc(&col).real() } res diff --git a/src/base/ops.rs b/src/base/ops.rs index bee584b6..f2a683c1 100644 --- a/src/base/ops.rs +++ b/src/base/ops.rs @@ -13,7 +13,7 @@ use base::constraint::{ }; use base::dimension::{Dim, DimMul, DimName, DimProd}; use base::storage::{ContiguousStorageMut, Storage, StorageMut}; -use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, MatrixSum, Scalar}; +use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, MatrixSum, Scalar, VectorSliceN}; /* * @@ -629,13 +629,28 @@ where res } - /// Equivalent to `self.transpose() * rhs` but stores the result into `out` to avoid - /// allocations. + /// Equivalent to `self.adjoint() * rhs`. #[inline] - pub fn tr_mul_to( + pub fn ad_mul(&self, rhs: &Matrix) -> MatrixMN + where + N: Complex, + SB: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = + unsafe { Matrix::new_uninitialized_generic(self.data.shape().1, rhs.data.shape().1) }; + + self.ad_mul_to(rhs, &mut res); + res + } + + #[inline(always)] + fn xx_mul_to( &self, rhs: &Matrix, out: &mut Matrix, + dot: impl Fn(&VectorSliceN, &VectorSliceN) -> N, ) where SB: Storage, SC: StorageMut, @@ -656,12 +671,43 @@ where for i in 0..ncols1 { for j in 0..ncols2 { - let dot = self.column(i).dot(&rhs.column(j)); + let dot = dot(&self.column(i), &rhs.column(j)); unsafe { *out.get_unchecked_mut((i, j)) = dot }; } } } + /// Equivalent to `self.transpose() * rhs` but stores the result into `out` to avoid + /// allocations. + #[inline] + pub fn tr_mul_to( + &self, + rhs: &Matrix, + out: &mut Matrix, + ) where + SB: Storage, + SC: StorageMut, + ShapeConstraint: SameNumberOfRows + DimEq + DimEq, + { + self.xx_mul_to(rhs, out, |a, b| a.dot(b)) + } + + /// Equivalent to `self.adjoint() * rhs` but stores the result into `out` to avoid + /// allocations. + #[inline] + pub fn ad_mul_to( + &self, + rhs: &Matrix, + out: &mut Matrix, + ) where + N: Complex, + SB: Storage, + SC: StorageMut, + ShapeConstraint: SameNumberOfRows + DimEq + DimEq, + { + self.xx_mul_to(rhs, out, |a, b| a.dotc(b)) + } + /// Equivalent to `self * rhs` but stores the result into `out` to avoid allocations. #[inline] pub fn mul_to( @@ -760,35 +806,16 @@ where } } -// XXX: avoid code duplication. -impl> Matrix { - /// Returns the absolute value of the component with the largest absolute value. - #[inline] - pub fn camax(&self) -> N::Real { - let mut max = N::Real::zero(); +impl> Matrix { + #[inline(always)] + fn xcmp(&self, abs: impl Fn(N) -> N2, cmp: impl Fn(N2, N2) -> bool) -> N2 + where N2: Scalar + PartialOrd + Zero { + let mut max = N2::zero(); for e in self.iter() { - let ae = e.asum(); + let ae = abs(*e); - if ae > max { - max = ae; - } - } - - max - } -} - -impl> Matrix { - /// Returns the absolute value of the component with the largest absolute value. - #[inline] - pub fn amax(&self) -> N { - let mut max = N::zero(); - - for e in self.iter() { - let ae = e.abs(); - - if ae > max { + if cmp(ae, max) { max = ae; } } @@ -796,61 +823,45 @@ impl> Matri max } - /// Returns the absolute value of the component with the smallest absolute value. + /// Returns the absolute value of the component with the largest absolute value. #[inline] - pub fn amin(&self) -> N { - let mut it = self.iter(); - let mut min = it - .next() - .expect("amin: empty matrices not supported.") - .abs(); + pub fn amax(&self) -> N + where N: PartialOrd + Signed { + self.xcmp(|e| e.abs(), |a, b| a > b) + } - for e in it { - let ae = e.abs(); - - if ae < min { - min = ae; - } - } - - min + /// Returns the the 1-norm of the complex component with the largest 1-norm. + #[inline] + pub fn camax(&self) -> N::Real + where N: Complex { + self.xcmp(|e| e.norm1(), |a, b| a > b) } /// Returns the component with the largest value. #[inline] - pub fn max(&self) -> N { - let mut it = self.iter(); - let mut max = it - .next() - .expect("max: empty matrices not supported."); + pub fn max(&self) -> N + where N: PartialOrd + Signed { + self.xcmp(|e| e, |a, b| a > b) + } - for e in it { - let ae = e; + /// Returns the absolute value of the component with the smallest absolute value. + #[inline] + pub fn amin(&self) -> N + where N: PartialOrd + Signed { + self.xcmp(|e| e.abs(), |a, b| a < b) + } - if ae > max { - max = ae; - } - } - - *max + /// Returns the the 1-norm of the complex component with the smallest 1-norm. + #[inline] + pub fn camin(&self) -> N::Real + where N: Complex { + self.xcmp(|e| e.norm1(), |a, b| a < b) } /// Returns the component with the smallest value. #[inline] - pub fn min(&self) -> N { - let mut it = self.iter(); - let mut min = it - .next() - .expect("min: empty matrices not supported."); - - for e in it { - let ae = e; - - if ae < min { - min = ae; - } - } - - *min + pub fn min(&self) -> N + where N: PartialOrd + Signed { + self.xcmp(|e| e, |a, b| a < b) } -} +} \ No newline at end of file diff --git a/src/base/properties.rs b/src/base/properties.rs index a511c0f8..de6be72a 100644 --- a/src/base/properties.rs +++ b/src/base/properties.rs @@ -97,8 +97,7 @@ impl> Matrix { N::Epsilon: Copy, DefaultAllocator: Allocator + Allocator, { - // FIXME: add a conjugate-transpose-mul - (self.conjugate().tr_mul(self)).is_identity(eps) + (self.ad_mul(self)).is_identity(eps) } } diff --git a/src/debug/random_sdp.rs b/src/debug/random_sdp.rs index 772eccc9..f8ee6d30 100644 --- a/src/debug/random_sdp.rs +++ b/src/debug/random_sdp.rs @@ -31,7 +31,7 @@ where DefaultAllocator: Allocator /// random reals generators. pub fn new N>(dim: D, mut rand: Rand) -> Self { let mut m = RandomOrthogonal::new(dim, || rand()).unwrap(); - let mt = m.conjugate_transpose(); + let mt = m.adjoint(); for i in 0..dim.value() { let mut col = m.column_mut(i); diff --git a/src/geometry/reflection.rs b/src/geometry/reflection.rs index b3c77bd9..653a9bdc 100644 --- a/src/geometry/reflection.rs +++ b/src/geometry/reflection.rs @@ -35,7 +35,7 @@ impl> Reflection { D: DimName, DefaultAllocator: Allocator, { - let bias = axis.cdot(&pt.coords); + let bias = axis.dotc(&pt.coords); Self::new(axis, bias) } @@ -56,7 +56,7 @@ impl> Reflection { // dot product, and then mutably. Somehow, this allows significantly // better optimizations of the dot product from the compiler. let m_two: N = ::convert(-2.0f64); - let factor = (self.axis.cdot(&rhs.column(i)) - self.bias) * m_two; + let factor = (self.axis.dotc(&rhs.column(i)) - self.bias) * m_two; rhs.column_mut(i).axpy(factor, &self.axis, N::one()); } } @@ -73,7 +73,7 @@ impl> Reflection { // dot product, and then mutably. Somehow, this allows significantly // better optimizations of the dot product from the compiler. let m_two = sign.scale(::convert(-2.0f64)); - let factor = (self.axis.cdot(&rhs.column(i)) - self.bias) * m_two; + let factor = (self.axis.dotc(&rhs.column(i)) - self.bias) * m_two; rhs.column_mut(i).axpy(factor, &self.axis, sign); } } diff --git a/src/lib.rs b/src/lib.rs index a718cf91..c2bb04f2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -170,6 +170,7 @@ pub use alga::general::{Id, Real, Complex}; /// Gets the ubiquitous multiplicative identity element. /// /// Same as `Id::new()`. +#[deprecated(note = "use `Id::new()` instead.")] #[inline] pub fn id() -> Id { Id::new() @@ -416,6 +417,7 @@ pub fn partial_sort2<'a, T: PartialOrd>(a: &'a T, b: &'a T) -> Option<(&'a T, &' /// # See also: /// /// * [`inverse`](fn.inverse.html) +#[deprecated(note = "use the `.try_inverse()` method instead")] #[inline] pub fn try_inverse(m: &M) -> Option { m.try_inverse() @@ -426,6 +428,7 @@ pub fn try_inverse(m: &M) -> Option { /// # See also: /// /// * [`try_inverse`](fn.try_inverse.html) +#[deprecated(note = "use the `.inverse()` method instead")] #[inline] pub fn inverse>(m: &M) -> M { m.two_sided_inverse() diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 682bd990..d4973df1 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -121,7 +121,7 @@ where DefaultAllocator: Allocator ShapeConstraint: SameNumberOfRows, { let _ = self.chol.solve_lower_triangular_mut(b); - let _ = self.chol.conjugate().tr_solve_lower_triangular_mut(b); + let _ = self.chol.ad_solve_lower_triangular_mut(b); } /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 6c7f7194..4418b283 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -30,6 +30,6 @@ pub use self::lu::*; pub use self::permutation_sequence::*; pub use self::qr::*; pub use self::schur::*; -//pub use self::svd::*; +pub use self::svd::*; pub use self::symmetric_eigen::*; pub use self::symmetric_tridiagonal::*; diff --git a/src/linalg/solve.rs b/src/linalg/solve.rs index b43b44df..fad85a81 100644 --- a/src/linalg/solve.rs +++ b/src/linalg/solve.rs @@ -4,7 +4,7 @@ use base::allocator::Allocator; use base::constraint::{SameNumberOfRows, ShapeConstraint}; use base::dimension::{Dim, U1}; use base::storage::{Storage, StorageMut}; -use base::{DefaultAllocator, Matrix, MatrixMN, SquareMatrix, Vector}; +use base::{DefaultAllocator, Matrix, MatrixMN, SquareMatrix, Vector, DVectorSlice}; impl> SquareMatrix { /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only @@ -180,10 +180,9 @@ impl> SquareMatrix { /* * - * Transpose versions + * Transpose and adjoint versions * */ - /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only /// the lower-triangular part of `self` (including the diagonal) is considered not-zero. #[inline] @@ -191,10 +190,10 @@ impl> SquareMatrix { &self, b: &Matrix, ) -> Option> - where - S2: StorageMut, - DefaultAllocator: Allocator, - ShapeConstraint: SameNumberOfRows, + where + S2: StorageMut, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, { let mut res = b.clone_owned(); if self.tr_solve_lower_triangular_mut(&mut res) { @@ -211,10 +210,10 @@ impl> SquareMatrix { &self, b: &Matrix, ) -> Option> - where - S2: StorageMut, - DefaultAllocator: Allocator, - ShapeConstraint: SameNumberOfRows, + where + S2: StorageMut, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, { let mut res = b.clone_owned(); if self.tr_solve_upper_triangular_mut(&mut res) { @@ -230,14 +229,14 @@ impl> SquareMatrix { &self, b: &mut Matrix, ) -> bool - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, { let cols = b.ncols(); for i in 0..cols { - if !self.tr_solve_lower_triangular_vector_mut(&mut b.column_mut(i)) { + if !self.xx_solve_lower_triangular_vector_mut(&mut b.column_mut(i), |e| e, |a, b| a.dot(b)) { return false; } } @@ -245,46 +244,20 @@ impl> SquareMatrix { true } - fn tr_solve_lower_triangular_vector_mut(&self, b: &mut Vector) -> bool - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, - { - let dim = self.nrows(); - - for i in (0..dim).rev() { - let dot = self.slice_range(i + 1.., i).dot(&b.slice_range(i + 1.., 0)); - - unsafe { - let b_i = b.vget_unchecked_mut(i); - - let diag = *self.get_unchecked((i, i)); - - if diag.is_zero() { - return false; - } - - *b_i = (*b_i - dot) / diag; - } - } - - true - } - /// Solves the linear system `self.transpose() . x = b` where `x` is the unknown and only the /// upper-triangular part of `self` (including the diagonal) is considered not-zero. pub fn tr_solve_upper_triangular_mut( &self, b: &mut Matrix, ) -> bool - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, { let cols = b.ncols(); for i in 0..cols { - if !self.tr_solve_upper_triangular_vector_mut(&mut b.column_mut(i)) { + if !self.xx_solve_upper_triangular_vector_mut(&mut b.column_mut(i), |e| e, |a, b| a.dot(b)) { return false; } } @@ -292,19 +265,140 @@ impl> SquareMatrix { true } - fn tr_solve_upper_triangular_vector_mut(&self, b: &mut Vector) -> bool - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, + /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only + /// the lower-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn ad_solve_lower_triangular( + &self, + b: &Matrix, + ) -> Option> + where + S2: StorageMut, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + if self.ad_solve_lower_triangular_mut(&mut res) { + Some(res) + } else { + None + } + } + + /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only + /// the upper-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn ad_solve_upper_triangular( + &self, + b: &Matrix, + ) -> Option> + where + S2: StorageMut, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + if self.ad_solve_upper_triangular_mut(&mut res) { + Some(res) + } else { + None + } + } + + /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the + /// lower-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn ad_solve_lower_triangular_mut( + &self, + b: &mut Matrix, + ) -> bool + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let cols = b.ncols(); + + for i in 0..cols { + if !self.xx_solve_lower_triangular_vector_mut(&mut b.column_mut(i), |e| e.conjugate(), |a, b| a.dotc(b)) { + return false; + } + } + + true + } + + /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the + /// upper-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn ad_solve_upper_triangular_mut( + &self, + b: &mut Matrix, + ) -> bool + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let cols = b.ncols(); + + for i in 0..cols { + if !self.xx_solve_upper_triangular_vector_mut(&mut b.column_mut(i), |e| e.conjugate(), |a, b| a.dotc(b)) { + return false; + } + } + + true + } + + + #[inline(always)] + fn xx_solve_lower_triangular_vector_mut( + &self, + b: &mut Vector, + conjugate: impl Fn(N) -> N, + dot: impl Fn(&DVectorSlice, &DVectorSlice) -> N, + ) -> bool + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let dim = self.nrows(); + + for i in (0..dim).rev() { + let dot = dot(&self.slice_range(i + 1.., i), &b.slice_range(i + 1.., 0)); + + unsafe { + let b_i = b.vget_unchecked_mut(i); + + let diag = conjugate(*self.get_unchecked((i, i))); + + if diag.is_zero() { + return false; + } + + *b_i = (*b_i - dot) / diag; + } + } + + true + } + + #[inline(always)] + fn xx_solve_upper_triangular_vector_mut( + &self, + b: &mut Vector, + conjugate: impl Fn(N) -> N, + dot: impl Fn(&DVectorSlice, &DVectorSlice) -> N, + ) -> bool + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, { let dim = self.nrows(); for i in 0..dim { - let dot = self.slice_range(..i, i).dot(&b.slice_range(..i, 0)); + let dot = dot(&self.slice_range(..i, i), &b.slice_range(..i, 0)); unsafe { let b_i = b.vget_unchecked_mut(i); - let diag = *self.get_unchecked((i, i)); + let diag = conjugate(*self.get_unchecked((i, i))); if diag.is_zero() { return false; diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index a455badd..12ee0469 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -496,7 +496,7 @@ where } } - self.recompose().map(|m| m.conjugate_transpose()) + self.recompose().map(|m| m.adjoint()) } } @@ -521,7 +521,7 @@ where else { match (&self.u, &self.v_t) { (Some(u), Some(v_t)) => { - let mut ut_b = u.conjugate().tr_mul(b); + let mut ut_b = u.ad_mul(b); for j in 0..ut_b.ncols() { let mut col = ut_b.column_mut(j); @@ -536,7 +536,7 @@ where } } - Ok(v_t.conjugate().tr_mul(&ut_b)) + Ok(v_t.ad_mul(&ut_b)) } (None, None) => Err("SVD solve: U and V^t have not been computed."), (None, _) => Err("SVD solve: U has not been computed."), diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index 29451d31..d652ce4c 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -75,13 +75,13 @@ where DefaultAllocator: Allocator + Allocator> if not_zero { let mut p = p.rows_range_mut(i..); - p.cgemv_symm(::convert(2.0), &m, &axis, N::zero()); - let dot = axis.cdot(&p); + p.hegemv(::convert(2.0), &m, &axis, N::zero()); + let dot = axis.dotc(&p); // p.axpy(-dot, &axis.conjugate(), N::one()); - m.ger_symm(-N::one(), &p, &axis.conjugate(), N::one()); - m.ger_symm(-N::one(), &axis, &p.conjugate(), N::one()); - m.ger_symm(dot * ::convert(2.0), &axis, &axis.conjugate(), N::one()); + m.hegerc(-N::one(), &p, &axis, N::one()); + m.hegerc(-N::one(), &axis, &p, N::one()); + m.hegerc(dot * ::convert(2.0), &axis, &axis, N::one()); } } @@ -142,7 +142,7 @@ where DefaultAllocator: Allocator + Allocator> self.tri[(i, i + 1)] = val; } - &q * self.tri * q.conjugate_transpose() + &q * self.tri * q.adjoint() } } diff --git a/tests/core/blas.rs b/tests/core/blas.rs index 00eac3a3..38113c17 100644 --- a/tests/core/blas.rs +++ b/tests/core/blas.rs @@ -19,14 +19,14 @@ quickcheck! { let mut y2 = y1.clone(); y1.gemv(alpha, &a, &x, beta); - y2.gemv_symm(alpha, &a.lower_triangle(), &x, beta); + y2.sygemv(alpha, &a.lower_triangle(), &x, beta); if !relative_eq!(y1, y2, epsilon = 1.0e-10) { return false; } y1.gemv(alpha, &a, &x, 0.0); - y2.gemv_symm(alpha, &a.lower_triangle(), &x, 0.0); + y2.sygemv(alpha, &a.lower_triangle(), &x, 0.0); relative_eq!(y1, y2, epsilon = 1.0e-10) } @@ -61,14 +61,14 @@ quickcheck! { let y = DVector::new_random(n); a1.ger(alpha, &x, &y, beta); - a2.ger_symm(alpha, &x, &y, beta); + a2.syger(alpha, &x, &y, beta); if !relative_eq!(a1.lower_triangle(), a2) { return false; } a1.ger(alpha, &x, &y, 0.0); - a2.ger_symm(alpha, &x, &y, 0.0); + a2.syger(alpha, &x, &y, 0.0); relative_eq!(a1.lower_triangle(), a2) } diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 5bc91dc2..774b76b0 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -16,7 +16,7 @@ macro_rules! gen_tests( fn cholesky(n: usize) -> bool { let m = RandomSDP::new(Dynamic::new(n.max(1).min(50)), || random::<$scalar>().0).unwrap(); let l = m.clone().cholesky().unwrap().unpack(); - relative_eq!(m, &l * l.conjugate_transpose(), epsilon = 1.0e-7) + relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7) } fn cholesky_static(_m: RandomSDP) -> bool { @@ -24,7 +24,7 @@ macro_rules! gen_tests( let chol = m.cholesky().unwrap(); let l = chol.unpack(); - if !relative_eq!(m, &l * l.conjugate_transpose(), epsilon = 1.0e-7) { + if !relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7) { false } else { diff --git a/tests/linalg/hessenberg.rs b/tests/linalg/hessenberg.rs index dfabd29d..03cba7d1 100644 --- a/tests/linalg/hessenberg.rs +++ b/tests/linalg/hessenberg.rs @@ -27,21 +27,21 @@ macro_rules! gen_tests( let hess = m.clone().hessenberg(); let (p, h) = hess.unpack(); - relative_eq!(m, &p * h * p.conjugate_transpose(), epsilon = 1.0e-7) + relative_eq!(m, &p * h * p.adjoint(), epsilon = 1.0e-7) } fn hessenberg_static_mat2(m: Matrix2<$scalar>) -> bool { let m = m.map(|e| e.0); let hess = m.hessenberg(); let (p, h) = hess.unpack(); - relative_eq!(m, p * h * p.conjugate_transpose(), epsilon = 1.0e-7) + relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7) } fn hessenberg_static(m: Matrix4<$scalar>) -> bool { let m = m.map(|e| e.0); let hess = m.hessenberg(); let (p, h) = hess.unpack(); - relative_eq!(m, p * h * p.conjugate_transpose(), epsilon = 1.0e-7) + relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7) } } } diff --git a/tests/linalg/schur.rs b/tests/linalg/schur.rs index 16767c94..f607ce06 100644 --- a/tests/linalg/schur.rs +++ b/tests/linalg/schur.rs @@ -31,21 +31,21 @@ mod quickcheck_tests { let (vecs, vals) = m.clone().schur().unpack(); - if !relative_eq!(&vecs * &vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7) { - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose()); + if !relative_eq!(&vecs * &vals * vecs.adjoint(), m, epsilon = 1.0e-7) { + println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint()); } - relative_eq!(&vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7) + relative_eq!(&vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7) } fn schur_static_mat2(m: Matrix2<$scalar>) -> bool { let m = m.map(|e| e.0); let (vecs, vals) = m.clone().schur().unpack(); - let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7); + let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7); if !ok { println!("Vecs: {:.5} Vals: {:.5}", vecs, vals); - println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.conjugate_transpose()); + println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.adjoint()); } ok } @@ -54,10 +54,10 @@ mod quickcheck_tests { let m = m.map(|e| e.0); let (vecs, vals) = m.clone().schur().unpack(); - let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7); + let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7); if !ok { println!("Vecs: {:.5} Vals: {:.5}", vecs, vals); - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose()); + println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint()); } ok } @@ -66,9 +66,9 @@ mod quickcheck_tests { let m = m.map(|e| e.0); let (vecs, vals) = m.clone().schur().unpack(); - let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7); + let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7); if !ok { - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose()); + println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint()); } ok diff --git a/tests/linalg/solve.rs b/tests/linalg/solve.rs index e2f8a3f8..3a9ff143 100644 --- a/tests/linalg/solve.rs +++ b/tests/linalg/solve.rs @@ -10,7 +10,7 @@ macro_rules! gen_tests( fn unzero_diagonal(a: &mut Matrix4) { for i in 0..4 { - if a[(i, i)].asum() < na::convert(1.0e-7) { + if a[(i, i)].norm1() < na::convert(1.0e-7) { a[(i, i)] = N::one(); } }