nalgebra/base/
ops.rs

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

nalgebra/base/
ops.rs
