refactor: removed vector and only use Row/Col Matrices. Added a rebracketing module. STill uncertain on how to use it. TODO: the computing the matrices product sould take into acount if the matrix are unit matrices and the types used

Author: Alexandre Hoffmann

include/FSLinalg/BasicLinalg.hpp

@@ -10,10 +10,14 @@
 #include <FSLinalg/Matrix/VectorCross.hpp>
 #include <FSLinalg/Matrix/StripSymbolsAndEvalMatrix.hpp>
 #include <FSLinalg/Matrix/UnitMatrix.hpp>
+#include <FSLinalg/Matrix/MatrixProductAnalyzer.hpp>
+#include <FSLinalg/Matrix/MatrixProductChain.hpp>
 
 #include <FSLinalg/Matrix/MatrixBase_impl.hpp>
 #include <FSLinalg/Matrix/MatrixConj_impl.hpp>
 #include <FSLinalg/Matrix/MatrixProduct_impl.hpp>
 #include <FSLinalg/Matrix/Matrix_impl.hpp>
 #include <FSLinalg/Matrix/MatrixTransposed_impl.hpp>
 #include <FSLinalg/Matrix/VectorCross_impl.hpp>
+#include <FSLinalg/Matrix/MatrixProductAnalyzer_impl.hpp>
+#include <FSLinalg/Matrix/MatrixProductChain_impl.hpp>

include/FSLinalg/Matrix.hpp

@@ -33,6 +33,13 @@ class Matrix : public MatrixBase< Matrix<T, Nrows, Ncols> >
 	using Self = Matrix<T, Nrows, Ncols>;
 	FSLINALG_DEFINE_MATRIX
 
+	template<class Dst>
+	struct CanBeAlisaedTo : std::bool_constant< 
+		    IsMatrix<Dst>::value 
+		and Base::nRows == Dst::nRows
+		and Base::nCols == Dst::nCols
+		and std::is_same<Scalar, typename Dst::Scalar>::value > {};
+	
 	static constexpr bool isScalarComplex = IsComplexScalar<Scalar>::value;
 	static constexpr bool isVector = isRowVector or isColVector;
 
@@ -46,11 +53,11 @@ class Matrix : public MatrixBase< Matrix<T, Nrows, Ncols> >
 
 	Matrix(const Matrix& other) : m_data(other.m_data) {}
 
-	template<class Expr> Matrix(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.assignTo(1., *this, std::false_type{}); }
+	template<class Expr> Matrix(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.assignTo(Scalar(1), *this, std::false_type{}); }
 
-	template<class Expr> Matrix& operator= (const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.assignTo  (1., *this, std::true_type{}); return *this; }
-	template<class Expr> Matrix& operator+=(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.increment (1., *this, std::true_type{}); return *this; }
-	template<class Expr> Matrix& operator-=(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.decrement (1., *this, std::true_type{}); return *this; }
+	template<class Expr> Matrix& operator= (const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.assignTo  (Scalar(1), *this, std::true_type{}); return *this; }
+	template<class Expr> Matrix& operator+=(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.increment (Scalar(1), *this, std::true_type{}); return *this; }
+	template<class Expr> Matrix& operator-=(const MatrixBase<Expr>& expr) requires(IsConstructibleFrom<Expr>::value) { expr.decrement (Scalar(1), *this, std::true_type{}); return *this; }
 
 	Matrix& operator*=(const RealScalar& alpha) requires(isScalarComplex) { for (Size i=0; i!=size; ++i) { m_data[i] *= alpha; } return *this; }
 	Matrix& operator/=(const RealScalar& alpha) requires(isScalarComplex) { for (Size i=0; i!=size; ++i) { m_data[i] /= alpha; } return *this; }
@@ -64,8 +71,8 @@ class Matrix : public MatrixBase< Matrix<T, Nrows, Ncols> >
 	const_ReturnType getImpl(const Size i, const Size j) const { return m_data[i*nCols + j]; }
 	      ReturnType getImpl(const Size i, const Size j)       { return m_data[i*nCols + j]; }
 
