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<'a, T, R: Dim, C: Dim, S> Neg for &'a 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        if let Some(first) = iter.next() {
402            iter.fold(first, |acc, x| acc + x)
403        } else {
404            panic!("Cannot compute `sum` of empty iterator.")
405        }
406    }
407}
408
409impl<'a, T, R: DimName, C: DimName> iter::Sum<&'a OMatrix<T, R, C>> for OMatrix<T, R, C>
410where
411    T: Scalar + ClosedAddAssign + Zero,
412    DefaultAllocator: Allocator<R, C>,
413{
414    fn sum<I: Iterator<Item = &'a OMatrix<T, R, C>>>(iter: I) -> OMatrix<T, R, C> {
415        iter.fold(Matrix::zero(), |acc, x| acc + x)
416    }
417}
418
419impl<'a, T, C: Dim> iter::Sum<&'a OMatrix<T, Dyn, C>> for OMatrix<T, Dyn, C>
420where
421    T: Scalar + ClosedAddAssign + Zero,
422    DefaultAllocator: Allocator<Dyn, C>,
423{
424    /// # Example
425    /// ```
426    /// # use nalgebra::DVector;
427    /// let v = &DVector::repeat(3, 1.0f64);
428    ///
429    /// assert_eq!(vec![v, v, v].into_iter().sum::<DVector<f64>>(),
430    ///            v + v + v);
431    /// ```
432    ///
433    /// # Panics
434    /// Panics if the iterator is empty:
435    /// ```should_panic
436    /// # use std::iter;
437    /// # use nalgebra::DMatrix;
438    /// iter::empty::<&DMatrix<f64>>().sum::<DMatrix<f64>>(); // panics!
439    /// ```
440    fn sum<I: Iterator<Item = &'a OMatrix<T, Dyn, C>>>(mut iter: I) -> OMatrix<T, Dyn, C> {
441        if let Some(first) = iter.next() {
442            iter.fold(first.clone(), |acc, x| acc + x)
443        } else {
444            panic!("Cannot compute `sum` of empty iterator.")
445        }
446    }
447}
448
449/*
450 *
451 * Multiplication
452 *
453 */
454
455// Matrix × Scalar
456// Matrix / Scalar
457macro_rules! componentwise_scalarop_impl(
458    ($Trait: ident, $method: ident, $bound: ident;
459     $TraitAssign: ident, $method_assign: ident) => {
460        impl<T, R: Dim, C: Dim, S> $Trait<T> for Matrix<T, R, C, S>
461            where T: Scalar + $bound,
462                  S: Storage<T, R, C>,
463                  DefaultAllocator: Allocator<R, C> {
464            type Output = OMatrix<T, R, C>;
465
466            #[inline]
467            fn $method(self, rhs: T) -> Self::Output {
468                let mut res = self.into_owned();
469
470                // XXX: optimize our iterator!
471                //
472                // Using our own iterator prevents loop unrolling, which breaks some optimization
473                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
474
475                // for left in res.iter_mut() {
476                for left in res.as_mut_slice().iter_mut() {
477                    *left = left.clone().$method(rhs.clone())
478                }
479
480                res
481            }
482        }
483
484        impl<'a, T, R: Dim, C: Dim, S> $Trait<T> for &'a Matrix<T, R, C, S>
485            where T: Scalar + $bound,
486                  S: Storage<T, R, C>,
487                  DefaultAllocator: Allocator<R, C> {
488            type Output = OMatrix<T, R, C>;
489
490            #[inline]
491            fn $method(self, rhs: T) -> Self::Output {
492                self.clone_owned().$method(rhs)
493            }
494        }
495
496        impl<T, R: Dim, C: Dim, S> $TraitAssign<T> for Matrix<T, R, C, S>
497            where T: Scalar + $bound,
498                  S: StorageMut<T, R, C> {
499            #[inline]
500            fn $method_assign(&mut self, rhs: T) {
501                for j in 0 .. self.ncols() {
502                    for i in 0 .. self.nrows() {
503                        unsafe { self.get_unchecked_mut((i, j)).$method_assign(rhs.clone()) };
504                    }
505                }
506            }
507        }
508    }
509);
510
511componentwise_scalarop_impl!(Mul, mul, ClosedMulAssign; MulAssign, mul_assign);
512componentwise_scalarop_impl!(Div, div, ClosedDivAssign; DivAssign, div_assign);
513
514macro_rules! left_scalar_mul_impl(
515    ($($T: ty),* $(,)*) => {$(
516        impl<R: Dim, C: Dim, S: Storage<$T, R, C>> Mul<Matrix<$T, R, C, S>> for $T
517            where DefaultAllocator: Allocator<R, C> {
518            type Output = OMatrix<$T, R, C>;
519
520            #[inline]
521            fn mul(self, rhs: Matrix<$T, R, C, S>) -> Self::Output {
522                let mut res = rhs.into_owned();
523
524                // XXX: optimize our iterator!
525                //
526                // Using our own iterator prevents loop unrolling, which breaks some optimization
527                // (like SIMD). On the other hand, using the slice iterator is 4x faster.
528
529                // for rhs in res.iter_mut() {
530                for rhs in res.as_mut_slice().iter_mut() {
531                    *rhs *= self
532                }
533
534                res
535            }
536        }
537
538        impl<'b, R: Dim, C: Dim, S: Storage<$T, R, C>> Mul<&'b Matrix<$T, R, C, S>> for $T
539            where DefaultAllocator: Allocator<R, C> {
540            type Output = OMatrix<$T, R, C>;
541
542            #[inline]
543            fn mul(self, rhs: &'b Matrix<$T, R, C, S>) -> Self::Output {
544                self * rhs.clone_owned()
545            }
546        }
547    )*}
548);
549
550left_scalar_mul_impl!(u8, u16, u32, u64, usize, i8, i16, i32, i64, isize, f32, f64);
551
552// Matrix × Matrix
553impl<'a, 'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<&'b Matrix<T, R2, C2, SB>>
554    for &'a Matrix<T, R1, C1, SA>
555where
556    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
557    SA: Storage<T, R1, C1>,
558    SB: Storage<T, R2, C2>,
559    DefaultAllocator: Allocator<R1, C2>,
560    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
561{
562    type Output = OMatrix<T, R1, C2>;
563
564    #[inline]
565    fn mul(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
566        let mut res = Matrix::uninit(self.shape_generic().0, rhs.shape_generic().1);
567        unsafe {
568            // SAFETY: this is OK because status = Uninit && bevy == 0
569            gemm_uninit(Uninit, &mut res, T::one(), self, rhs, T::zero());
570            res.assume_init()
571        }
572    }
573}
574
575impl<'a, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<Matrix<T, R2, C2, SB>>
576    for &'a Matrix<T, R1, C1, SA>
577where
578    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
579    SB: Storage<T, R2, C2>,
580    SA: Storage<T, R1, C1>,
581    DefaultAllocator: Allocator<R1, C2>,
582    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
583{
584    type Output = OMatrix<T, R1, C2>;
585
586    #[inline]
587    fn mul(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
588        self * &rhs
589    }
590}
591
592impl<'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<&'b Matrix<T, R2, C2, SB>>
593    for Matrix<T, R1, C1, SA>
594where
595    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
596    SB: Storage<T, R2, C2>,
597    SA: Storage<T, R1, C1>,
598    DefaultAllocator: Allocator<R1, C2>,
599    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
600{
601    type Output = OMatrix<T, R1, C2>;
602
603    #[inline]
604    fn mul(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
605        &self * rhs
606    }
607}
608
609impl<T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> Mul<Matrix<T, R2, C2, SB>>
610    for Matrix<T, R1, C1, SA>
611where
612    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
613    SB: Storage<T, R2, C2>,
614    SA: Storage<T, R1, C1>,
615    DefaultAllocator: Allocator<R1, C2>,
616    ShapeConstraint: AreMultipliable<R1, C1, R2, C2>,
617{
618    type Output = OMatrix<T, R1, C2>;
619
620    #[inline]
621    fn mul(self, rhs: Matrix<T, R2, C2, SB>) -> Self::Output {
622        &self * &rhs
623    }
624}
625
626// TODO: this is too restrictive:
627//    − we can't use `a *= b` when `a` is a mutable slice.
628//    − we can't use `a *= b` when C2 is not equal to C1.
629impl<T, R1, C1, R2, SA, SB> MulAssign<Matrix<T, R2, C1, SB>> for Matrix<T, R1, C1, SA>
630where
631    R1: Dim,
632    C1: Dim,
633    R2: Dim,
634    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
635    SB: Storage<T, R2, C1>,
636    SA: StorageMut<T, R1, C1> + IsContiguous + Clone, // TODO: get rid of the IsContiguous
637    ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
638    DefaultAllocator: Allocator<R1, C1, Buffer<T> = SA>,
639{
640    #[inline]
641    fn mul_assign(&mut self, rhs: Matrix<T, R2, C1, SB>) {
642        *self = &*self * rhs
643    }
644}
645
646impl<'b, T, R1, C1, R2, SA, SB> MulAssign<&'b Matrix<T, R2, C1, SB>> for Matrix<T, R1, C1, SA>
647where
648    R1: Dim,
649    C1: Dim,
650    R2: Dim,
651    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
652    SB: Storage<T, R2, C1>,
653    SA: StorageMut<T, R1, C1> + IsContiguous + Clone, // TODO: get rid of the IsContiguous
654    ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
655    // TODO: this is too restrictive. See comments for the non-ref version.
656    DefaultAllocator: Allocator<R1, C1, Buffer<T> = SA>,
657{
658    #[inline]
659    fn mul_assign(&mut self, rhs: &'b Matrix<T, R2, C1, SB>) {
660        *self = &*self * rhs
661    }
662}
663
664/// # Special multiplications.
665impl<T, R1: Dim, C1: Dim, SA> Matrix<T, R1, C1, SA>
666where
667    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
668    SA: Storage<T, R1, C1>,
669{
670    /// Equivalent to `self.transpose() * rhs`.
671    #[inline]
672    #[must_use]
673    pub fn tr_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
674    where
675        SB: Storage<T, R2, C2>,
676        DefaultAllocator: Allocator<C1, C2>,
677        ShapeConstraint: SameNumberOfRows<R1, R2>,
678    {
679        let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1);
680        self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dot(b));
681        // SAFETY: this is OK because the result is now initialized.
682        unsafe { res.assume_init() }
683    }
684
685    /// Equivalent to `self.adjoint() * rhs`.
686    #[inline]
687    #[must_use]
688    pub fn ad_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
689    where
690        T: SimdComplexField,
691        SB: Storage<T, R2, C2>,
692        DefaultAllocator: Allocator<C1, C2>,
693        ShapeConstraint: SameNumberOfRows<R1, R2>,
694    {
695        let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1);
696        self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dotc(b));
697        // SAFETY: this is OK because the result is now initialized.
698        unsafe { res.assume_init() }
699    }
700
701    #[inline(always)]
702    fn xx_mul_to_uninit<Status, R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
703        &self,
704        _status: Status,
705        rhs: &Matrix<T, R2, C2, SB>,
706        out: &mut Matrix<Status::Value, R3, C3, SC>,
707        dot: impl Fn(
708            &VectorView<'_, T, R1, SA::RStride, SA::CStride>,
709            &VectorView<'_, T, R2, SB::RStride, SB::CStride>,
710        ) -> T,
711    ) where
712        Status: InitStatus<T>,
713        SB: RawStorage<T, R2, C2>,
714        SC: RawStorageMut<Status::Value, R3, C3>,
715        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
716    {
717        let (nrows1, ncols1) = self.shape();
718        let (nrows2, ncols2) = rhs.shape();
719        let (nrows3, ncols3) = out.shape();
720
721        assert!(
722            nrows1 == nrows2,
723            "Matrix multiplication dimensions mismatch {:?} and {:?}: left rows != right rows.",
724            self.shape(),
725            rhs.shape()
726        );
727        assert!(
728            ncols1 == nrows3,
729            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right rows.",
730            self.shape(),
731            out.shape()
732        );
733        assert!(
734            ncols2 == ncols3,
735            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right cols",
736            rhs.shape(),
737            out.shape()
738        );
739
740        for i in 0..ncols1 {
741            for j in 0..ncols2 {
742                let dot = dot(&self.column(i), &rhs.column(j));
743                let elt = unsafe { out.get_unchecked_mut((i, j)) };
744                Status::init(elt, dot);
745            }
746        }
747    }
748
749    /// Equivalent to `self.transpose() * rhs` but stores the result into `out` to avoid
750    /// allocations.
751    #[inline]
752    pub fn tr_mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
753        &self,
754        rhs: &Matrix<T, R2, C2, SB>,
755        out: &mut Matrix<T, R3, C3, SC>,
756    ) where
757        SB: Storage<T, R2, C2>,
758        SC: StorageMut<T, R3, C3>,
759        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
760    {
761        self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dot(b))
762    }
763
764    /// Equivalent to `self.adjoint() * rhs` but stores the result into `out` to avoid
765    /// allocations.
766    #[inline]
767    pub fn ad_mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
768        &self,
769        rhs: &Matrix<T, R2, C2, SB>,
770        out: &mut Matrix<T, R3, C3, SC>,
771    ) where
772        T: SimdComplexField,
773        SB: Storage<T, R2, C2>,
774        SC: StorageMut<T, R3, C3>,
775        ShapeConstraint: SameNumberOfRows<R1, R2> + DimEq<C1, R3> + DimEq<C2, C3>,
776    {
777        self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dotc(b))
778    }
779
780    /// Equivalent to `self * rhs` but stores the result into `out` to avoid allocations.
781    #[inline]
782    pub fn mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
783        &self,
784        rhs: &Matrix<T, R2, C2, SB>,
785        out: &mut Matrix<T, R3, C3, SC>,
786    ) where
787        SB: Storage<T, R2, C2>,
788        SC: StorageMut<T, R3, C3>,
789        ShapeConstraint: SameNumberOfRows<R3, R1>
790            + SameNumberOfColumns<C3, C2>
791            + AreMultipliable<R1, C1, R2, C2>,
792    {
793        out.gemm(T::one(), self, rhs, T::zero());
794    }
795
796    /// The kronecker product of two matrices (aka. tensor product of the corresponding linear
797    /// maps).
798    #[must_use]
799    pub fn kronecker<R2: Dim, C2: Dim, SB>(
800        &self,
801        rhs: &Matrix<T, R2, C2, SB>,
802    ) -> OMatrix<T, DimProd<R1, R2>, DimProd<C1, C2>>
803    where
804        T: ClosedMulAssign,
805        R1: DimMul<R2>,
806        C1: DimMul<C2>,
807        SB: Storage<T, R2, C2>,
808        DefaultAllocator: Allocator<DimProd<R1, R2>, DimProd<C1, C2>>,
809    {
810        let (nrows1, ncols1) = self.shape_generic();
811        let (nrows2, ncols2) = rhs.shape_generic();
812
813        let mut res = Matrix::uninit(nrows1.mul(nrows2), ncols1.mul(ncols2));
814        let mut data_res = res.data.ptr_mut();
815
816        unsafe {
817            for j1 in 0..ncols1.value() {
818                for j2 in 0..ncols2.value() {
819                    for i1 in 0..nrows1.value() {
820                        let coeff = self.get_unchecked((i1, j1)).clone();
821
822                        for i2 in 0..nrows2.value() {
823                            *data_res = MaybeUninit::new(
824                                coeff.clone() * rhs.get_unchecked((i2, j2)).clone(),
825                            );
826                            data_res = data_res.offset(1);
827                        }
828                    }
829                }
830            }
831
832            // SAFETY: the result matrix has been initialized by the loop above.
833            res.assume_init()
834        }
835    }
836}
837
838impl<T, D: DimName> iter::Product for OMatrix<T, D, D>
839where
840    T: Scalar + Zero + One + ClosedMulAssign + ClosedAddAssign,
841    DefaultAllocator: Allocator<D, D>,
842{
843    fn product<I: Iterator<Item = OMatrix<T, D, D>>>(iter: I) -> OMatrix<T, D, D> {
844        iter.fold(Matrix::one(), |acc, x| acc * x)
845    }
846}
847
848impl<'a, T, D: DimName> iter::Product<&'a OMatrix<T, D, D>> for OMatrix<T, D, D>
849where
850    T: Scalar + Zero + One + ClosedMulAssign + ClosedAddAssign,
851    DefaultAllocator: Allocator<D, D>,
852{
853    fn product<I: Iterator<Item = &'a OMatrix<T, D, D>>>(iter: I) -> OMatrix<T, D, D> {
854        iter.fold(Matrix::one(), |acc, x| acc * x)
855    }
856}
nalgebra/base/ops.rs

nalgebra/base/
ops.rs
