nalgebra/base/
ops.rs

1use num::{One, Zero};
2use std::iter;
3use std::ops::{
4    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
5};
6
7use simba::scalar::{
8    ClosedAddAssign, ClosedDivAssign, ClosedMulAssign, ClosedNeg, ClosedSubAssign,
9};
10
11use crate::base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR};
12use crate::base::blas_uninit::gemm_uninit;
13use crate::base::constraint::{
14    AreMultipliable, DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint,
15};
16use crate::base::dimension::{Dim, DimMul, DimName, DimProd, Dyn};
17use crate::base::storage::{Storage, StorageMut};
18use crate::base::uninit::Uninit;
19use crate::base::{DefaultAllocator, Matrix, MatrixSum, OMatrix, Scalar, VectorView};
20use crate::storage::IsContiguous;
21use crate::uninit::{Init, InitStatus};
22use crate::{RawStorage, RawStorageMut, SimdComplexField};
23use std::mem::MaybeUninit;
24
25/*
26 *
27 * Indexing.
28 *
29 */
30impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Index<usize> for Matrix<T, R, C, S> {
31    type Output = T;
32
33    #[inline]
34    fn index(&self, i: usize) -> &Self::Output {
35        let ij = self.vector_to_matrix_index(i);
36        &self[ij]
37    }
38}
39
40impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Index<(usize, usize)> for Matrix<T, R, C, S> {
41    type Output = T;
42
43    #[inline]
44    fn index(&self, ij: (usize, usize)) -> &Self::Output {
45        let shape = self.shape();
46        assert!(
47            ij.0 < shape.0 && ij.1 < shape.1,
48            "Matrix index out of bounds."
49        );
50
51        unsafe { self.get_unchecked((ij.0, ij.1)) }
52    }
53}
54
55// Mutable versions.
56impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> IndexMut<usize> for Matrix<T, R, C, S> {
57    #[inline]
58    fn index_mut(&mut self, i: usize) -> &mut T {
59        let ij = self.vector_to_matrix_index(i);
60        &mut self[ij]
61    }
62}
63
64impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> IndexMut<(usize, usize)> for Matrix<T, R, C, S> {
65    #[inline]
66    fn index_mut(&mut self, ij: (usize, usize)) -> &mut T {
67        let shape = self.shape();
68        assert!(
69            ij.0 < shape.0 && ij.1 < shape.1,
70            "Matrix index out of bounds."
71        );
72
73        unsafe { self.get_unchecked_mut((ij.0, ij.1)) }
74    }
75}
76
77/*
78 *
79 * Neg
80 *
81 */
82impl<T, R: Dim, C: Dim, S> Neg for Matrix<T, R, C, S>
83where
84    T: Scalar + ClosedNeg,
85    S: Storage<T, R, C>,
86    DefaultAllocator: Allocator<R, C>,
87{
88    type Output = OMatrix<T, R, C>;
89
90    #[inline]
91    fn neg(self) -> Self::Output {
92        let mut res = self.into_owned();
93        res.neg_mut();
94        res
95    }
96}
97
98impl<T, R: Dim, C: Dim, S> Neg for &Matrix<T, R, C, S>
99where
100    T: Scalar + ClosedNeg,
101    S: Storage<T, R, C>,
102    DefaultAllocator: Allocator<R, C>,
103{
104    type Output = OMatrix<T, R, C>;
105
106    #[inline]
107    fn neg(self) -> Self::Output {
108        -self.clone_owned()
109    }
110}
111
112impl<T, R: Dim, C: Dim, S> Matrix<T, R, C, S>
113where
114    T: Scalar + ClosedNeg,
115    S: StorageMut<T, R, C>,
116{
117    /// Negates `self` in-place.
118    #[inline]
119    pub fn neg_mut(&mut self) {
120        for e in self.iter_mut() {
121            *e = -e.clone()
122        }
123    }
124}
125
126/*
127 *
128 * Addition & Subtraction
129 *
130 */
131
132macro_rules! componentwise_binop_impl(
133    ($Trait: ident, $method: ident, $bound: ident;
134     $TraitAssign: ident, $method_assign: ident, $method_assign_statically_unchecked: ident,
135     $method_assign_statically_unchecked_rhs: ident;
136     $method_to: ident, $method_to_statically_unchecked_uninit: ident) => {
137
138        impl<T, R1: Dim, C1: Dim, SA: Storage<T, R1, C1>> Matrix<T, R1, C1, SA>
139            where T: Scalar + $bound {
140
141            /*
142             *
143             * Methods without dimension checking at compile-time.
144             * This is useful for code reuse because the sum representative system does not plays
145             * easily with static checks.
146             *
147             */
148            #[inline]
149            fn $method_to_statically_unchecked_uninit<Status, R2: Dim, C2: Dim, SB,
150                                                       R3: Dim, C3: Dim, SC>(&self,
151                                                                     _status: Status,
152                                                                     rhs: &Matrix<T, R2, C2, SB>,
153                                                                     out: &mut Matrix<Status::Value, R3, C3, SC>)
154                where Status: InitStatus<T>,
155                      SB: RawStorage<T, R2, C2>,
156                      SC: RawStorageMut<Status::Value, R3, C3> {
157                assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
158                assert_eq!(self.shape(), out.shape(), "Matrix addition/subtraction output dimensions mismatch.");
159
160                // This is the most common case and should be deduced at compile-time.
161                // TODO: use specialization instead?
162                unsafe {
163                    if self.data.is_contiguous() && rhs.data.is_contiguous() && out.data.is_contiguous() {
164                        let arr1 = self.data.as_slice_unchecked();
165                        let arr2 = rhs.data.as_slice_unchecked();
166                        let out  = out.data.as_mut_slice_unchecked();
167                        for i in 0 .. arr1.len() {
168                            Status::init(out.get_unchecked_mut(i), arr1.get_unchecked(i).clone().$method(arr2.get_unchecked(i).clone()));
169                        }
170                    } else {
171                        for j in 0 .. self.ncols() {
172                            for i in 0 .. self.nrows() {
173                                let val = self.get_unchecked((i, j)).clone().$method(rhs.get_unchecked((i, j)).clone());
174                                Status::init(out.get_unchecked_mut((i, j)), val);
175                            }
176                        }
177                    }
178                }
179            }
180
181
182            #[inline]
183            fn $method_assign_statically_unchecked<R2, C2, SB>(&mut self, rhs: &Matrix<T, R2, C2, SB>)
184                where R2: Dim,
185                      C2: Dim,
186                      SA: StorageMut<T, R1, C1>,
187                      SB: Storage<T, R2, C2> {
188                assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
189
190                // This is the most common case and should be deduced at compile-time.
191                // TODO: use specialization instead?
192                unsafe {
193                    if self.data.is_contiguous() && rhs.data.is_contiguous() {
194                        let arr1 = self.data.as_mut_slice_unchecked();
195                        let arr2 = rhs.data.as_slice_unchecked();
196
197                        for i in 0 .. arr2.len() {
198                            arr1.get_unchecked_mut(i).$method_assign(arr2.get_unchecked(i).clone());
199                        }
200                    } else {
201                        for j in 0 .. rhs.ncols() {
202                            for i in 0 .. rhs.nrows() {
203                                self.get_unchecked_mut((i, j)).$method_assign(rhs.get_unchecked((i, j)).clone())
204                            }
205                        }
206                    }
207                }
208            }
209
210
211            #[inline]
212            fn $method_assign_statically_unchecked_rhs<R2, C2, SB>(&self, rhs: &mut Matrix<T, R2, C2, SB>)
213                where R2: Dim,
214                      C2: Dim,
215                      SB: StorageMut<T, R2, C2> {
216                assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
217
218                // This is the most common case and should be deduced at compile-time.
219                // TODO: use specialization instead?
220                unsafe {
221                    if self.data.is_contiguous() && rhs.data.is_contiguous() {
222                        let arr1 = self.data.as_slice_unchecked();
223                        let arr2 = rhs.data.as_mut_slice_unchecked();
224
225                        for i in 0 .. arr1.len() {
226                            let res = arr1.get_unchecked(i).clone().$method(arr2.get_unchecked(i).clone());
227                            *arr2.get_unchecked_mut(i) = res;
228                        }
229                    } else {
230                        for j in 0 .. self.ncols() {
231                            for i in 0 .. self.nrows() {
232                                let r = rhs.get_unchecked_mut((i, j));
233                                *r = self.get_unchecked((i, j)).clone().$method(r.clone())
234                            }
235                        }
236                    }
237                }
238            }
239
240
241            /*
242             *
243             * Methods without dimension checking at compile-time.
244             * This is useful for code reuse because the sum representative system does not plays
245             * easily with static checks.
246             *
247             */
248            /// Equivalent to `self + rhs` but stores the result into `out` to avoid allocations.
249            #[inline]
250            pub fn $method_to<R2: Dim, C2: Dim, SB,
251                              R3: Dim, C3: Dim, SC>(&self,
252                                                    rhs: &Matrix<T, R2, C2, SB>,
253                                                    out: &mut Matrix<T, R3, C3, SC>)
254                where SB: Storage<T, R2, C2>,
255                      SC: StorageMut<T, R3, C3>,
256                      ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> +
257                                       SameNumberOfRows<R1, R3> + SameNumberOfColumns<C1, C3> {
258                self.$method_to_statically_unchecked_uninit(Init, rhs, out)
259            }
260        }
261
262        impl<'b, T, R1, C1, R2, C2, SA, SB> $Trait<&'b Matrix<T, R2, C2, SB>> for Matrix<T, R1, C1, SA>
263            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
264                  T: Scalar + $bound,
265                  SA: Storage<T, R1, C1>,
266                  SB: Storage<T, R2, C2>,
267                  DefaultAllocator: SameShapeAllocator<R1, C1, R2, C2>,
268                  ShapeConstraint:  SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
269            type Output = MatrixSum<T, R1, C1, R2, C2>;
270
271            #[inline]
272            fn $method(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
273                assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
274                let mut res = self.into_owned_sum::<R2, C2>();
275                res.$method_assign_statically_unchecked(rhs);
276                res
277            }
278        }
279
280        impl<'a, T, R1, C1, R2, C2, SA, SB> $Trait<Matrix<T, R2, C2, SB>> for &'a Matrix<T, R1, C1, SA>
281            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
282                  T: Scalar + $bound,
283                  SA: Storage<T, R1, C1>,
284                  SB: Storage<T, R2, C2>,
285                  DefaultAllocator: SameShapeAllocator<R2, C2, R1, C1>,
286                  ShapeConstraint:  SameNumberOfRows<R2, R1> + SameNumberOfColumns<C2, C1> {
287            type Output = MatrixSum<T, R2, C2, R1, C1>;
288
289            #[inline]
290            fn $method(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
291                let mut rhs = rhs.into_owned_sum::<R1, C1>();
292                assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch.");
293                self.$method_assign_statically_unchecked_rhs(&mut rhs);
294                rhs
295            }
296        }
297
298        impl<T, R1, C1, R2, C2, SA, SB> $Trait<Matrix<T, R2, C2, SB>> for Matrix<T, R1, C1, SA>
299            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
300                  T: Scalar + $bound,
301                  SA: Storage<T, R1, C1>,
302                  SB: Storage<T, R2, C2>,
303                  DefaultAllocator: SameShapeAllocator<R1, C1, R2, C2>,
304                  ShapeConstraint:  SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
305            type Output = MatrixSum<T, R1, C1, R2, C2>;
306
307            #[inline]
308            fn $method(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
309                self.$method(&rhs)
310            }
311        }
312
313        impl<'a, 'b, T, R1, C1, R2, C2, SA, SB> $Trait<&'b Matrix<T, R2, C2, SB>> for &'a Matrix<T, R1, C1, SA>
314            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
315                  T: Scalar + $bound,
316                  SA: Storage<T, R1, C1>,
317                  SB: Storage<T, R2, C2>,
318                  DefaultAllocator: SameShapeAllocator<R1, C1, R2, C2>,
319                  ShapeConstraint:  SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
320            type Output = MatrixSum<T, R1, C1, R2, C2>;
321
322            #[inline]
323            fn $method(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
324                let (nrows, ncols) = self.shape();
325                let nrows: SameShapeR<R1, R2> = Dim::from_usize(nrows);
326                let ncols: SameShapeC<C1, C2> = Dim::from_usize(ncols);
327                let mut res = Matrix::uninit(nrows, ncols);
328                self.$method_to_statically_unchecked_uninit(Uninit, rhs, &mut res);
329                // SAFETY: the output has been initialized above.
330                unsafe { res.assume_init() }
331            }
332        }
333
334        impl<'b, T, R1, C1, R2, C2, SA, SB> $TraitAssign<&'b Matrix<T, R2, C2, SB>> for Matrix<T, R1, C1, SA>
335            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
336                  T: Scalar + $bound,
337                  SA: StorageMut<T, R1, C1>,
338                  SB: Storage<T, R2, C2>,
339                  ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
340
341            #[inline]
342            fn $method_assign(&mut self, rhs: &'b Matrix<T, R2, C2, SB>) {
343                self.$method_assign_statically_unchecked(rhs)
344            }
345        }
346
347        impl<T, R1, C1, R2, C2, SA, SB> $TraitAssign<Matrix<T, R2, C2, SB>> for Matrix<T, R1, C1, SA>
348            where R1: Dim, C1: Dim, R2: Dim, C2: Dim,
349                  T: Scalar + $bound,
350                  SA: StorageMut<T, R1, C1>,
351                  SB: Storage<T, R2, C2>,
352                  ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
353
354            #[inline]
355            fn $method_assign(&mut self, rhs: Matrix<T, R2, C2, SB>) {
356                self.$method_assign(&rhs)
357            }
358        }
359    }
360);
361
362componentwise_binop_impl!(Add, add, ClosedAddAssign;
363                          AddAssign, add_assign, add_assign_statically_unchecked, add_assign_statically_unchecked_mut;
364                          add_to, add_to_statically_unchecked_uninit);
365componentwise_binop_impl!(Sub, sub, ClosedSubAssign;
366                          SubAssign, sub_assign, sub_assign_statically_unchecked, sub_assign_statically_unchecked_mut;
367                          sub_to, sub_to_statically_unchecked_uninit);
368
369impl<T, R: DimName, C: DimName> iter::Sum for OMatrix<T, R, C>
370where
371    T: Scalar + ClosedAddAssign + Zero,
372    DefaultAllocator: Allocator<R, C>,
373{
374    fn sum<I: Iterator<Item = OMatrix<T, R, C>>>(iter: I) -> OMatrix<T, R, C> {
375        iter.fold(Matrix::zero(), |acc, x| acc + x)
376    }
377}
378
379impl<T, C: Dim> iter::Sum for OMatrix<T, Dyn, C>
380where
381    T: Scalar + ClosedAddAssign + Zero,
382    DefaultAllocator: Allocator<Dyn, C>,
383{
384    /// # Example
385    /// ```
386    /// # use nalgebra::DVector;
387    /// assert_eq!(vec![DVector::repeat(3, 1.0f64),
388    ///                 DVector::repeat(3, 1.0f64),
389    ///                 DVector::repeat(3, 1.0f64)].into_iter().sum::<DVector<f64>>(),
390    ///            DVector::repeat(3, 1.0f64) + DVector::repeat(3, 1.0f64) + DVector::repeat(3, 1.0f64));
391    /// ```
392    ///
393    /// # Panics
394    /// Panics if the iterator is empty:
395    /// ```should_panic
396    /// # use std::iter;
397    /// # use nalgebra::DMatrix;
398    /// iter::empty::<DMatrix<f64>>().sum::<DMatrix<f64>>(); // panics!
399    /// ```
400    fn sum<I: Iterator<Item = OMatrix<T, Dyn, C>>>(mut iter: I) -> OMatrix<T, Dyn, C> {
401        match iter.next() {
402            Some(first) => iter.fold(first, |acc, x| acc + x),
403            None => {
404                panic!("Cannot compute `sum` of empty iterator.")
405            }
406        }
407    }
408}
409
410impl<'a, T, R: DimName, C: DimName> iter::Sum<&'a OMatrix<T, R, C>> for OMatrix<T, R, C>
411where
412    T: Scalar + ClosedAddAssign + Zero,
413    DefaultAllocator: Allocator<R, C>,
414{
415    fn sum<I: Iterator<Item = &'a OMatrix<T, R, C>>>(iter: I) -> OMatrix<T, R, C> {
416        iter.fold(Matrix::zero(), |acc, x| acc + x)
417    }
418}
419
420impl<'a, T, C: Dim> iter::Sum<&'a OMatrix<T, Dyn, C>> for OMatrix<T, Dyn, C>
421where
422    T: Scalar + ClosedAddAssign + Zero,
423    DefaultAllocator: Allocator<Dyn, C>,
424{
425    /// # Example
426    /// ```
427    /// # use nalgebra::DVector;
428    /// let v = &DVector::repeat(3, 1.0f64);
429    ///
430    /// assert_eq!(vec![v, v, v].into_iter().sum::<DVector<f64>>(),
431    ///            v + v + v);
432    /// ```
433    ///
434    /// # Panics
435    /// Panics if the iterator is empty:
436    /// ```should_panic
437    /// # use std::iter;
438    /// # use nalgebra::DMatrix;
439    /// iter::empty::<&DMatrix<f64>>().sum::<DMatrix<f64>>(); // panics!
440    /// ```
441    fn sum<I: Iterator<Item = &'a OMatrix<T, Dyn, C>>>(mut iter: I) -> OMatrix<T, Dyn, C> {
442        if let Some(first) = iter.next() {
443            iter.fold(first.clone(), |acc, x| acc + x)
444        } else {
445            panic!("Cannot compute `sum` of empty iterator.")
446        }
447    }
448}
449
450/*
451 *
452 * Multiplication
453 *
454 */
455
456// Matrix × Scalar
457// Matrix / Scalar
458macro_rules! componentwise_scalarop_impl(
459    ($Trait: ident, $method: ident, $bound: ident;
460     $TraitAssign: ident, $method_assign: ident) => {
461        impl<T, R: Dim, C: Dim, S> $Trait<T> for Matrix<T, R, C, S>
462            where T: Scalar + $bound,
463                  S: Storage<T, R, C>,
464                  DefaultAllocator: Allocator<R, C> {
465            type Output = OMatrix<T, R, C>;
466
467            #[inline]
468            fn $method(self, rhs: T) -> Self::Output {
469                let mut res = self.into_owned();
470
471                // XXX: optimize our iterator!
472                //
473                // Using our own iterator prevents loop unrolling, which breaks some optimization
474                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
475
476                // for left in res.iter_mut() {
477                for left in res.as_mut_slice().iter_mut() {
478                    *left = left.clone().$method(rhs.clone())
479                }
480
481                res
482            }
483        }
484
485        impl<'a, T, R: Dim, C: Dim, S> $Trait<T> for &'a Matrix<T, R, C, S>
486            where T: Scalar + $bound,
487                  S: Storage<T, R, C>,
488                  DefaultAllocator: Allocator<R, C> {
489            type Output = OMatrix<T, R, C>;
490
491            #[inline]
492            fn $method(self, rhs: T) -> Self::Output {
493                self.clone_owned().$method(rhs)
494            }
495        }
496
497        impl<T, R: Dim, C: Dim, S> $TraitAssign<T> for Matrix<T, R, C, S>
498            where T: Scalar + $bound,
499                  S: StorageMut<T, R, C> {
500            #[inline]
501            fn $method_assign(&mut self, rhs: T) {
502                for j in 0 .. self.ncols() {
503                    for i in 0 .. self.nrows() {
504                        unsafe { self.get_unchecked_mut((i, j)).$method_assign(rhs.clone()) };
505                    }
506                }
507            }
508        }
509    }
510);
511
512componentwise_scalarop_impl!(Mul, mul, ClosedMulAssign; MulAssign, mul_assign);
513componentwise_scalarop_impl!(Div, div, ClosedDivAssign; DivAssign, div_assign);
514
515macro_rules! left_scalar_mul_impl(
516    ($($T: ty),* $(,)*) => {$(
517        impl<R: Dim, C: Dim, S: Storage<$T, R, C>> Mul<Matrix<$T, R, C, S>> for $T
518            where DefaultAllocator: Allocator<R, C> {
519            type Output = OMatrix<$T, R, C>;
520
521            #[inline]
522            fn mul(self, rhs: Matrix<$T, R, C, S>) -> Self::Output {
523                let mut res = rhs.into_owned();
524
525                // XXX: optimize our iterator!
526                //
527                // Using our own iterator prevents loop unrolling, which breaks some optimization
528                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
529
530                // for rhs in res.iter_mut() {
531                for rhs in res.as_mut_slice().iter_mut() {
532                    *rhs *= self
533                }
534
535                res
536            }
537        }
538
539        impl<'b, R: Dim, C: Dim, S: Storage<$T, R, C>> Mul<&'b Matrix<$T, R, C, S>> for $T
540            where DefaultAllocator: Allocator<R, C> {
541            type Output = OMatrix<$T, R, C>;
542
543            #[inline]
544            fn mul(self, rhs: &'b Matrix<$T, R, C, S>) -> Self::Output {
545                self * rhs.clone_owned()
546            }
547        }
548    )*}
549);
550
551left_scalar_mul_impl!(u8, u16, u32, u64, usize, i8, i16, i32, i64, isize, f32, f64);
552
553// Matrix × Matrix
554impl<'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<&'b Matrix<T, R2, C2, SB>>
555    for &Matrix<T, R1, C1, SA>
556where
557    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
558    SA: Storage<T, R1, C1>,
559    SB: Storage<T, R2, C2>,
560    DefaultAllocator: Allocator<R1, C2>,
561    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
562{
563    type Output = OMatrix<T, R1, C2>;
564
565    #[inline]
566    fn mul(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
567        let mut res = Matrix::uninit(self.shape_generic().0, rhs.shape_generic().1);
568        unsafe {
569            // SAFETY: this is OK because status = Uninit && bevy == 0
570            gemm_uninit(Uninit, &mut res, T::one(), self, rhs, T::zero());
571            res.assume_init()
572        }
573    }
574}
575
576impl<T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<Matrix<T, R2, C2, SB>>
577    for &Matrix<T, R1, C1, SA>
578where
579    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
580    SB: Storage<T, R2, C2>,
581    SA: Storage<T, R1, C1>,
582    DefaultAllocator: Allocator<R1, C2>,
583    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
584{
585    type Output = OMatrix<T, R1, C2>;
586
587    #[inline]
588    fn mul(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
589        self * &rhs
590    }
591}
592
593impl<'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<&'b Matrix<T, R2, C2, SB>>
594    for Matrix<T, R1, C1, SA>
595where
596    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
597    SB: Storage<T, R2, C2>,
598    SA: Storage<T, R1, C1>,
599    DefaultAllocator: Allocator<R1, C2>,
600    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
601{
602    type Output = OMatrix<T, R1, C2>;
603
604    #[inline]
605    fn mul(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
606        &self * rhs
607    }
608}
609
610impl<T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<Matrix<T, R2, C2, SB>>
611    for Matrix<T, R1, C1, SA>
612where
613    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
614    SB: Storage<T, R2, C2>,
615    SA: Storage<T, R1, C1>,
616    DefaultAllocator: Allocator<R1, C2>,
617    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
618{
619    type Output = OMatrix<T, R1, C2>;
620
621    #[inline]
622    fn mul(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
623        &self * &rhs
624    }
625}
626
627// TODO: this is too restrictive:
628//    − we can't use `a *= b` when `a` is a mutable slice.
629//    − we can't use `a *= b` when C2 is not equal to C1.
630impl<T, R1, C1, R2, SA, SB> MulAssign<Matrix<T, R2, C1, SB>> for Matrix<T, R1, C1, SA>
631where
632    R1: Dim,
633    C1: Dim,
634    R2: Dim,
635    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
636    SB: Storage<T, R2, C1>,
637    SA: StorageMut<T, R1, C1> + IsContiguous + Clone, // TODO: get rid of the IsContiguous
638    ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
639    DefaultAllocator: Allocator<R1, C1, Buffer<T> = SA>,
640{
641    #[inline]
642    fn mul_assign(&mut self, rhs: Matrix<T, R2, C1, SB>) {
643        *self = &*self * rhs
644    }
645}
646
647impl<'b, T, R1, C1, R2, SA, SB> MulAssign<&'b Matrix<T, R2, C1, SB>> for Matrix<T, R1, C1, SA>
648where
649    R1: Dim,
650    C1: Dim,
651    R2: Dim,
652    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
653    SB: Storage<T, R2, C1>,
654    SA: StorageMut<T, R1, C1> + IsContiguous + Clone, // TODO: get rid of the IsContiguous
655    ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
656    // TODO: this is too restrictive. See comments for the non-ref version.
657    DefaultAllocator: Allocator<R1, C1, Buffer<T> = SA>,
658{
659    #[inline]
660    fn mul_assign(&mut self, rhs: &'b Matrix<T, R2, C1, SB>) {
661        *self = &*self * rhs
662    }
663}
664
665/// # Special multiplications.
666impl<T, R1: Dim, C1: Dim, SA> Matrix<T, R1, C1, SA>
667where
668    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
669    SA: Storage<T, R1, C1>,
670{
671    /// Equivalent to `self.transpose() * rhs`.
672    #[inline]
673    #[must_use]
674    pub fn tr_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
675    where
676        SB: Storage<T, R2, C2>,
677        DefaultAllocator: Allocator<C1, C2>,
678        ShapeConstraint: SameNumberOfRows<R1, R2>,
679    {
680        let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1);
681        self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dot(b));
682        // SAFETY: this is OK because the result is now initialized.
683        unsafe { res.assume_init() }
684    }
685
686    /// Equivalent to `self.adjoint() * rhs`.
687    #[inline]
688    #[must_use]
689    pub fn ad_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
690    where
691        T: SimdComplexField,
692        SB: Storage<T, R2, C2>,
693        DefaultAllocator: Allocator<C1, C2>,
694        ShapeConstraint: SameNumberOfRows<R1, R2>,
695    {
696        let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1);
697        self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dotc(b));
698        // SAFETY: this is OK because the result is now initialized.
699        unsafe { res.assume_init() }
700    }
701
702    #[inline(always)]
703    fn xx_mul_to_uninit<Status, R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
704        &self,
705        _status: Status,
706        rhs: &Matrix<T, R2, C2, SB>,
707        out: &mut Matrix<Status::Value, R3, C3, SC>,
708        dot: impl Fn(
709            &VectorView<'_, T, R1, SA::RStride, SA::CStride>,
710            &VectorView<'_, T, R2, SB::RStride, SB::CStride>,
711        ) -> T,
712    ) where
713        Status: InitStatus<T>,
714        SB: RawStorage<T, R2, C2>,
715        SC: RawStorageMut<Status::Value, R3, C3>,
716        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
717    {
718        let (nrows1, ncols1) = self.shape();
719        let (nrows2, ncols2) = rhs.shape();
720        let (nrows3, ncols3) = out.shape();
721
722        assert!(
723            nrows1 == nrows2,
724            "Matrix multiplication dimensions mismatch {:?} and {:?}: left rows != right rows.",
725            self.shape(),
726            rhs.shape()
727        );
728        assert!(
729            ncols1 == nrows3,
730            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right rows.",
731            self.shape(),
732            out.shape()
733        );
734        assert!(
735            ncols2 == ncols3,
736            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right cols",
737            rhs.shape(),
738            out.shape()
739        );
740
741        for i in 0..ncols1 {
742            for j in 0..ncols2 {
743                let dot = dot(&self.column(i), &rhs.column(j));
744                let elt = unsafe { out.get_unchecked_mut((i, j)) };
745                Status::init(elt, dot);
746            }
747        }
748    }
749
750    /// Equivalent to `self.transpose() * rhs` but stores the result into `out` to avoid
751    /// allocations.
752    #[inline]
753    pub fn tr_mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
754        &self,
755        rhs: &Matrix<T, R2, C2, SB>,
756        out: &mut Matrix<T, R3, C3, SC>,
757    ) where
758        SB: Storage<T, R2, C2>,
759        SC: StorageMut<T, R3, C3>,
760        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
761    {
762        self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dot(b))
763    }
764
765    /// Equivalent to `self.adjoint() * rhs` but stores the result into `out` to avoid
766    /// allocations.
767    #[inline]
768    pub fn ad_mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
769        &self,
770        rhs: &Matrix<T, R2, C2, SB>,
771        out: &mut Matrix<T, R3, C3, SC>,
772    ) where
773        T: SimdComplexField,
774        SB: Storage<T, R2, C2>,
775        SC: StorageMut<T, R3, C3>,
776        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
777    {
778        self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dotc(b))
779    }
780
781    /// Equivalent to `self * rhs` but stores the result into `out` to avoid allocations.
782    #[inline]
783    pub fn mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
784        &self,
785        rhs: &Matrix<T, R2, C2, SB>,
786        out: &mut Matrix<T, R3, C3, SC>,
787    ) where
788        SB: Storage<T, R2, C2>,
789        SC: StorageMut<T, R3, C3>,
790        ShapeConstraint: SameNumberOfRows<R3, R1>
791            + SameNumberOfColumns<C3, C2>
792            + AreMultipliable<R1, C1, R2, C2>,
793    {
794        out.gemm(T::one(), self, rhs, T::zero());
795    }
796
797    /// The kronecker product of two matrices (aka. tensor product of the corresponding linear
798    /// maps).
799    #[must_use]
800    pub fn kronecker<R2: Dim, C2: Dim, SB>(
801        &self,
802        rhs: &Matrix<T, R2, C2, SB>,
803    ) -> OMatrix<T, DimProd<R1, R2>, DimProd<C1, C2>>
804    where
805        T: ClosedMulAssign,
806        R1: DimMul<R2>,
807        C1: DimMul<C2>,
808        SB: Storage<T, R2, C2>,
809        DefaultAllocator: Allocator<DimProd<R1, R2>, DimProd<C1, C2>>,
810    {
811        let (nrows1, ncols1) = self.shape_generic();
812        let (nrows2, ncols2) = rhs.shape_generic();
813
814        let mut res = Matrix::uninit(nrows1.mul(nrows2), ncols1.mul(ncols2));
815        let mut data_res = res.data.ptr_mut();
816
817        unsafe {
818            for j1 in 0..ncols1.value() {
819                for j2 in 0..ncols2.value() {
820                    for i1 in 0..nrows1.value() {
821                        let coeff = self.get_unchecked((i1, j1)).clone();
822
823                        for i2 in 0..nrows2.value() {
824                            *data_res = MaybeUninit::new(
825                                coeff.clone() * rhs.get_unchecked((i2, j2)).clone(),
826                            );
827                            data_res = data_res.offset(1);
828                        }
829                    }
830                }
831            }
832
833            // SAFETY: the result matrix has been initialized by the loop above.
834            res.assume_init()
835        }
836    }
837}
838
839impl<T, D: DimName> iter::Product for OMatrix<T, D, D>
840where
841    T: Scalar + Zero + One + ClosedMulAssign + ClosedAddAssign,
842    DefaultAllocator: Allocator<D, D>,
843{
844    fn product<I: Iterator<Item = OMatrix<T, D, D>>>(iter: I) -> OMatrix<T, D, D> {
845        iter.fold(Matrix::one(), |acc, x| acc * x)
846    }
847}
848
849impl<'a, T, D: DimName> iter::Product<&'a OMatrix<T, D, D>> for OMatrix<T, D, D>
850where
851    T: Scalar + Zero + One + ClosedMulAssign + ClosedAddAssign,
852    DefaultAllocator: Allocator<D, D>,
853{
854    fn product<I: Iterator<Item = &'a OMatrix<T, D, D>>>(iter: I) -> OMatrix<T, D, D> {
855        iter.fold(Matrix::one(), |acc, x| acc * x)
856    }
857}
nalgebra/base/ops.rs

nalgebra/base/
ops.rs