-	template<class Dst>           bool isAliasedToImpl(const MatrixBase<Dst>& dst) const requires(nRows == Dst::nRows and nCols == Dst::nCols) { return std::addressof(dst.derived()) == this; }
-	template<class Dst> constexpr bool isAliasedToImpl(const MatrixBase<Dst>&    ) const requires(nRows != Dst::nRows or  nCols != Dst::nCols) { return false; }
+	template<class Dst>           bool isAliasedToImpl(const MatrixBase<Dst>& dst) const requires(    CanBeAlisaedTo<Dst>::value) { return std::addressof(dst.derived()) == this; }
+	template<class Dst> constexpr bool isAliasedToImpl(const MatrixBase<Dst>&    ) const requires(not CanBeAlisaedTo<Dst>::value) { return false; }
 private:
 	std::array<Scalar, size> m_data;
 };

include/FSLinalg/Matrix/Matrix.hpp

@@ -25,13 +25,13 @@ class MatrixBase : public CRTPBase<Derived>
 
 	template<class Dst> 
 	struct IsConvertibleTo : std::bool_constant<
-		std::is_base_of<MatrixBase<Dst>, Dst>::value 
+		    std::is_base_of<MatrixBase<Dst>, Dst>::value 
 		and Dst::hasWriteRandomAccess
 		and std::is_convertible<Scalar, typename Dst::Scalar>::value> {};
 
 	template<class Src> 
 	struct IsConstructibleFrom : std::bool_constant<
-		std::is_base_of<MatrixBase<Src>, Src>::value 
+		    std::is_base_of<MatrixBase<Src>, Src>::value 
 		and DerivedTraits::hasWriteRandomAccess
 		and std::is_convertible<typename Src::Scalar, Scalar>::value> {};
 

include/FSLinalg/Matrix/MatrixBase.hpp

@@ -1,13 +1,21 @@
-#ifndef FSLINALG_MATRIX_MATRIX_PRODUCT_HPP
-#define FSLINALG_MATRIX_MATRIX_PRODUCT_HPP
+#ifndef FSLINALG_MATRIX_PRODUCT_HPP
+#define FSLINALG_MATRIX_PRODUCT_HPP
 
 #include <FSLinalg/Matrix/MatrixBase.hpp>
 #include <FSLinalg/Matrix/MatrixTransposed.hpp>
 #include <FSLinalg/Matrix/StripSymbolsAndEvalMatrix.hpp>
