proof of concept for mul_tr

src/base/matrix_view.rs

@@ -546,17 +546,17 @@ macro_rules! matrix_view_impl (
             $me.$generic_view_with_steps(start, shape, steps)
         }
 
-        /// Slices this matrix starting at its component `(irow, icol)` and with `(R::dim(),
-        /// CView::dim())` consecutive components.
+        /// Slices this matrix starting at its component `(irow, icol)` and with `(RVIEW, CVIEW)`
+        /// consecutive components.
         #[inline]
         #[deprecated = slice_deprecation_note!($fixed_view)]
         pub fn $fixed_slice<const RVIEW: usize, const CVIEW: usize>($me: $Me, irow: usize, icol: usize)
             -> $MatrixView<'_, T, Const<RVIEW>, Const<CVIEW>, S::RStride, S::CStride> {
             $me.$fixed_view(irow, icol)
         }
 
-        /// Return a view of this matrix starting at its component `(irow, icol)` and with `(R::dim(),
-        /// CView::dim())` consecutive components.
+        /// Return a view of this matrix starting at its component `(irow, icol)` and with
+        /// `(RVIEW, CVIEW)` consecutive components.
         #[inline]
         pub fn $fixed_view<const RVIEW: usize, const CVIEW: usize>($me: $Me, irow: usize, icol: usize)
             -> $MatrixView<'_, T, Const<RVIEW>, Const<CVIEW>, S::RStride, S::CStride> {

src/base/ops.rs

@@ -679,6 +679,20 @@ where
         // SAFETY: this is OK because the result is now initialized.
         unsafe { res.assume_init() }
     }
+
+    #[inline]
+    #[must_use]
+    pub fn mul_tr<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
+    where
+        SB: Storage<T, R2, C2>,
+        DefaultAllocator: Allocator<T, R1, R2>,
+        ShapeConstraint: SameNumberOfColumns<C1, C2>,
+    {
+        let mut res = Matrix::uninit(self.shape_generic().0, rhs.shape_generic().0);
+        self.yy_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dot(b)); //note: this was changed
+        // SAFETY: this is OK because the result is now initialized.
+        unsafe { res.assume_init() }
+    }
 
     /// Equivalent to `self.adjoint() * rhs`.
     #[inline]
@@ -744,6 +758,54 @@ where
         }
     }
 
+     #[inline(always)]
+    fn yy_mul_to_uninit<Status, R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
+        &self,
+        _status: Status,
+        rhs: &Matrix<T, R2, C2, SB>,
+        out: &mut Matrix<Status::Value, R3, C3, SC>,
+        dot: impl Fn(
+            &VectorView<'_, T, R1, SA::RStride, SA::CStride>,
+            &VectorView<'_, T, R2, SB::RStride, SB::CStride>,
+        ) -> T,
+    ) where
+        Status: InitStatus<T>,
+        SB: RawStorage<T, R2, C2>,
+        SC: RawStorageMut<Status::Value, R3, C3>,
+        ShapeConstraint: SameNumberOfColumns<C1, C2> + DimEq<R1, R3> + DimEq<R2, C3>,
+    {
+        let (nrows1, ncols1) = self.shape();
+        let (nrows2, ncols2) = rhs.shape();
+        let (nrows3, ncols3) = out.shape();
+
+        assert!(
+            ncols1 == ncols2,
+            "Matrix multiplication dimensions mismatch {:?} and {:?}: left rows != right rows.",
+            self.shape(),
+            rhs.shape()
+        );
+        assert!(
+            nrows1 == nrows3,
+            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right rows.",
+            self.shape(),
+            out.shape()
+        );
+        assert!(
+            nrows2 == ncols3,
+            "Matrix multiplication output dimensions mismatch {:?} and {:?}: left cols != right cols",
+            rhs.shape(),
+            out.shape()
+        );
+
+        for i in 0..nrows1 {
+            for j in 0..nrows2 {
+                let dot = dot(&self.row(i), &rhs.row(j)); 
+                let elt = unsafe { out.get_unchecked_mut((i, j)) };
+                Status::init(elt, dot)
+            }
+        }
+    }
+
     /// Equivalent to `self.transpose() * rhs` but stores the result into `out` to avoid
     /// allocations.
     #[inline]

tests/core/matrix.rs

@@ -920,9 +920,13 @@ mod transposition_tests {
         }
 
         #[test]
-        fn tr_mul_is_transpose_then_mul(m in matrix(PROPTEST_F64, Const::<4>, Const::<6>), v in vector4()) {
+        fn tr_mul_is_transpose_lhs_then_mul(m in matrix(PROPTEST_F64, Const::<4>, Const::<6>), v in vector4()) {
             prop_assert!(relative_eq!(m.transpose() * v, m.tr_mul(&v), epsilon = 1.0e-7))
         }
+        #[test]
+        fn mul_tr_is_transpose_rhs_then_mul(m in matrix(PROPTEST_F64, Const::<4>, Const::<6>), v in vector4()) {
+            prop_assert!(relative_eq!(m*v.transpose(), m.mul_tr(&v), epsilon = 1.0e-7))
+        }
     }
 }