+#include <FSLinalg/Matrix/MatrixProductAnalyzer.hpp>
 
 namespace FSLinalg
 {
 
+namespace detail
+{
+
+template<class Expr> struct MatrixProductAnalyzerImpl;
+
+} // namespace detail
+
 template<class Lhs, class Rhs> class MatrixProduct;
 
 template<class Lhs, class Rhs>
@@ -19,7 +27,7 @@ struct MatrixTraits< MatrixProduct<Lhs, Rhs> >
 	using Scalar = decltype(std::declval<typename Lhs::Scalar>() * std::declval<typename Rhs::Scalar>());
 	using Size   = std::common_type_t<typename Lhs::Size, typename Rhs::Size>;
 
-	static constexpr bool hasReadRandomAccess  = Lhs::isRowVector and Rhs::isColVector and Lhs::hasFlatRandomAccess and Rhs::hasFlatRandomAccess;
+	static constexpr bool hasReadRandomAccess  = (Lhs::isRowVector and Lhs::hasFlatRandomAccess) or (Rhs::isColVector and Rhs::hasFlatRandomAccess);
 	static constexpr bool hasWriteRandomAccess = false;
 	static constexpr bool hasFlatRandomAccess  = false;
 	static constexpr bool causesAliasingIssues = true;
@@ -36,10 +44,15 @@ class MatrixProduct : public MatrixBase< MatrixProduct<Lhs,Rhs> >
 	using Self = MatrixProduct<Lhs,Rhs>;
 	FSLINALG_DEFINE_MATRIX
 
+	using OptimallyBracketedSelf = MatrixProductAnalyzer<Self>;
+	
+	friend struct detail::MatrixProductAnalyzerImpl< Self >;
+	
 	MatrixProduct(const MatrixBase<Lhs>& lhs, const MatrixBase<Rhs>& rhs) : m_lhs(lhs.derived()), m_rhs(rhs.derived()) {}
 
 	static constexpr bool createTemporaryLhs = StripSymbolsAndEvalMatrix<Lhs>::createsTemporary;
 	static constexpr bool createTemporaryRhs = StripSymbolsAndEvalMatrix<Rhs>::createsTemporary;
+	static constexpr bool isOptimallyBracked = std::is_same<Self, OptimallyBracketedSelf>::value;
 
 	/**
 	 * @brief Only awailable when multiplying row-vector and col-vector
@@ -86,4 +99,4 @@ MatrixProduct< Lhs,MatrixTransposed<Rhs> > outer(const MatrixBase<Lhs>& lhs, con
 
 } // namespace FSLinalg
 
-#endif // FSLINALG_MATRIX_MATRIX_PRODUCT_HPP
+#endif // FSLINALG_MATRIX_PRODUCT_HPP

include/FSLinalg/Matrix/MatrixProduct.hpp

@@ -0,0 +1,125 @@
+#ifndef FSLINALG_MATRIX_PRODUCT_TRAITS_HPP
+#define FSLINALG_MATRIX_PRODUCT_TRAITS_HPP
+
+#include <cstddef>
+#include <type_traits>
+
+#include <FSLinalg/Matrix/MatrixProductChain.hpp>
+
+namespace FSLinalg
+{
+
+template<class Lhs, class Rhs> class MatrixProduct;
+	
+namespace detail
+{
+
+template<class Expr>
+struct MatrixProductAnalyzerImpl
+{
+	static constexpr size_t length = 1;
+	
+	template<size_t n> requires(n == 0)
+	using NthMatrix = Expr;
+	
+	template<size_t n>
+	static constexpr const NthMatrix<n>& getMatrix(const Expr& expr, std::integral_constant<size_t, n>) { return expr; }
+	
+	template<size_t idx, size_t chainLenP1> requires (idx+1 <chainLenP1) 
+	static constexpr void fillDims(std::integral_constant<size_t, idx>, std::array<size_t, chainLenP1>& dims) { dims[idx] = Expr::nRows; dims[idx+1] = Expr::nCols; }
+};
+
+template<class Lhs, class Rhs>
+struct MatrixProductAnalyzerImpl< MatrixProduct<Lhs, Rhs> >
+{
+	static constexpr size_t lhsLength = MatrixProductAnalyzerImpl<Lhs>::length;
+	static constexpr size_t rhsLength = MatrixProductAnalyzerImpl<Rhs>::length;
+	static constexpr size_t length    = lhsLength + rhsLength;
+	
+	template<size_t n, class Enable = void>
+	struct NthMatrixHelper;
+	
+	template<size_t n> 
+	struct NthMatrixHelper<n, std::enable_if_t<n < lhsLength>>
+	{
+		using Type = typename MatrixProductAnalyzerImpl<Lhs>::template NthMatrix<n>;
+	};
+	
+	template<size_t n> 
+	struct NthMatrixHelper<n, std::enable_if_t<lhsLength <= n and n < length>>
+	{
+		using Type = typename MatrixProductAnalyzerImpl<Rhs>::template NthMatrix<n - lhsLength>;
+	};
+
+	template<size_t n> 
+	using NthMatrix = typename NthMatrixHelper<n>::Type;
+	
+	template<size_t n>
+	static constexpr const NthMatrix<n>& getMatrix(const MatrixProduct<Lhs, Rhs>& expr, std::integral_constant<size_t, n>);
+
+	template<size_t idx, size_t chainLenP1> requires (idx+1 <chainLenP1) 
+	static constexpr void fillDims(std::integral_constant<size_t, idx>, std::array<size_t, chainLenP1>& dims);
+};
+
+} // namespace detail
+
+template<class Expr>
+struct MatrixProductAnalyzer
+{
+	static_assert(IsMatrix<Expr>::value, "Expr must be a matrix");
+	
+	using Impl = detail::MatrixProductAnalyzerImpl<Expr>;
+	using DimArray = std::array<size_t, Impl::length+1>;
+	
+	template<size_t n> using NthMatrix = typename Impl::template NthMatrix<n>;
+	
+	template<size_t n>
+	static const NthMatrix<n>& getMatrix(const Expr& expr, std::integral_constant<size_t, n> ic) { return Impl::getMatrix(expr, ic); }
+	
+	static constexpr size_t getLength() { return Impl::length; }
+	
+	static constexpr DimArray getDims();
+	
+	// we use an external template class as to not recompute everything for every product
+	// if the optimal splitting for an chain with the same dims has already been computed 
+	// we can re-use it.
+	static constexpr size_t getMinalCost()    { return MatrixProductChain<getDims()>::template minCostAndSplit<0, getLength()>.first;  }
+	static constexpr size_t getOptimalSplit() { return MatrixProductChain<getDims()>::template minCostAndSplit<0, getLength()>.second; }
+private:
+	template<size_t start, size_t end> requires(start <= end and end <= getLength())
+	struct OptimalBracketingHelper
+	{
+		static constexpr size_t split = MatrixProductChain<getDims()>::template minCostAndSplit<start, end>.second;
+		
+		static_assert(start <= split and split+1 < end+1, "invalid split");
+		
+		using LhsBracketing = OptimalBracketingHelper<start, split>;
+		using RhsBracketing = OptimalBracketingHelper<split, end>;
+		
+		using Lhs = typename LhsBracketing::Type;
+		using Rhs = typename RhsBracketing::Type;
+		
+		using Type          = MatrixProduct<Lhs, Rhs>;
+		using ReBracketType = Type; 
+		
+		static ReBracketType reBracket(const Expr& expr) { return ReBracketType(LhsBracketing::reBracket(expr), RhsBracketing::reBracket(expr)); }
+	};
+	
+	template<size_t idx> requires(idx < getLength())
+	struct OptimalBracketingHelper<idx,idx+1>
+	{
+		using Type          = NthMatrix<idx>;
+		using ReBracketType = const Type&;
+		
+		static ReBracketType reBracket(const Expr& expr) { return MatrixProductAnalyzer<Expr>::getMatrix(expr, std::integral_constant<size_t, idx>{}); }
+	};
+public:
+	using OptimalBracketing = typename OptimalBracketingHelper<0, getLength()>::Type;
+	using ReBracketType     = typename OptimalBracketingHelper<0, getLength()>::ReBracketType;
+	
+	static ReBracketType reBracket(const Expr& expr) { return OptimalBracketingHelper<0, getLength()>::reBracket(expr); }
+};
+
+} // namespace FSLinalg
+
+#endif // FSLINALG_MATRIX_PRODUCT_TRAITS_HPP

include/FSLinalg/Matrix/MatrixProductAnalyzer.hpp

@@ -0,0 +1,62 @@
+#ifndef FSLINALG_MATRIX_PRODUCT_TRAITS_IMPL_HPP
+#define FSLINALG_MATRIX_PRODUCT_TRAITS_IMPL_HPP
+
+#include <FSLinalg/Matrix/MatrixProductAnalyzer.hpp>
+
+namespace FSLinalg
+{
+namespace detail
+{
+
+template<class Lhs, class Rhs> template<size_t idx, size_t chainLenP1> requires (idx+1 <chainLenP1)
+constexpr void MatrixProductAnalyzerImpl< MatrixProduct<Lhs, Rhs> >::fillDims(std::integral_constant<size_t, idx> ic, std::array<size_t, chainLenP1>& dims)
+{
+	if constexpr (IsMatrixProduct<Lhs>::value and IsMatrixProduct<Rhs>::value)
+	{		
+		MatrixProductAnalyzerImpl<Lhs>::fillDims(ic, dims);
+		MatrixProductAnalyzerImpl<Rhs>::fillDims(std::integral_constant<size_t, idx + lhsLength>{}, dims);
+	}
+	else if constexpr (IsMatrixProduct<Lhs>::value)
+	{	
+		MatrixProductAnalyzerImpl<Lhs>::fillDims(ic, dims);	
+		dims[idx + lhsLength + 1] = Rhs::nCols;
+	}
+	else if constexpr (IsMatrixProduct<Rhs>::value)
+	{		
+		dims[idx] = Lhs::nRows;
+		MatrixProductAnalyzerImpl<Rhs>::fillDims(std::integral_constant<size_t, idx + 1>{}, dims);
+	}
+	else
+	{		
+		dims[idx]     = Lhs::nRows;
+		dims[idx + 1] = Lhs::nCols;
+		dims[idx + 2] = Rhs::nCols;	
+	}
+}
+
+template<class Lhs, class Rhs> template<size_t n>
+constexpr auto MatrixProductAnalyzerImpl< MatrixProduct<Lhs, Rhs> >::getMatrix(const MatrixProduct<Lhs, Rhs>& expr, std::integral_constant<size_t, n> ic) -> const NthMatrix<n>&
+{
+	if constexpr (n < lhsLength)
+	{
+		return MatrixProductAnalyzerImpl<Lhs>::getMatrix(expr.m_lhs, ic);
+	}
+	else
+	{
+		return MatrixProductAnalyzerImpl<Rhs>::getMatrix(expr.m_rhs, std::integral_constant<size_t, n - lhsLength>{});
+	}
+}
+
+} // namespace detail
+
+template<class Expr>
+constexpr auto MatrixProductAnalyzer<Expr>::getDims() -> DimArray
+{ 
+	DimArray dims; 
+	Impl::fillDims(std::integral_constant<size_t, 0>{}, dims); 
+	return dims; 
+}
+
+} // namespace FSLinalg
+
+#endif // FSLINALG_MATRIX_PRODUCT_TRAITS_IMPL_HPP

include/FSLinalg/Matrix/MatrixProductAnalyzer_impl.hpp

@@ -0,0 +1,21 @@
+#ifndef FSLINALG_MATRIX_PRODUCT_CHAIN_HPP
+#define FSLINALG_MATRIX_PRODUCT_CHAIN_HPP
+
+#include <array>
+
+namespace FSLinalg
+{
+
+template<std::array dims> 
+struct MatrixProductChain
+{
+	template<size_t i, size_t j>
+	static constexpr std::pair<size_t, size_t> minMulCostAndSplitRec(std::integral_constant<size_t, i>, std::integral_constant<size_t, j>);
+	
+	template<size_t i, size_t j> 
+	static constexpr std::pair<size_t, size_t> minCostAndSplit = minMulCostAndSplitRec(std::integral_constant<size_t, i>{}, std::integral_constant<size_t, j>{});
+};
+
+} // namespace FSLinalg
+
+#endif // FSLINALG_MATRIX_PRODUCT_CHAIN_HPP

include/FSLinalg/Matrix/MatrixProductChain.hpp

@@ -0,0 +1,38 @@
+#ifndef FSLINALG_MATRIX_PRODUCT_CHAIN_IMPL_HPP
+#define FSLINALG_MATRIX_PRODUCT_CHAIN_IMPL_HPP
+
+#include <FSLinalg/Matrix/MatrixProductChain.hpp>
+#include <FSLinalg/misc/StaticFor.hpp>
+
+namespace FSLinalg
+{
+
+template<std::array dims> template<size_t i, size_t j>
+constexpr std::pair<size_t, size_t> MatrixProductChain<dims>::minMulCostAndSplitRec(std::integral_constant<size_t, i>, std::integral_constant<size_t, j>)
+{
+	if constexpr (i + 1 == j) 
+	{ 
+		return std::make_pair(0u, i + 1);
+	}
+	else
+	{
+		size_t minCost     = std::numeric_limits<size_t>::max();
+		size_t optSpliting = i + 1; 
+		
+		misc::StaticFor<i + 1, j>::run([&minCost, &optSpliting](const auto k) -> void
+		{
+			constexpr size_t curr = minCostAndSplit<i, k>.first + minCostAndSplit<k, j>.first + dims[i]*dims[k]*dims[j];
+			if (curr < minCost)
+			{
+				minCost = curr;
+				optSpliting = k;
+			}
+		});
+		
+		return std::make_pair(minCost, optSpliting);
+	}
+}
+
+} //namespace FSLinalg
+
+#endif // FSLINALG_MATRIX_PRODUCT_CHAIN_IMPL_HPP

include/FSLinalg/Matrix/MatrixProductChain_impl.hpp

@@ -1,5 +1,5 @@
-#ifndef FSLINALG_MATRIX_MATRIX_PRODUCT_IMPL_HPP
-#define FSLINALG_MATRIX_MATRIX_PRODUCT_IMPL_HPP
+#ifndef FSLINALG_MATRIX_PRODUCT_IMPL_HPP
+#define FSLINALG_MATRIX_PRODUCT_IMPL_HPP
 
 #include <FSLinalg/Matrix/MatrixProduct.hpp>
 #include <FSLinalg/BasicLinalg/GeneralMatrixMatrixProduct.hpp>
@@ -86,4 +86,4 @@ void MatrixProduct<Lhs,Rhs>::decrementImpl(const Alpha& alpha, MatrixBase<Dst>&
 
 } // namespace FSLinalg
 
-#endif // FSLINALG_MATRIX_MATRIX_PRODUCT_IMPL_HPP
+#endif // FSLINALG_MATRIX_PRODUCT_IMPL_HPP