1 //=================================================================================================
2 /*!
3 //  \file blaze/math/expressions/TDMatSMatMultExpr.h
4 //  \brief Header file for the transpose dense matrix/sparse matrix multiplication expression
5 //
6 //  Copyright (C) 2012-2020 Klaus Iglberger - All Rights Reserved
7 //
8 //  This file is part of the Blaze library. You can redistribute it and/or modify it under
9 //  the terms of the New (Revised) BSD License. Redistribution and use in source and binary
10 //  forms, with or without modification, are permitted provided that the following conditions
11 //  are met:
12 //
13 //  1. Redistributions of source code must retain the above copyright notice, this list of
14 //     conditions and the following disclaimer.
15 //  2. Redistributions in binary form must reproduce the above copyright notice, this list
16 //     of conditions and the following disclaimer in the documentation and/or other materials
17 //     provided with the distribution.
18 //  3. Neither the names of the Blaze development group nor the names of its contributors
19 //     may be used to endorse or promote products derived from this software without specific
20 //     prior written permission.
21 //
22 //  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
23 //  EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
24 //  OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
25 //  SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 //  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 //  TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
28 //  BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 //  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 //  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
31 //  DAMAGE.
32 */
33 //=================================================================================================
34 
35 #ifndef _BLAZE_MATH_EXPRESSIONS_TDMATSMATMULTEXPR_H_
36 #define _BLAZE_MATH_EXPRESSIONS_TDMATSMATMULTEXPR_H_
37 
38 
39 //*************************************************************************************************
40 // Includes
41 //*************************************************************************************************
42 
43 #include <blaze/math/Aliases.h>
44 #include <blaze/math/constraints/ColumnMajorMatrix.h>
45 #include <blaze/math/constraints/DenseMatrix.h>
46 #include <blaze/math/constraints/MatMatMultExpr.h>
47 #include <blaze/math/constraints/RequiresEvaluation.h>
48 #include <blaze/math/constraints/RowMajorMatrix.h>
49 #include <blaze/math/constraints/SparseMatrix.h>
50 #include <blaze/math/constraints/StorageOrder.h>
51 #include <blaze/math/constraints/Symmetric.h>
52 #include <blaze/math/constraints/Zero.h>
53 #include <blaze/math/Exception.h>
54 #include <blaze/math/expressions/Computation.h>
55 #include <blaze/math/expressions/DenseMatrix.h>
56 #include <blaze/math/expressions/Forward.h>
57 #include <blaze/math/expressions/MatMatMultExpr.h>
58 #include <blaze/math/functors/DeclDiag.h>
59 #include <blaze/math/functors/DeclHerm.h>
60 #include <blaze/math/functors/DeclLow.h>
61 #include <blaze/math/functors/DeclSym.h>
62 #include <blaze/math/functors/DeclUpp.h>
63 #include <blaze/math/functors/Noop.h>
64 #include <blaze/math/shims/Conjugate.h>
65 #include <blaze/math/shims/IsDefault.h>
66 #include <blaze/math/shims/PrevMultiple.h>
67 #include <blaze/math/shims/Reset.h>
68 #include <blaze/math/shims/Serial.h>
69 #include <blaze/math/traits/DeclDiagTrait.h>
70 #include <blaze/math/traits/DeclHermTrait.h>
71 #include <blaze/math/traits/DeclLowTrait.h>
72 #include <blaze/math/traits/DeclSymTrait.h>
73 #include <blaze/math/traits/DeclUppTrait.h>
74 #include <blaze/math/traits/MultTrait.h>
75 #include <blaze/math/typetraits/IsAligned.h>
76 #include <blaze/math/typetraits/IsColumnMajorMatrix.h>
77 #include <blaze/math/typetraits/IsComputation.h>
78 #include <blaze/math/typetraits/IsDiagonal.h>
79 #include <blaze/math/typetraits/IsExpression.h>
80 #include <blaze/math/typetraits/IsIdentity.h>
81 #include <blaze/math/typetraits/IsLower.h>
82 #include <blaze/math/typetraits/IsResizable.h>
83 #include <blaze/math/typetraits/IsStrictlyLower.h>
84 #include <blaze/math/typetraits/IsStrictlyUpper.h>
85 #include <blaze/math/typetraits/IsSymmetric.h>
86 #include <blaze/math/typetraits/IsTriangular.h>
87 #include <blaze/math/typetraits/IsUpper.h>
88 #include <blaze/math/typetraits/IsZero.h>
89 #include <blaze/math/typetraits/RequiresEvaluation.h>
90 #include <blaze/math/typetraits/Size.h>
91 #include <blaze/math/views/Check.h>
92 #include <blaze/system/Optimizations.h>
93 #include <blaze/system/Thresholds.h>
94 #include <blaze/util/algorithms/Max.h>
95 #include <blaze/util/algorithms/Min.h>
96 #include <blaze/util/Assert.h>
97 #include <blaze/util/EnableIf.h>
98 #include <blaze/util/FunctionTrace.h>
99 #include <blaze/util/IntegralConstant.h>
100 #include <blaze/util/MaybeUnused.h>
101 #include <blaze/util/mpl/If.h>
102 #include <blaze/util/Types.h>
103 #include <blaze/util/typetraits/IsBuiltin.h>
104 
105 
106 namespace blaze {
107 
108 //=================================================================================================
109 //
110 //  CLASS TDMATSMATMULTEXPR
111 //
112 //=================================================================================================
113 
114 //*************************************************************************************************
115 /*!\brief Expression object for transpose dense matrix-sparse matrix multiplications.
116 // \ingroup dense_matrix_expression
117 //
118 // The TDMatSMatMultExpr class represents the compile time expression for multiplications between
119 // a column-major dense matrix and a row-major sparse matrix.
120 */
121 template< typename MT1  // Type of the left-hand side dense matrix
122         , typename MT2  // Type of the right-hand side sparse matrix
123         , bool SF       // Symmetry flag
124         , bool HF       // Hermitian flag
125         , bool LF       // Lower flag
126         , bool UF >     // Upper flag
127 class TDMatSMatMultExpr
128    : public MatMatMultExpr< DenseMatrix< TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>, true > >
129    , private Computation
130 {
131  private:
132    //**Type definitions****************************************************************************
133    using RT1 = ResultType_t<MT1>;     //!< Result type of the left-hand side dense matrix expression.
134    using RT2 = ResultType_t<MT2>;     //!< Result type of the right-hand side sparse matrix expression.
135    using ET1 = ElementType_t<RT1>;    //!< Element type of the left-hand side dense matrix expression.
136    using ET2 = ElementType_t<RT2>;    //!< Element type of the right-hand side sparse matrix expression.
137    using CT1 = CompositeType_t<MT1>;  //!< Composite type of the left-hand side dense matrix expression.
138    using CT2 = CompositeType_t<MT2>;  //!< Composite type of the right-hand side sparse matrix expression.
139    //**********************************************************************************************
140 
141    //**********************************************************************************************
142    //! Compilation switch for the composite type of the left-hand side dense matrix expression.
143    static constexpr bool evaluateLeft = ( IsComputation_v<MT1> || RequiresEvaluation_v<MT1> );
144    //**********************************************************************************************
145 
146    //**********************************************************************************************
147    //! Compilation switch for the composite type of the right-hand side sparse matrix expression.
148    static constexpr bool evaluateRight = ( IsComputation_v<MT2> || RequiresEvaluation_v<MT2> );
149    //**********************************************************************************************
150 
151    //**********************************************************************************************
152    static constexpr bool SYM  = ( SF && !( HF || LF || UF )    );  //!< Flag for symmetric matrices.
153    static constexpr bool HERM = ( HF && !( LF || UF )          );  //!< Flag for Hermitian matrices.
154    static constexpr bool LOW  = ( LF || ( ( SF || HF ) && UF ) );  //!< Flag for lower matrices.
155    static constexpr bool UPP  = ( UF || ( ( SF || HF ) && LF ) );  //!< Flag for upper matrices.
156    //**********************************************************************************************
157 
158    //**********************************************************************************************
159    /*! \cond BLAZE_INTERNAL */
160    //! Helper variable template for the explicit application of the SFINAE principle.
161    /*! This variable template is a helper for the selection of the optimal evaluation strategy.
162        In case the right-hand side matrix operands of type \a T3 is symmetric, the variable is
163        set to 1 and an optimized evaluation strategy is selected. Otherwise the variable is set
164        to 0 and the default strategy is chosen. */
165    template< typename T1, typename T2, typename T3 >
166    static constexpr bool CanExploitSymmetry_v = IsSymmetric_v<T3>;
167    /*! \endcond */
168    //**********************************************************************************************
169 
170    //**********************************************************************************************
171    /*! \cond BLAZE_INTERNAL */
172    //! Helper variable template for the explicit application of the SFINAE principle.
173    /*! This variable template is a helper for the selection of the parallel evaluation strategy.
174        In case either of the two matrix operands requires an intermediate evaluation, the variable
175        will be set to 1, otherwise it will be 0. */
176    template< typename T1, typename T2, typename T3 >
177    static constexpr bool IsEvaluationRequired_v =
178       ( ( evaluateLeft || evaluateRight ) && !CanExploitSymmetry_v<T1,T2,T3> );
179    /*! \endcond */
180    //**********************************************************************************************
181 
182    //**********************************************************************************************
183    /*! \cond BLAZE_INTERNAL */
184    //! Helper variable template for the explicit application of the SFINAE principle.
185    /*! In case no SMP assignment is required and the element type of the target matrix has a
186        fixed size (i.e. is not resizable), the variable will be set to 1, otherwise it will be 0. */
187    template< typename T1, typename T2, typename T3 >
188    static constexpr bool UseOptimizedKernel_v =
189       ( useOptimizedKernels &&
190         !IsDiagonal_v<T2> &&
191         !IsResizable_v< ElementType_t<T1> > &&
192         !( IsColumnMajorMatrix_v<T1> && IsResizable_v<ET2> ) );
193    /*! \endcond */
194    //**********************************************************************************************
195 
196    //**********************************************************************************************
197    /*! \cond BLAZE_INTERNAL */
198    //! Helper variable template for the explicit application of the SFINAE principle.
199    /*! In case no SMP assignment is required and the element type of the target matrix is
200        resizable, the variable will be set to 1, otherwise it will be 0. */
201    template< typename T1, typename T2, typename T3 >
202    static constexpr bool UseDefaultKernel_v = !UseOptimizedKernel_v<T1,T2,T3>;
203    /*! \endcond */
204    //**********************************************************************************************
205 
206    //**********************************************************************************************
207    /*! \cond BLAZE_INTERNAL */
208    //! Type of the functor for forwarding an expression to another assign kernel.
209    /*! In case a temporary matrix needs to be created, this functor is used to forward the
210        resulting expression to another assign kernel. */
211    using ForwardFunctor = If_t< HERM
212                               , DeclHerm
213                               , If_t< SYM
214                                     , DeclSym
215                                     , If_t< LOW
216                                           , If_t< UPP
217                                                 , DeclDiag
218                                                 , DeclLow >
219                                           , If_t< UPP
220                                                 , DeclUpp
221                                                 , Noop > > > >;
222    /*! \endcond */
223    //**********************************************************************************************
224 
225  public:
226    //**Type definitions****************************************************************************
227    //! Type of this TDMatSMatMultExpr instance.
228    using This = TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>;
229 
230    //! Base type of this TDMatSMatMultExpr instance.
231    using BaseType = MatMatMultExpr< DenseMatrix<This,true> >;
232 
233    //! Result type for expression template evaluations.
234    using ResultType = typename If_t< HERM
235                                    , DeclHermTrait< MultTrait_t<RT1,RT2> >
236                                    , If_t< SYM
237                                          , DeclSymTrait< MultTrait_t<RT1,RT2> >
238                                          , If_t< LOW
239                                                , If_t< UPP
240                                                      , DeclDiagTrait< MultTrait_t<RT1,RT2> >
241                                                      , DeclLowTrait< MultTrait_t<RT1,RT2> > >
242                                                , If_t< UPP
243                                                      , DeclUppTrait< MultTrait_t<RT1,RT2> >
244                                                      , MultTrait<RT1,RT2> > > > >::Type;
245 
246    using OppositeType  = OppositeType_t<ResultType>;   //!< Result type with opposite storage order for expression template evaluations.
247    using TransposeType = TransposeType_t<ResultType>;  //!< Transpose type for expression template evaluations.
248    using ElementType   = ElementType_t<ResultType>;    //!< Resulting element type.
249    using ReturnType    = const ElementType;            //!< Return type for expression template evaluations.
250    using CompositeType = const ResultType;             //!< Data type for composite expression templates.
251 
252    //! Composite type of the left-hand side dense matrix expression.
253    using LeftOperand = If_t< IsExpression_v<MT1>, const MT1, const MT1& >;
254 
255    //! Composite type of the right-hand side sparse matrix expression.
256    using RightOperand = If_t< IsExpression_v<MT2>, const MT2, const MT2& >;
257 
258    //! Type for the assignment of the left-hand side dense matrix operand.
259    using LT = If_t< evaluateLeft, const RT1, CT1 >;
260 
261    //! Type for the assignment of the right-hand side sparse matrix operand.
262    using RT = If_t< evaluateRight, const RT2, CT2 >;
263    //**********************************************************************************************
264 
265    //**Compilation flags***************************************************************************
266    //! Compilation switch for the expression template evaluation strategy.
267    static constexpr bool simdEnabled = false;
268 
269    //! Compilation switch for the expression template assignment strategy.
270    static constexpr bool smpAssignable =
271       ( !evaluateLeft && MT1::smpAssignable && !evaluateRight && MT2::smpAssignable );
272    //**********************************************************************************************
273 
274    //**Constructor*********************************************************************************
275    /*!\brief Constructor for the TDMatSMatMultExpr class.
276    //
277    // \param lhs The left-hand side dense matrix operand of the multiplication expression.
278    // \param rhs The right-hand side sparse matrix operand of the multiplication expression.
279    */
TDMatSMatMultExpr(const MT1 & lhs,const MT2 & rhs)280    inline TDMatSMatMultExpr( const MT1& lhs, const MT2& rhs ) noexcept
281       : lhs_( lhs )  // Left-hand side dense matrix of the multiplication expression
282       , rhs_( rhs )  // Right-hand side sparse matrix of the multiplication expression
283    {
284       BLAZE_INTERNAL_ASSERT( lhs.columns() == rhs.rows(), "Invalid matrix sizes" );
285    }
286    //**********************************************************************************************
287 
288    //**Access operator*****************************************************************************
289    /*!\brief 2D-access to the matrix elements.
290    //
291    // \param i Access index for the row. The index has to be in the range \f$[0..M-1]\f$.
292    // \param j Access index for the column. The index has to be in the range \f$[0..N-1]\f$.
293    // \return The resulting value.
294    */
operator()295    inline ReturnType operator()( size_t i, size_t j ) const {
296       BLAZE_INTERNAL_ASSERT( i < lhs_.rows()   , "Invalid row access index"    );
297       BLAZE_INTERNAL_ASSERT( j < rhs_.columns(), "Invalid column access index" );
298 
299       if( IsDiagonal_v<MT1> ) {
300          return lhs_(i,i) * rhs_(i,j);
301       }
302       else if( IsDiagonal_v<MT2> ) {
303          return lhs_(i,j) * rhs_(j,j);
304       }
305       else if( IsTriangular_v<MT1> || IsTriangular_v<MT2> ) {
306          const size_t begin( ( IsUpper_v<MT1> )
307                              ?( ( IsLower_v<MT2> )
308                                 ?( max( ( IsStrictlyUpper_v<MT1> ? i+1UL : i )
309                                       , ( IsStrictlyLower_v<MT2> ? j+1UL : j ) ) )
310                                 :( IsStrictlyUpper_v<MT1> ? i+1UL : i ) )
311                              :( ( IsLower_v<MT2> )
312                                 ?( IsStrictlyLower_v<MT2> ? j+1UL : j )
313                                 :( 0UL ) ) );
314          const size_t end( ( IsLower_v<MT1> )
315                            ?( ( IsUpper_v<MT2> )
316                               ?( min( ( IsStrictlyLower_v<MT1> ? i : i+1UL )
317                                     , ( IsStrictlyUpper_v<MT2> ? j : j+1UL ) ) )
318                               :( IsStrictlyLower_v<MT1> ? i : i+1UL ) )
319                            :( ( IsUpper_v<MT2> )
320                               ?( IsStrictlyUpper_v<MT2> ? j : j+1UL )
321                               :( lhs_.columns() ) ) );
322 
323          if( begin >= end ) return ElementType();
324 
325          const size_t n( end - begin );
326 
327          return subvector( row( lhs_, i, unchecked ), begin, n, unchecked ) *
328                 subvector( column( rhs_, j, unchecked ), begin, n, unchecked );
329       }
330       else {
331          return row( lhs_, i, unchecked ) * column( rhs_, j, unchecked );
332       }
333    }
334    //**********************************************************************************************
335 
336    //**At function*********************************************************************************
337    /*!\brief Checked access to the matrix elements.
338    //
339    // \param i Access index for the row. The index has to be in the range \f$[0..M-1]\f$.
340    // \param j Access index for the column. The index has to be in the range \f$[0..N-1]\f$.
341    // \return The resulting value.
342    // \exception std::out_of_range Invalid matrix access index.
343    */
at(size_t i,size_t j)344    inline ReturnType at( size_t i, size_t j ) const {
345       if( i >= lhs_.rows() ) {
346          BLAZE_THROW_OUT_OF_RANGE( "Invalid row access index" );
347       }
348       if( j >= rhs_.columns() ) {
349          BLAZE_THROW_OUT_OF_RANGE( "Invalid column access index" );
350       }
351       return (*this)(i,j);
352    }
353    //**********************************************************************************************
354 
355    //**Rows function*******************************************************************************
356    /*!\brief Returns the current number of rows of the matrix.
357    //
358    // \return The number of rows of the matrix.
359    */
rows()360    inline size_t rows() const noexcept {
361       return lhs_.rows();
362    }
363    //**********************************************************************************************
364 
365    //**Columns function****************************************************************************
366    /*!\brief Returns the current number of columns of the matrix.
367    //
368    // \return The number of columns of the matrix.
369    */
columns()370    inline size_t columns() const noexcept {
371       return rhs_.columns();
372    }
373    //**********************************************************************************************
374 
375    //**Left operand access*************************************************************************
376    /*!\brief Returns the left-hand side transpose dense matrix operand.
377    //
378    // \return The left-hand side transpose dense matrix operand.
379    */
leftOperand()380    inline LeftOperand leftOperand() const noexcept {
381       return lhs_;
382    }
383    //**********************************************************************************************
384 
385    //**Right operand access************************************************************************
386    /*!\brief Returns the right-hand side sparse matrix operand.
387    //
388    // \return The right-hand side sparse matrix operand.
389    */
rightOperand()390    inline RightOperand rightOperand() const noexcept {
391       return rhs_;
392    }
393    //**********************************************************************************************
394 
395    //**********************************************************************************************
396    /*!\brief Returns whether the expression can alias with the given address \a alias.
397    //
398    // \param alias The alias to be checked.
399    // \return \a true in case the expression can alias, \a false otherwise.
400    */
401    template< typename T >
canAlias(const T * alias)402    inline bool canAlias( const T* alias ) const noexcept {
403       return ( lhs_.isAliased( alias ) || rhs_.isAliased( alias ) );
404    }
405    //**********************************************************************************************
406 
407    //**********************************************************************************************
408    /*!\brief Returns whether the expression is aliased with the given address \a alias.
409    //
410    // \param alias The alias to be checked.
411    // \return \a true in case an alias effect is detected, \a false otherwise.
412    */
413    template< typename T >
isAliased(const T * alias)414    inline bool isAliased( const T* alias ) const {
415       return ( lhs_.isAliased( alias ) || rhs_.isAliased( alias ) );
416    }
417    //**********************************************************************************************
418 
419    //**********************************************************************************************
420    /*!\brief Returns whether the operands of the expression are properly aligned in memory.
421    //
422    // \return \a true in case the operands are aligned, \a false if not.
423    */
isAligned()424    inline bool isAligned() const noexcept {
425       return lhs_.isAligned();
426    }
427    //**********************************************************************************************
428 
429    //**********************************************************************************************
430    /*!\brief Returns whether the expression can be used in SMP assignments.
431    //
432    // \return \a true in case the expression can be used in SMP assignments, \a false if not.
433    */
canSMPAssign()434    inline bool canSMPAssign() const noexcept {
435       return ( rows() * columns() >= SMP_TDMATSMATMULT_THRESHOLD ) && !IsDiagonal_v<MT1>;
436    }
437    //**********************************************************************************************
438 
439  private:
440    //**Member variables****************************************************************************
441    LeftOperand  lhs_;  //!< Left-hand side dense matrix of the multiplication expression.
442    RightOperand rhs_;  //!< Right-hand side sparse matrix of the multiplication expression.
443    //**********************************************************************************************
444 
445    //**Assignment to dense matrices****************************************************************
446    /*! \cond BLAZE_INTERNAL */
447    /*!\brief Assignment of a transpose dense matrix-sparse matrix multiplication to a dense matrix
448    //        (\f$ C=A*B \f$).
449    // \ingroup dense_matrix
450    //
451    // \param lhs The target left-hand side dense matrix.
452    // \param rhs The right-hand side multiplication expression to be assigned.
453    // \return void
454    //
455    // This function implements the performance optimized assignment of a transpose dense matrix-
456    // sparse matrix multiplication expression to a dense matrix.
457    */
458    template< typename MT  // Type of the target dense matrix
459            , bool SO >    // Storage order of the target dense matrix
460    friend inline auto assign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
461       -> DisableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
462    {
463       BLAZE_FUNCTION_TRACE;
464 
465       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
466       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
467 
468       LT A( serial( rhs.lhs_ ) );  // Evaluation of the left-hand side dense matrix operand
469       RT B( serial( rhs.rhs_ ) );  // Evaluation of the right-hand side sparse matrix operand
470 
471       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
472       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
473       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
474       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
475       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
476       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
477 
478       TDMatSMatMultExpr::selectAssignKernel( *lhs, A, B );
479    }
480    /*! \endcond */
481    //**********************************************************************************************
482 
483    //**Assignment to dense matrices (kernel selection)*********************************************
484    /*! \cond BLAZE_INTERNAL */
485    /*!\brief Selection of the kernel for an assignment of a transpose dense matrix-sparse matrix
486    //        multiplication to a dense matrix (\f$ C=A*B \f$).
487    // \ingroup dense_matrix
488    //
489    // \param C The target left-hand side dense matrix.
490    // \param A The left-hand side multiplication operand.
491    // \param B The right-hand side multiplication operand.
492    // \return void
493    */
494    template< typename MT3    // Type of the left-hand side target matrix
495            , typename MT4    // Type of the left-hand side matrix operand
496            , typename MT5 >  // Type of the right-hand side matrix operand
selectAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)497    static inline void selectAssignKernel( MT3& C, const MT4& A, const MT5& B )
498    {
499       if( C.rows() * C.columns() < TDMATSMATMULT_THRESHOLD )
500          selectSmallAssignKernel( C, A, B );
501       else
502          selectLargeAssignKernel( C, A, B );
503    }
504    /*! \endcond */
505    //**********************************************************************************************
506 
507    //**Default assignment to dense matrices********************************************************
508    /*! \cond BLAZE_INTERNAL */
509    /*!\brief Default assignment of a transpose dense matrix-sparse matrix multiplication
510    //        (\f$ C=A*B \f$).
511    // \ingroup dense_matrix
512    //
513    // \param C The target left-hand side dense matrix.
514    // \param A The left-hand side multiplication operand.
515    // \param B The right-hand side multiplication operand.
516    // \return void
517    //
518    // This function implements the default assignment of a transpose dense matrix-sparse matrix
519    // multiplication expression to a dense matrix. This assign function is used in case the
520    // element type of the target matrix is resizable.
521    */
522    template< typename MT3    // Type of the left-hand side target matrix
523            , typename MT4    // Type of the left-hand side matrix operand
524            , typename MT5 >  // Type of the right-hand side matrix operand
selectDefaultAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)525    static inline void selectDefaultAssignKernel( MT3& C, const MT4& A, const MT5& B )
526    {
527       reset( C );
528 
529       for( size_t j=0UL; j<B.rows(); ++j )
530       {
531          auto element( B.begin(j) );
532          const auto end( B.end(j) );
533 
534          if( IsDiagonal_v<MT4> )
535          {
536             for( ; element!=end; ++element ) {
537                C(j,element->index()) = A(j,j) * element->value();
538             }
539          }
540          else
541          {
542             for( ; element!=end; ++element )
543             {
544                const size_t j1( element->index() );
545 
546                const size_t ibegin( ( IsLower_v<MT4> )
547                                     ?( ( IsStrictlyLower_v<MT4> )
548                                        ?( SYM || HERM || LOW ? max(j1,j+1UL) : j+1UL )
549                                        :( SYM || HERM || LOW ? max(j1,j) : j ) )
550                                     :( SYM || HERM || LOW ? j1 : 0UL ) );
551                const size_t iend( ( IsUpper_v<MT4> )
552                                   ?( ( IsStrictlyUpper_v<MT4> )
553                                      ?( UPP ? min(j1+1UL,j) : j )
554                                      :( UPP ? min(j1,j)+1UL : j+1UL ) )
555                                   :( UPP ? j1+1UL : A.rows() ) );
556 
557                if( ( SYM || HERM || LOW || UPP ) && ( ibegin >= iend ) ) continue;
558                BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
559 
560                for( size_t i=ibegin; i<iend; ++i ) {
561                   if( isDefault( C(i,j1) ) )
562                      C(i,j1) = A(i,j) * element->value();
563                   else
564                      C(i,j1) += A(i,j) * element->value();
565                }
566             }
567          }
568       }
569 
570       if( SYM || HERM ) {
571          for( size_t j=1UL; j<B.columns(); ++j ) {
572             for( size_t i=0UL; i<j; ++i ) {
573                C(i,j) = HERM ? conj( C(j,i) ) : C(j,i);
574             }
575          }
576       }
577    }
578    /*! \endcond */
579    //**********************************************************************************************
580 
581    //**Default assignment to dense matrices (small matrices)***************************************
582    /*! \cond BLAZE_INTERNAL */
583    /*!\brief Default assignment of a small transpose dense matrix-sparse matrix multiplication
584    //        (\f$ C=A*B \f$).
585    // \ingroup dense_matrix
586    //
587    // \param C The target left-hand side dense matrix.
588    // \param A The left-hand side multiplication operand.
589    // \param B The right-hand side multiplication operand.
590    // \return void
591    //
592    // This function relays to the default implementation of the assignment of a transpose dense
593    // matrix-sparse matrix multiplication expression to a dense matrix.
594    */
595    template< typename MT3    // Type of the left-hand side target matrix
596            , typename MT4    // Type of the left-hand side matrix operand
597            , typename MT5 >  // Type of the right-hand side matrix operand
598    static inline auto selectSmallAssignKernel( MT3& C, const MT4& A, const MT5& B )
599       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
600    {
601       selectDefaultAssignKernel( C, A, B );
602    }
603    /*! \endcond */
604    //**********************************************************************************************
605 
606    //**Optimized assignment to dense matrices (small matrices)*************************************
607    /*! \cond BLAZE_INTERNAL */
608    /*!\brief Optimized assignment of a small transpose dense matrix-sparse matrix multiplication
609    //        (\f$ C=A*B \f$).
610    // \ingroup dense_matrix
611    //
612    // \param C The target left-hand side dense matrix.
613    // \param A The left-hand side multiplication operand.
614    // \param B The right-hand side multiplication operand.
615    // \return void
616    //
617    // This function implements the performance optimized assignment of a transpose dense matrix-
618    // sparse matrix multiplication expression to a dense matrix. This assign function is used in
619    // case the element type of the target matrix is not resizable.
620    */
621    template< typename MT3    // Type of the left-hand side target matrix
622            , typename MT4    // Type of the left-hand side matrix operand
623            , typename MT5 >  // Type of the right-hand side matrix operand
624    static inline auto selectSmallAssignKernel( MT3& C, const MT4& A, const MT5& B )
625       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
626    {
627       reset( C );
628 
629       for( size_t j=0UL; j<B.rows(); ++j )
630       {
631          auto element( B.begin(j) );
632          const auto end( B.end(j) );
633 
634          const size_t nonzeros( B.nonZeros(j) );
635          const size_t kpos( prevMultiple( nonzeros, 4UL ) );
636          BLAZE_INTERNAL_ASSERT( kpos <= nonzeros, "Invalid end calculation" );
637 
638          for( size_t k=0UL; k<kpos; k+=4UL )
639          {
640             const size_t j1( element->index() );
641             const ET2    v1( element->value() );
642             ++element;
643             const size_t j2( element->index() );
644             const ET2    v2( element->value() );
645             ++element;
646             const size_t j3( element->index() );
647             const ET2    v3( element->value() );
648             ++element;
649             const size_t j4( element->index() );
650             const ET2    v4( element->value() );
651             ++element;
652 
653             BLAZE_INTERNAL_ASSERT( j1 < j2 && j2 < j3 && j3 < j4, "Invalid sparse matrix index detected" );
654 
655             const size_t ibegin( ( IsLower_v<MT4> )
656                                  ?( ( IsStrictlyLower_v<MT4> )
657                                     ?( SYM || HERM || LOW ? max(j1,j+1UL) : j+1UL )
658                                     :( SYM || HERM || LOW ? max(j1,j) : j ) )
659                                  :( SYM || HERM || LOW ? j1 : 0UL ) );
660             const size_t iend( ( IsUpper_v<MT4> )
661                                ?( ( IsStrictlyUpper_v<MT4> )
662                                   ?( UPP ? min(j4+1UL,j) : j )
663                                   :( UPP ? min(j4,j)+1UL : j+1UL ) )
664                                :( UPP ? j4+1UL : A.rows() ) );
665 
666             if( ( SYM || HERM || LOW || UPP ) && ( ibegin >= iend ) ) continue;
667             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
668 
669             const size_t inum( iend - ibegin );
670             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
671             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
672 
673             size_t i( ibegin );
674 
675             for( ; i<ipos; i+=4UL ) {
676                C(i    ,j1) += A(i    ,j) * v1;
677                C(i+1UL,j1) += A(i+1UL,j) * v1;
678                C(i+2UL,j1) += A(i+2UL,j) * v1;
679                C(i+3UL,j1) += A(i+3UL,j) * v1;
680                C(i    ,j2) += A(i    ,j) * v2;
681                C(i+1UL,j2) += A(i+1UL,j) * v2;
682                C(i+2UL,j2) += A(i+2UL,j) * v2;
683                C(i+3UL,j2) += A(i+3UL,j) * v2;
684                C(i    ,j3) += A(i    ,j) * v3;
685                C(i+1UL,j3) += A(i+1UL,j) * v3;
686                C(i+2UL,j3) += A(i+2UL,j) * v3;
687                C(i+3UL,j3) += A(i+3UL,j) * v3;
688                C(i    ,j4) += A(i    ,j) * v4;
689                C(i+1UL,j4) += A(i+1UL,j) * v4;
690                C(i+2UL,j4) += A(i+2UL,j) * v4;
691                C(i+3UL,j4) += A(i+3UL,j) * v4;
692             }
693             for( ; i<iend; ++i ) {
694                C(i,j1) += A(i,j) * v1;
695                C(i,j2) += A(i,j) * v2;
696                C(i,j3) += A(i,j) * v3;
697                C(i,j4) += A(i,j) * v4;
698             }
699          }
700 
701          for( ; element!=end; ++element )
702          {
703             const size_t j1( element->index() );
704             const ET2    v1( element->value() );
705 
706             const size_t ibegin( ( IsLower_v<MT4> )
707                                  ?( ( IsStrictlyLower_v<MT4> )
708                                     ?( SYM || HERM || LOW ? max(j1,j+1UL) : j+1UL )
709                                     :( SYM || HERM || LOW ? max(j1,j) : j ) )
710                                  :( SYM || HERM || LOW ? j1 : 0UL ) );
711             const size_t iend( ( IsUpper_v<MT4> )
712                                ?( ( IsStrictlyUpper_v<MT4> )
713                                   ?( UPP ? min(j1+1UL,j) : j )
714                                   :( UPP ? min(j1,j)+1UL : j+1UL ) )
715                                :( UPP ? j1+1UL : A.rows() ) );
716 
717             if( ( SYM || HERM || LOW || UPP ) && ( ibegin >= iend ) ) continue;
718             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
719 
720             const size_t inum( iend - ibegin );
721             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
722             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
723 
724             size_t i( ibegin );
725 
726             for( ; i<ipos; i+=4UL ) {
727                C(i    ,j1) += A(i    ,j) * v1;
728                C(i+1UL,j1) += A(i+1UL,j) * v1;
729                C(i+2UL,j1) += A(i+2UL,j) * v1;
730                C(i+3UL,j1) += A(i+3UL,j) * v1;
731             }
732             for( ; i<iend; ++i ) {
733                C(i,j1) += A(i,j) * v1;
734             }
735          }
736       }
737 
738       if( SYM || HERM ) {
739          for( size_t j=1UL; j<B.columns(); ++j ) {
740             for( size_t i=0UL; i<j; ++i ) {
741                C(i,j) = HERM ? conj( C(j,i) ) : C(j,i);
742             }
743          }
744       }
745    }
746    /*! \endcond */
747    //**********************************************************************************************
748 
749    //**Default assignment to dense matrices (large matrices)***************************************
750    /*! \cond BLAZE_INTERNAL */
751    /*!\brief Default assignment of a large transpose dense matrix-sparse matrix multiplication
752    //        (\f$ C=A*B \f$).
753    // \ingroup dense_matrix
754    //
755    // \param C The target left-hand side dense matrix.
756    // \param A The left-hand side multiplication operand.
757    // \param B The right-hand side multiplication operand.
758    // \return void
759    //
760    // This function relays to the default implementation of the assignment of a transpose dense
761    // matrix-sparse matrix multiplication expression to a dense matrix.
762    */
763    template< typename MT3    // Type of the left-hand side target matrix
764            , typename MT4    // Type of the left-hand side matrix operand
765            , typename MT5 >  // Type of the right-hand side matrix operand
766    static inline auto selectLargeAssignKernel( MT3& C, const MT4& A, const MT5& B )
767       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
768    {
769       selectDefaultAssignKernel( C, A, B );
770    }
771    /*! \endcond */
772    //**********************************************************************************************
773 
774    //**Optimized assignment to dense matrices (large matrices)*************************************
775    /*! \cond BLAZE_INTERNAL */
776    /*!\brief Default assignment of a large transpose dense matrix-sparse matrix multiplication
777    //        (\f$ C=A*B \f$).
778    // \ingroup dense_matrix
779    //
780    // \param C The target left-hand side dense matrix.
781    // \param A The left-hand side multiplication operand.
782    // \param B The right-hand side multiplication operand.
783    // \return void
784    //
785    // This function implements the performance optimized assignment of a transpose dense matrix-
786    // sparse matrix multiplication expression to a dense matrix. This assign function is used in
787    // case the element type of the target matrix is not resizable.
788    */
789    template< typename MT3    // Type of the left-hand side target matrix
790            , typename MT4    // Type of the left-hand side matrix operand
791            , typename MT5 >  // Type of the right-hand side matrix operand
792    static inline auto selectLargeAssignKernel( MT3& C, const MT4& A, const MT5& B )
793       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
794    {
795       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( OppositeType_t<MT5> );
796 
797       const ForwardFunctor fwd;
798 
799       const OppositeType_t<MT5> tmp( serial( B ) );
800       assign( C, fwd( A * tmp ) );
801    }
802    /*! \endcond */
803    //**********************************************************************************************
804 
805    //**Assignment to sparse matrices***************************************************************
806    /*! \cond BLAZE_INTERNAL */
807    /*!\brief Assignment of a transpose dense matrix-sparse matrix multiplication to a sparse matrix
808    //        (\f$ A=B*C \f$).
809    // \ingroup dense_matrix
810    //
811    // \param lhs The target left-hand side sparse matrix.
812    // \param rhs The right-hand side multiplication expression to be assigned.
813    // \return void
814    //
815    // This function implements the performance optimized assignment of a transpose dense matrix-
816    // sparse matrix multiplication expression to a sparse matrix.
817    */
818    template< typename MT  // Type of the target sparse matrix
819            , bool SO >    // Storage order of the target sparse matrix
820    friend inline auto assign( SparseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
821       -> DisableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
822    {
823       BLAZE_FUNCTION_TRACE;
824 
825       using TmpType = If_t< SO, ResultType, OppositeType >;
826 
827       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( ResultType );
828       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( OppositeType );
829       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( ResultType );
830       BLAZE_CONSTRAINT_MUST_BE_ROW_MAJOR_MATRIX_TYPE( OppositeType );
831       BLAZE_CONSTRAINT_MATRICES_MUST_HAVE_SAME_STORAGE_ORDER( MT, TmpType );
832       BLAZE_CONSTRAINT_MUST_NOT_REQUIRE_EVALUATION( TmpType );
833 
834       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
835       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
836 
837       const ForwardFunctor fwd;
838 
839       const TmpType tmp( serial( rhs ) );
840       assign( *lhs, fwd( tmp ) );
841    }
842    /*! \endcond */
843    //**********************************************************************************************
844 
845    //**Restructuring assignment********************************************************************
846    /*! \cond BLAZE_INTERNAL */
847    /*!\brief Restructuring assignment of a transpose dense matrix-sparse matrix multiplication
848    //        (\f$ C=A*B \f$).
849    // \ingroup dense_matrix
850    //
851    // \param lhs The target left-hand side matrix.
852    // \param rhs The right-hand side multiplication expression to be assigned.
853    // \return void
854    //
855    // This function implements the symmetry-based restructuring assignment of a transpose dense
856    // matrix-sparse matrix multiplication expression. Due to the explicit application of the
857    // SFINAE principle this function can only be selected by the compiler in case the symmetry
858    // of either of the two matrix operands can be exploited.
859    */
860    template< typename MT  // Type of the target matrix
861            , bool SO >    // Storage order of the target matrix
862    friend inline auto assign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
863       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
864    {
865       BLAZE_FUNCTION_TRACE;
866 
867       BLAZE_CONSTRAINT_MUST_NOT_BE_SYMMETRIC_MATRIX_TYPE( MT );
868 
869       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
870       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
871 
872       const ForwardFunctor fwd;
873 
874       assign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
875    }
876    /*! \endcond */
877    //**********************************************************************************************
878 
879    //**Addition assignment to dense matrices*******************************************************
880    /*! \cond BLAZE_INTERNAL */
881    /*!\brief Addition assignment of a transpose dense matrix-sparse matrix multiplication to
882    //        a dense matrix (\f$ C+=A*B \f$).
883    // \ingroup dense_matrix
884    //
885    // \param lhs The target left-hand side dense matrix.
886    // \param rhs The right-hand side multiplication expression to be added.
887    // \return void
888    //
889    // This function implements the performance optimized addition assignment of a transpose dense
890    // matrix-sparse matrix multiplication expression to a dense matrix.
891    */
892    template< typename MT  // Type of the target dense matrix
893            , bool SO >    // Storage order of the target dense matrix
894    friend inline auto addAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
895       -> DisableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
896    {
897       BLAZE_FUNCTION_TRACE;
898 
899       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
900       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
901 
902       LT A( serial( rhs.lhs_ ) );  // Evaluation of the left-hand side dense matrix operand
903       RT B( serial( rhs.rhs_ ) );  // Evaluation of the right-hand side sparse matrix operand
904 
905       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
906       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
907       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
908       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
909       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
910       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
911 
912       TDMatSMatMultExpr::selectAddAssignKernel( *lhs, A, B );
913    }
914    /*! \endcond */
915    //**********************************************************************************************
916 
917    //**Addition assignment to dense matrices (kernel selection)************************************
918    /*! \cond BLAZE_INTERNAL */
919    /*!\brief Selection of the kernel for an addition assignment of a transpose dense matrix-sparse
920    //        matrix multiplication to a dense matrix (\f$ C+=A*B \f$).
921    // \ingroup dense_matrix
922    //
923    // \param C The target left-hand side dense matrix.
924    // \param A The left-hand side multiplication operand.
925    // \param B The right-hand side multiplication operand.
926    // \return void
927    */
928    template< typename MT3    // Type of the left-hand side target matrix
929            , typename MT4    // Type of the left-hand side matrix operand
930            , typename MT5 >  // Type of the right-hand side matrix operand
selectAddAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)931    static inline void selectAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
932    {
933       if( C.rows() * C.columns() < TDMATSMATMULT_THRESHOLD )
934          selectSmallAddAssignKernel( C, A, B );
935       else
936          selectLargeAddAssignKernel( C, A, B );
937    }
938    /*! \endcond */
939    //**********************************************************************************************
940 
941    //**Default addition assignment to dense matrices***********************************************
942    /*! \cond BLAZE_INTERNAL */
943    /*!\brief Default addition assignment of a transpose dense matrix-sparse matrix multiplication
944    //        (\f$ C+=A*B \f$).
945    // \ingroup dense_matrix
946    //
947    // \param C The target left-hand side dense matrix.
948    // \param A The left-hand side multiplication operand.
949    // \param B The right-hand side multiplication operand.
950    // \return void
951    //
952    // This function implements the default addition assignment of a transpose dense matrix-sparse
953    // matrix multiplication expression to a dense matrix. This assign function is use in case the
954    // element type of the target matrix is resizable.
955    */
956    template< typename MT3    // Type of the left-hand side target matrix
957            , typename MT4    // Type of the left-hand side matrix operand
958            , typename MT5 >  // Type of the right-hand side matrix operand
selectDefaultAddAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)959    static inline void selectDefaultAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
960    {
961       size_t i( 0UL );
962 
963       for( size_t j=0UL; j<B.rows(); ++j )
964       {
965          auto element( B.begin(j) );
966          const auto end( B.end(j) );
967 
968          if( IsDiagonal_v<MT4> )
969          {
970             for( ; element!=end; ++element ) {
971                C(j,element->index()) += A(j,j) * element->value();
972             }
973          }
974          else
975          {
976             for( ; element!=end; ++element )
977             {
978                const size_t j1( element->index() );
979 
980                const size_t ibegin( ( IsLower_v<MT4> )
981                                     ?( ( IsStrictlyLower_v<MT4> )
982                                        ?( LOW ? max(j1,j+1UL) : j+1UL )
983                                        :( LOW ? max(j1,j) : j ) )
984                                     :( LOW ? j1 : 0UL ) );
985                const size_t iend( ( IsUpper_v<MT4> )
986                                   ?( ( IsStrictlyUpper_v<MT4> )
987                                      ?( UPP ? min(j1+1UL,j) : j )
988                                      :( UPP ? min(j1,j)+1UL : j+1UL ) )
989                                   :( UPP ? j1+1UL : A.rows() ) );
990 
991                if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
992                BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
993 
994                const size_t inum( iend - ibegin );
995                const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
996                BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
997 
998                for( i=ibegin; i<ipos; i+=4UL ) {
999                   C(i    ,j1) += A(i    ,j) * element->value();
1000                   C(i+1UL,j1) += A(i+1UL,j) * element->value();
1001                   C(i+2UL,j1) += A(i+2UL,j) * element->value();
1002                   C(i+3UL,j1) += A(i+3UL,j) * element->value();
1003                }
1004                for( ; i<iend; ++i ) {
1005                   C(i,j1) += A(i,j) * element->value();
1006                }
1007             }
1008          }
1009       }
1010    }
1011    /*! \endcond */
1012    //**********************************************************************************************
1013 
1014    //**Default addition assignment to dense matrices (small matrices)******************************
1015    /*! \cond BLAZE_INTERNAL */
1016    /*!\brief Default addition assignment of a small transpose dense matrix-sparse matrix
1017    //        multiplication (\f$ C+=A*B \f$).
1018    // \ingroup dense_matrix
1019    //
1020    // \param C The target left-hand side dense matrix.
1021    // \param A The left-hand side multiplication operand.
1022    // \param B The right-hand side multiplication operand.
1023    // \return void
1024    //
1025    // This function relays to the default implementation of the addition assignment of a transpose
1026    // dense matrix-sparse matrix multiplication expression to a dense matrix.
1027    */
1028    template< typename MT3    // Type of the left-hand side target matrix
1029            , typename MT4    // Type of the left-hand side matrix operand
1030            , typename MT5 >  // Type of the right-hand side matrix operand
1031    static inline auto selectSmallAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
1032       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
1033    {
1034       selectDefaultAddAssignKernel( C, A, B );
1035    }
1036    /*! \endcond */
1037    //**********************************************************************************************
1038 
1039    //**Optimized addition assignment to dense matrices (small matrices)****************************
1040    /*! \cond BLAZE_INTERNAL */
1041    /*!\brief Optimized addition assignment of a small transpose dense matrix-sparse matrix
1042    //        multiplication (\f$ C+=A*B \f$).
1043    // \ingroup dense_matrix
1044    //
1045    // \param C The target left-hand side dense matrix.
1046    // \param A The left-hand side multiplication operand.
1047    // \param B The right-hand side multiplication operand.
1048    // \return void
1049    //
1050    // This function implements the performance optimized addition assignment of a transpose dense
1051    // matrix-sparse matrix multiplication expression to a dense matrix. This assign function is
1052    // used in case the element type of the target matrix is not resizable.
1053    */
1054    template< typename MT3    // Type of the left-hand side target matrix
1055            , typename MT4    // Type of the left-hand side matrix operand
1056            , typename MT5 >  // Type of the right-hand side matrix operand
1057    static inline auto selectSmallAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
1058       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
1059    {
1060       for( size_t j=0UL; j<B.rows(); ++j )
1061       {
1062          auto element( B.begin(j) );
1063          const auto end( B.end(j) );
1064 
1065          const size_t nonzeros( B.nonZeros(j) );
1066          const size_t kpos( prevMultiple( nonzeros, 4UL ) );
1067          BLAZE_INTERNAL_ASSERT( kpos <= nonzeros, "Invalid end calculation" );
1068 
1069          for( size_t k=0UL; k<kpos; k+=4UL )
1070          {
1071             const size_t j1( element->index() );
1072             const ET2    v1( element->value() );
1073             ++element;
1074             const size_t j2( element->index() );
1075             const ET2    v2( element->value() );
1076             ++element;
1077             const size_t j3( element->index() );
1078             const ET2    v3( element->value() );
1079             ++element;
1080             const size_t j4( element->index() );
1081             const ET2    v4( element->value() );
1082             ++element;
1083 
1084             BLAZE_INTERNAL_ASSERT( j1 < j2 && j2 < j3 && j3 < j4, "Invalid sparse matrix index detected" );
1085 
1086             const size_t ibegin( ( IsLower_v<MT4> )
1087                                  ?( ( IsStrictlyLower_v<MT4> )
1088                                     ?( LOW ? max(j1,j+1UL) : j+1UL )
1089                                     :( LOW ? max(j1,j) : j ) )
1090                                  :( LOW ? j1 : 0UL ) );
1091             const size_t iend( ( IsUpper_v<MT4> )
1092                                ?( ( IsStrictlyUpper_v<MT4> )
1093                                   ?( UPP ? min(j4+1UL,j) : j )
1094                                   :( UPP ? min(j4,j)+1UL : j+1UL ) )
1095                                :( UPP ? j4+1UL : A.rows() ) );
1096 
1097             if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
1098             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
1099 
1100             const size_t inum( iend - ibegin );
1101             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
1102             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
1103 
1104             size_t i( ibegin );
1105 
1106             for( i=ibegin; i<ipos; i+=4UL ) {
1107                C(i    ,j1) += A(i    ,j) * v1;
1108                C(i+1UL,j1) += A(i+1UL,j) * v1;
1109                C(i+2UL,j1) += A(i+2UL,j) * v1;
1110                C(i+3UL,j1) += A(i+3UL,j) * v1;
1111                C(i    ,j2) += A(i    ,j) * v2;
1112                C(i+1UL,j2) += A(i+1UL,j) * v2;
1113                C(i+2UL,j2) += A(i+2UL,j) * v2;
1114                C(i+3UL,j2) += A(i+3UL,j) * v2;
1115                C(i    ,j3) += A(i    ,j) * v3;
1116                C(i+1UL,j3) += A(i+1UL,j) * v3;
1117                C(i+2UL,j3) += A(i+2UL,j) * v3;
1118                C(i+3UL,j3) += A(i+3UL,j) * v3;
1119                C(i    ,j4) += A(i    ,j) * v4;
1120                C(i+1UL,j4) += A(i+1UL,j) * v4;
1121                C(i+2UL,j4) += A(i+2UL,j) * v4;
1122                C(i+3UL,j4) += A(i+3UL,j) * v4;
1123             }
1124             for( ; i<iend; ++i ) {
1125                C(i,j1) += A(i,j) * v1;
1126                C(i,j2) += A(i,j) * v2;
1127                C(i,j3) += A(i,j) * v3;
1128                C(i,j4) += A(i,j) * v4;
1129             }
1130          }
1131 
1132          for( ; element!=end; ++element )
1133          {
1134             const size_t j1( element->index() );
1135             const ET2    v1( element->value() );
1136 
1137             const size_t ibegin( ( IsLower_v<MT4> )
1138                                  ?( ( IsStrictlyLower_v<MT4> )
1139                                     ?( LOW ? max(j1,j+1UL) : j+1UL )
1140                                     :( LOW ? max(j1,j) : j ) )
1141                                  :( LOW ? j1 : 0UL ) );
1142             const size_t iend( ( IsUpper_v<MT4> )
1143                                ?( ( IsStrictlyUpper_v<MT4> )
1144                                   ?( UPP ? min(j1+1UL,j) : j )
1145                                   :( UPP ? min(j1,j)+1UL : j+1UL ) )
1146                                :( UPP ? j1+1UL : A.rows() ) );
1147 
1148             if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
1149             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
1150 
1151             const size_t inum( iend - ibegin );
1152             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
1153             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
1154 
1155             size_t i( ibegin );
1156 
1157             for( ; i<ipos; i+=4UL ) {
1158                C(i    ,j1) += A(i    ,j) * v1;
1159                C(i+1UL,j1) += A(i+1UL,j) * v1;
1160                C(i+2UL,j1) += A(i+2UL,j) * v1;
1161                C(i+3UL,j1) += A(i+3UL,j) * v1;
1162             }
1163             for( ; i<iend; ++i ) {
1164                C(i,j1) += A(i,j) * v1;
1165             }
1166          }
1167       }
1168    }
1169    /*! \endcond */
1170    //**********************************************************************************************
1171 
1172    //**Default addition assignment to dense matrices (large matrices)******************************
1173    /*! \cond BLAZE_INTERNAL */
1174    /*!\brief Default addition assignment of a large transpose dense matrix-sparse matrix
1175    //        multiplication (\f$ C+=A*B \f$).
1176    // \ingroup dense_matrix
1177    //
1178    // \param C The target left-hand side dense matrix.
1179    // \param A The left-hand side multiplication operand.
1180    // \param B The right-hand side multiplication operand.
1181    // \return void
1182    //
1183    // This function relays to the default implementation of the addition assignment of a transpose
1184    // dense matrix-sparse matrix multiplication expression to a dense matrix.
1185    */
1186    template< typename MT3    // Type of the left-hand side target matrix
1187            , typename MT4    // Type of the left-hand side matrix operand
1188            , typename MT5 >  // Type of the right-hand side matrix operand
1189    static inline auto selectLargeAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
1190       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
1191    {
1192       selectDefaultAddAssignKernel( C, A, B );
1193    }
1194    /*! \endcond */
1195    //**********************************************************************************************
1196 
1197    //**Optimized addition assignment to dense matrices (large matrices)****************************
1198    /*! \cond BLAZE_INTERNAL */
1199    /*!\brief Default addition assignment of a large transpose dense matrix-sparse matrix
1200    //        multiplication (\f$ C+=A*B \f$).
1201    // \ingroup dense_matrix
1202    //
1203    // \param C The target left-hand side dense matrix.
1204    // \param A The left-hand side multiplication operand.
1205    // \param B The right-hand side multiplication operand.
1206    // \return void
1207    //
1208    // This function implements the performance optimized addition assignment of a transpose dense
1209    // matrix-sparse matrix multiplication expression to a dense matrix. This assign function is
1210    // used in case the element type of the target matrix is not resizable.
1211    */
1212    template< typename MT3    // Type of the left-hand side target matrix
1213            , typename MT4    // Type of the left-hand side matrix operand
1214            , typename MT5 >  // Type of the right-hand side matrix operand
1215    static inline auto selectLargeAddAssignKernel( MT3& C, const MT4& A, const MT5& B )
1216       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
1217    {
1218       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( OppositeType_t<MT5> );
1219 
1220       const ForwardFunctor fwd;
1221 
1222       const OppositeType_t<MT5> tmp( serial( B ) );
1223       addAssign( C, fwd( A * tmp ) );
1224    }
1225    /*! \endcond */
1226    //**********************************************************************************************
1227 
1228    //**Restructuring addition assignment***********************************************************
1229    /*! \cond BLAZE_INTERNAL */
1230    /*!\brief Restructuring addition assignment of a transpose dense matrix-sparse matrix
1231    //        multiplication (\f$ C+=A*B \f$).
1232    // \ingroup dense_matrix
1233    //
1234    // \param lhs The target left-hand side matrix.
1235    // \param rhs The right-hand side multiplication expression to be added.
1236    // \return void
1237    //
1238    // This function implements the symmetry-based restructuring addition assignment of a transpose
1239    // dense matrix-sparse matrix multiplication expression. Due to the explicit application of the
1240    // SFINAE principle this function can only be selected by the compiler in case the symmetry of
1241    // either of the two matrix operands can be exploited.
1242    */
1243    template< typename MT  // Type of the target matrix
1244            , bool SO >    // Storage order of the target matrix
1245    friend inline auto addAssign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1246       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1247    {
1248       BLAZE_FUNCTION_TRACE;
1249 
1250       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1251       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1252 
1253       const ForwardFunctor fwd;
1254 
1255       addAssign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
1256    }
1257    /*! \endcond */
1258    //**********************************************************************************************
1259 
1260    //**Addition assignment to sparse matrices******************************************************
1261    // No special implementation for the addition assignment to sparse matrices.
1262    //**********************************************************************************************
1263 
1264    //**Subtraction assignment to dense matrices****************************************************
1265    /*! \cond BLAZE_INTERNAL */
1266    /*!\brief Subtraction assignment of a transpose dense matrix-sparse matrix multiplication to a
1267    //        dense matrix (\f$ C-=A*B \f$).
1268    // \ingroup dense_matrix
1269    //
1270    // \param lhs The target left-hand side dense matrix.
1271    // \param rhs The right-hand side multiplication expression to be subtracted.
1272    // \return void
1273    //
1274    // This function implements the performance optimized subtraction assignment of a transpose
1275    // dense matrix-sparse matrix multiplication expression to a dense matrix.
1276    */
1277    template< typename MT  // Type of the target dense matrix
1278            , bool SO >    // Storage order of the target dense matrix
1279    friend inline auto subAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1280       -> DisableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1281    {
1282       BLAZE_FUNCTION_TRACE;
1283 
1284       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1285       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1286 
1287       LT A( serial( rhs.lhs_ ) );  // Evaluation of the left-hand side dense matrix operand
1288       RT B( serial( rhs.rhs_ ) );  // Evaluation of the right-hand side sparse matrix operand
1289 
1290       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
1291       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
1292       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
1293       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
1294       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
1295       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
1296 
1297       TDMatSMatMultExpr::selectSubAssignKernel( *lhs, A, B );
1298    }
1299    /*! \endcond */
1300    //**********************************************************************************************
1301 
1302    //**Subtraction assignment to dense matrices (kernel selection)*********************************
1303    /*! \cond BLAZE_INTERNAL */
1304    /*!\brief Selection of the kernel for a subtraction assignment of a transpose dense matrix-
1305    //        sparse matrix multiplication to a dense matrix (\f$ C-=A*B \f$).
1306    // \ingroup dense_matrix
1307    //
1308    // \param C The target left-hand side dense matrix.
1309    // \param A The left-hand side multiplication operand.
1310    // \param B The right-hand side multiplication operand.
1311    // \return void
1312    */
1313    template< typename MT3    // Type of the left-hand side target matrix
1314            , typename MT4    // Type of the left-hand side matrix operand
1315            , typename MT5 >  // Type of the right-hand side matrix operand
selectSubAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)1316    static inline void selectSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1317    {
1318       if( C.rows() * C.columns() < TDMATSMATMULT_THRESHOLD )
1319          selectSmallSubAssignKernel( C, A, B );
1320       else
1321          selectLargeSubAssignKernel( C, A, B );
1322    }
1323    /*! \endcond */
1324    //**********************************************************************************************
1325 
1326    //**Default subtraction assignment to dense matrices********************************************
1327    /*! \cond BLAZE_INTERNAL */
1328    /*!\brief Default subtraction assignment of a transpose dense matrix-sparse matrix
1329    //        multiplication (\f$ C-=A*B \f$).
1330    // \ingroup dense_matrix
1331    //
1332    // \param C The target left-hand side dense matrix.
1333    // \param A The left-hand side multiplication operand.
1334    // \param B The right-hand side multiplication operand.
1335    // \return void
1336    //
1337    // This function implements the default subtraction assignment of a transpose dense matrix-
1338    // sparse matrix multiplication expression to a dense matrix. This assign function is used
1339    // in case the element type of the target matrix is resizable.
1340    */
1341    template< typename MT3    // Type of the left-hand side target matrix
1342            , typename MT4    // Type of the left-hand side matrix operand
1343            , typename MT5 >  // Type of the right-hand side matrix operand
selectDefaultSubAssignKernel(MT3 & C,const MT4 & A,const MT5 & B)1344    static inline void selectDefaultSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1345    {
1346       size_t i( 0UL );
1347 
1348       for( size_t j=0UL; j<B.rows(); ++j )
1349       {
1350          auto element( B.begin(j) );
1351          const auto end( B.end(j) );
1352 
1353          if( IsDiagonal_v<MT4> )
1354          {
1355             for( ; element!=end; ++element ) {
1356                C(j,element->index()) -= A(j,j) * element->value();
1357             }
1358          }
1359          else
1360          {
1361             for( ; element!=end; ++element )
1362             {
1363                const size_t j1( element->index() );
1364 
1365                const size_t ibegin( ( IsLower_v<MT4> )
1366                                     ?( ( IsStrictlyLower_v<MT4> )
1367                                        ?( LOW ? max(j1,j+1UL) : j+1UL )
1368                                        :( LOW ? max(j1,j) : j ) )
1369                                     :( LOW ? j1 : 0UL ) );
1370                const size_t iend( ( IsUpper_v<MT4> )
1371                                   ?( ( IsStrictlyUpper_v<MT4> )
1372                                      ?( UPP ? min(j1+1UL,j) : j )
1373                                      :( UPP ? min(j1,j)+1UL : j+1UL ) )
1374                                   :( UPP ? j1+1UL : A.rows() ) );
1375 
1376                if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
1377                BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
1378 
1379                const size_t inum( iend - ibegin );
1380                const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
1381                BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
1382 
1383                for( i=ibegin; i<ipos; i+=4UL ) {
1384                   C(i    ,j1) -= A(i    ,j) * element->value();
1385                   C(i+1UL,j1) -= A(i+1UL,j) * element->value();
1386                   C(i+2UL,j1) -= A(i+2UL,j) * element->value();
1387                   C(i+3UL,j1) -= A(i+3UL,j) * element->value();
1388                }
1389                for( ; i<iend; ++i ) {
1390                   C(i,j1) -= A(i,j) * element->value();
1391                }
1392             }
1393          }
1394       }
1395    }
1396    /*! \endcond */
1397    //**********************************************************************************************
1398 
1399    //**Default subtraction assignment to dense matrices (small matrices)***************************
1400    /*! \cond BLAZE_INTERNAL */
1401    /*!\brief Default subtraction assignment of a small transpose dense matrix-sparse matrix
1402    //        multiplication (\f$ C-=A*B \f$).
1403    // \ingroup dense_matrix
1404    //
1405    // \param C The target left-hand side dense matrix.
1406    // \param A The left-hand side multiplication operand.
1407    // \param B The right-hand side multiplication operand.
1408    // \return void
1409    //
1410    // This function relays to the default implementation of the subtraction assignment of a
1411    // transpose dense matrix-sparse matrix multiplication expression to a dense matrix.
1412    */
1413    template< typename MT3    // Type of the left-hand side target matrix
1414            , typename MT4    // Type of the left-hand side matrix operand
1415            , typename MT5 >  // Type of the right-hand side matrix operand
1416    static inline auto selectSmallSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1417       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
1418    {
1419       selectDefaultSubAssignKernel( C, A, B );
1420    }
1421    /*! \endcond */
1422    //**********************************************************************************************
1423 
1424    //**Optimized subtraction assignment to dense matrices (small matrices)*************************
1425    /*! \cond BLAZE_INTERNAL */
1426    /*!\brief Optimized subtraction assignment of a small transpose dense matrix-sparse matrix
1427    //        multiplication (\f$ C-=A*B \f$).
1428    // \ingroup dense_matrix
1429    //
1430    // \param C The target left-hand side dense matrix.
1431    // \param A The left-hand side multiplication operand.
1432    // \param B The right-hand side multiplication operand.
1433    // \return void
1434    //
1435    // This function implements the performance optimized subtraction assignment of a transpose
1436    // dense matrix-sparse matrix multiplication expression to a dense matrix. This assign function
1437    // is used in case the element type of the target matrix is not resizable.
1438    */
1439    template< typename MT3    // Type of the left-hand side target matrix
1440            , typename MT4    // Type of the left-hand side matrix operand
1441            , typename MT5 >  // Type of the right-hand side matrix operand
1442    static inline auto selectSmallSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1443       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
1444    {
1445       for( size_t j=0UL; j<B.rows(); ++j )
1446       {
1447          auto element( B.begin(j) );
1448          const auto end( B.end(j) );
1449 
1450          const size_t nonzeros( B.nonZeros(j) );
1451          const size_t kpos( prevMultiple( nonzeros, 4UL ) );
1452          BLAZE_INTERNAL_ASSERT( kpos <= nonzeros, "Invalid end calculation" );
1453 
1454          for( size_t k=0UL; k<kpos; k+=4UL )
1455          {
1456             const size_t j1( element->index() );
1457             const ET2    v1( element->value() );
1458             ++element;
1459             const size_t j2( element->index() );
1460             const ET2    v2( element->value() );
1461             ++element;
1462             const size_t j3( element->index() );
1463             const ET2    v3( element->value() );
1464             ++element;
1465             const size_t j4( element->index() );
1466             const ET2    v4( element->value() );
1467             ++element;
1468 
1469             BLAZE_INTERNAL_ASSERT( j1 < j2 && j2 < j3 && j3 < j4, "Invalid sparse matrix index detected" );
1470 
1471             const size_t ibegin( ( IsLower_v<MT4> )
1472                                  ?( ( IsStrictlyLower_v<MT4> )
1473                                     ?( LOW ? max(j1,j+1UL) : j+1UL )
1474                                     :( LOW ? max(j1,j) : j ) )
1475                                  :( LOW ? j1 : 0UL ) );
1476             const size_t iend( ( IsUpper_v<MT4> )
1477                                ?( ( IsStrictlyUpper_v<MT4> )
1478                                   ?( UPP ? min(j4+1UL,j) : j )
1479                                   :( UPP ? min(j4,j)+1UL : j+1UL ) )
1480                                :( UPP ? j4+1UL : A.rows() ) );
1481 
1482             if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
1483             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
1484 
1485             const size_t inum( iend - ibegin );
1486             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
1487             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
1488 
1489             size_t i( ibegin );
1490 
1491             for( ; i<ipos; i+=4UL ) {
1492                C(i    ,j1) -= A(i    ,j) * v1;
1493                C(i+1UL,j1) -= A(i+1UL,j) * v1;
1494                C(i+2UL,j1) -= A(i+2UL,j) * v1;
1495                C(i+3UL,j1) -= A(i+3UL,j) * v1;
1496                C(i    ,j2) -= A(i    ,j) * v2;
1497                C(i+1UL,j2) -= A(i+1UL,j) * v2;
1498                C(i+2UL,j2) -= A(i+2UL,j) * v2;
1499                C(i+3UL,j2) -= A(i+3UL,j) * v2;
1500                C(i    ,j3) -= A(i    ,j) * v3;
1501                C(i+1UL,j3) -= A(i+1UL,j) * v3;
1502                C(i+2UL,j3) -= A(i+2UL,j) * v3;
1503                C(i+3UL,j3) -= A(i+3UL,j) * v3;
1504                C(i    ,j4) -= A(i    ,j) * v4;
1505                C(i+1UL,j4) -= A(i+1UL,j) * v4;
1506                C(i+2UL,j4) -= A(i+2UL,j) * v4;
1507                C(i+3UL,j4) -= A(i+3UL,j) * v4;
1508             }
1509             for( ; i<iend; ++i ) {
1510                C(i,j1) -= A(i,j) * v1;
1511                C(i,j2) -= A(i,j) * v2;
1512                C(i,j3) -= A(i,j) * v3;
1513                C(i,j4) -= A(i,j) * v4;
1514             }
1515          }
1516 
1517          for( ; element!=end; ++element )
1518          {
1519             const size_t j1( element->index() );
1520             const ET2    v1( element->value() );
1521 
1522             const size_t ibegin( ( IsLower_v<MT4> )
1523                                  ?( ( IsStrictlyLower_v<MT4> )
1524                                     ?( LOW ? max(j1,j+1UL) : j+1UL )
1525                                     :( LOW ? max(j1,j) : j ) )
1526                                  :( LOW ? j1 : 0UL ) );
1527             const size_t iend( ( IsUpper_v<MT4> )
1528                                ?( ( IsStrictlyUpper_v<MT4> )
1529                                   ?( UPP ? min(j1+1UL,j) : j )
1530                                   :( UPP ? min(j1,j)+1UL : j+1UL ) )
1531                                :( UPP ? j1+1UL : A.rows() ) );
1532 
1533             if( ( LOW || UPP ) && ( ibegin >= iend ) ) continue;
1534             BLAZE_INTERNAL_ASSERT( ibegin <= iend, "Invalid loop indices detected" );
1535 
1536             const size_t inum( iend - ibegin );
1537             const size_t ipos( ibegin + prevMultiple( inum, 4UL ) );
1538             BLAZE_INTERNAL_ASSERT( ipos <= ibegin+inum, "Invalid end calculation" );
1539 
1540             size_t i( ibegin );
1541 
1542             for( ; i<ipos; i+=4UL ) {
1543                C(i    ,j1) -= A(i    ,j) * v1;
1544                C(i+1UL,j1) -= A(i+1UL,j) * v1;
1545                C(i+2UL,j1) -= A(i+2UL,j) * v1;
1546                C(i+3UL,j1) -= A(i+3UL,j) * v1;
1547             }
1548             for( ; i<iend; ++i ) {
1549                C(i,j1) -= A(i,j) * v1;
1550             }
1551          }
1552       }
1553    }
1554    /*! \endcond */
1555    //**********************************************************************************************
1556 
1557    //**Default subtraction assignment to dense matrices (large matrices)***************************
1558    /*! \cond BLAZE_INTERNAL */
1559    /*!\brief Default subtraction assignment of a large transpose dense matrix-sparse matrix
1560    //        multiplication (\f$ C-=A*B \f$).
1561    // \ingroup dense_matrix
1562    //
1563    // \param C The target left-hand side dense matrix.
1564    // \param A The left-hand side multiplication operand.
1565    // \param B The right-hand side multiplication operand.
1566    // \return void
1567    //
1568    // This function relays to the default implementation of the subtraction assignment of a
1569    // transpose dense matrix-sparse matrix multiplication expression to a dense matrix.
1570    */
1571    template< typename MT3    // Type of the left-hand side target matrix
1572            , typename MT4    // Type of the left-hand side matrix operand
1573            , typename MT5 >  // Type of the right-hand side matrix operand
1574    static inline auto selectLargeSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1575       -> EnableIf_t< UseDefaultKernel_v<MT3,MT4,MT5> >
1576    {
1577       selectDefaultSubAssignKernel( C, A, B );
1578    }
1579    /*! \endcond */
1580    //**********************************************************************************************
1581 
1582    //**Optimized subtraction assignment to dense matrices (large matrices)*************************
1583    /*! \cond BLAZE_INTERNAL */
1584    /*!\brief Default subtraction assignment of a large transpose dense matrix-sparse matrix
1585    //        multiplication (\f$ C-=A*B \f$).
1586    // \ingroup dense_matrix
1587    //
1588    // \param C The target left-hand side dense matrix.
1589    // \param A The left-hand side multiplication operand.
1590    // \param B The right-hand side multiplication operand.
1591    // \return void
1592    //
1593    // This function implements the performance optimized subtraction assignment of a transpose
1594    // dense matrix-sparse matrix multiplication expression to a dense matrix. This assign function
1595    // is used in case the element type of the target matrix is not resizable.
1596    */
1597    template< typename MT3    // Type of the left-hand side target matrix
1598            , typename MT4    // Type of the left-hand side matrix operand
1599            , typename MT5 >  // Type of the right-hand side matrix operand
1600    static inline auto selectLargeSubAssignKernel( MT3& C, const MT4& A, const MT5& B )
1601       -> EnableIf_t< UseOptimizedKernel_v<MT3,MT4,MT5> >
1602    {
1603       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( OppositeType_t<MT5> );
1604 
1605       const ForwardFunctor fwd;
1606 
1607       const OppositeType_t<MT5> tmp( serial( B ) );
1608       subAssign( C, fwd( A * tmp ) );
1609    }
1610    /*! \endcond */
1611    //**********************************************************************************************
1612 
1613    //**Restructuring subtraction assignment********************************************************
1614    /*! \cond BLAZE_INTERNAL */
1615    /*!\brief Restructuring subtraction assignment of a transpose dense matrix-sparse matrix
1616    //        multiplication (\f$ C-=A*B \f$).
1617    // \ingroup dense_matrix
1618    //
1619    // \param lhs The target left-hand side matrix.
1620    // \param rhs The right-hand side multiplication expression to be subtracted.
1621    // \return void
1622    //
1623    // This function implements the symmetry-based restructuring subtraction assignment of a
1624    // transpose dense matrix-sparse matrix multiplication expression. Due to the explicit
1625    // application of the SFINAE principle this function can only be selected by the compiler
1626    // in case the symmetry of either of the two matrix operands can be exploited.
1627    */
1628    template< typename MT  // Type of the target matrix
1629            , bool SO >    // Storage order of the target matrix
1630    friend inline auto subAssign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1631       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1632    {
1633       BLAZE_FUNCTION_TRACE;
1634 
1635       BLAZE_CONSTRAINT_MUST_NOT_BE_SYMMETRIC_MATRIX_TYPE( MT );
1636 
1637       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1638       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1639 
1640       const ForwardFunctor fwd;
1641 
1642       subAssign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
1643    }
1644    /*! \endcond */
1645    //**********************************************************************************************
1646 
1647    //**Subtraction assignment to sparse matrices***************************************************
1648    // No special implementation for the subtraction assignment to sparse matrices.
1649    //**********************************************************************************************
1650 
1651    //**Schur product assignment to dense matrices**************************************************
1652    /*! \cond BLAZE_INTERNAL */
1653    /*!\brief Schur product assignment of a transpose dense matrix-sparse matrix multiplication
1654    //        to a dense matrix (\f$ C\circ=A*B \f$).
1655    // \ingroup dense_matrix
1656    //
1657    // \param lhs The target left-hand side dense matrix.
1658    // \param rhs The right-hand side multiplication expression for the Schur product.
1659    // \return void
1660    //
1661    // This function implements the performance optimized Schur product assignment of a transpose
1662    // dense matrix-sparse matrix multiplication expression to a dense matrix.
1663    */
1664    template< typename MT  // Type of the target dense matrix
1665            , bool SO >    // Storage order of the target dense matrix
schurAssign(DenseMatrix<MT,SO> & lhs,const TDMatSMatMultExpr & rhs)1666    friend inline void schurAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1667    {
1668       BLAZE_FUNCTION_TRACE;
1669 
1670       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( ResultType );
1671       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( ResultType );
1672       BLAZE_CONSTRAINT_MUST_NOT_REQUIRE_EVALUATION( ResultType );
1673 
1674       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1675       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1676 
1677       const ResultType tmp( serial( rhs ) );
1678       schurAssign( *lhs, tmp );
1679    }
1680    /*! \endcond */
1681    //**********************************************************************************************
1682 
1683    //**Schur product assignment to sparse matrices*************************************************
1684    // No special implementation for the Schur product assignment to sparse matrices.
1685    //**********************************************************************************************
1686 
1687    //**Multiplication assignment to dense matrices*************************************************
1688    // No special implementation for the multiplication assignment to dense matrices.
1689    //**********************************************************************************************
1690 
1691    //**Multiplication assignment to sparse matrices************************************************
1692    // No special implementation for the multiplication assignment to sparse matrices.
1693    //**********************************************************************************************
1694 
1695    //**SMP assignment to dense matrices************************************************************
1696    /*! \cond BLAZE_INTERNAL */
1697    /*!\brief Assignment of a transpose dense matrix-sparse matrix multiplication to a dense matrix
1698    //        (\f$ C=A*B \f$).
1699    // \ingroup dense_matrix
1700    //
1701    // \param lhs The target left-hand side dense matrix.
1702    // \param rhs The right-hand side multiplication expression to be assigned.
1703    // \return void
1704    //
1705    // This function implements the performance optimized assignment of a transpose dense matrix-
1706    // sparse matrix multiplication expression to a dense matrix. Due to the explicit application
1707    // of the SFINAE principle this function can only be selected by the compiler in case either
1708    // of the two matrix operands requires an intermediate evaluation and no symmetry can be
1709    // exploited.
1710    */
1711    template< typename MT  // Type of the target dense matrix
1712            , bool SO >    // Storage order of the target dense matrix
1713    friend inline auto smpAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1714       -> EnableIf_t< IsEvaluationRequired_v<MT,MT1,MT2> >
1715    {
1716       BLAZE_FUNCTION_TRACE;
1717 
1718       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1719       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1720 
1721       LT A( rhs.lhs_ );  // Evaluation of the left-hand side dense matrix operand
1722       RT B( rhs.rhs_ );  // Evaluation of the right-hand side sparse matrix operand
1723 
1724       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
1725       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
1726       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
1727       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
1728       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
1729       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
1730 
1731       smpAssign( *lhs, A * B );
1732    }
1733    /*! \endcond */
1734    //**********************************************************************************************
1735 
1736    //**SMP assignment to sparse matrices***********************************************************
1737    /*! \cond BLAZE_INTERNAL */
1738    /*!\brief SMP assignment of a transpose dense matrix-sparse matrix multiplication to a sparse
1739    //        matrix (\f$ A=B*C \f$).
1740    // \ingroup dense_matrix
1741    //
1742    // \param lhs The target left-hand side sparse matrix.
1743    // \param rhs The right-hand side multiplication expression to be assigned.
1744    // \return void
1745    //
1746    // This function implements the performance optimized SMP assignment of a transpose dense
1747    // matrix-sparse matrix multiplication expression to a sparse matrix. Due to the explicit
1748    // application of the SFINAE principle this function can only be selected by the compiler
1749    // in case either of the two matrix operands requires an intermediate evaluation and no
1750    // symmetry can be exploited.
1751    */
1752    template< typename MT  // Type of the target sparse matrix
1753            , bool SO >    // Storage order of the target sparse matrix
1754    friend inline auto smpAssign( SparseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1755       -> EnableIf_t< IsEvaluationRequired_v<MT,MT1,MT2> >
1756    {
1757       BLAZE_FUNCTION_TRACE;
1758 
1759       using TmpType = If_t< SO, ResultType, OppositeType >;
1760 
1761       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( ResultType );
1762       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( OppositeType );
1763       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( ResultType );
1764       BLAZE_CONSTRAINT_MUST_BE_ROW_MAJOR_MATRIX_TYPE( OppositeType );
1765       BLAZE_CONSTRAINT_MATRICES_MUST_HAVE_SAME_STORAGE_ORDER( MT, TmpType );
1766       BLAZE_CONSTRAINT_MUST_NOT_REQUIRE_EVALUATION( TmpType );
1767 
1768       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1769       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1770 
1771       const ForwardFunctor fwd;
1772 
1773       const TmpType tmp( rhs );
1774       smpAssign( *lhs, fwd( tmp ) );
1775    }
1776    /*! \endcond */
1777    //**********************************************************************************************
1778 
1779    //**Restructuring SMP assignment****************************************************************
1780    /*! \cond BLAZE_INTERNAL */
1781    /*!\brief Restructuring SMP assignment of a transpose dense matrix-sparse matrix multiplication
1782    //        (\f$ C=A*B \f$).
1783    // \ingroup dense_matrix
1784    //
1785    // \param lhs The target left-hand side matrix.
1786    // \param rhs The right-hand side multiplication expression to be assigned.
1787    // \return void
1788    //
1789    // This function implements the symmetry-based restructuring SMP assignment of a transpose
1790    // dense matrix-sparse matrix multiplication expression. Due to the explicit application
1791    // of the SFINAE principle this function can only be selected by the compiler in case the
1792    // symmetry of either of the two matrix operands can be exploited.
1793    */
1794    template< typename MT  // Type of the target matrix
1795            , bool SO >    // Storage order of the target matrix
1796    friend inline auto smpAssign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1797       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1798    {
1799       BLAZE_FUNCTION_TRACE;
1800 
1801       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1802       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1803 
1804       const ForwardFunctor fwd;
1805 
1806       smpAssign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
1807    }
1808    /*! \endcond */
1809    //**********************************************************************************************
1810 
1811    //**SMP addition assignment to dense matrices***************************************************
1812    /*! \cond BLAZE_INTERNAL */
1813    /*!\brief SMP addition assignment of a transpose dense matrix-sparse matrix multiplication to
1814    //        a dense matrix (\f$ C+=A*B \f$).
1815    // \ingroup dense_matrix
1816    //
1817    // \param lhs The target left-hand side dense matrix.
1818    // \param rhs The right-hand side multiplication expression to be added.
1819    // \return void
1820    //
1821    // This function implements the performance optimized addition assignment of a transpose dense
1822    // matrix-sparse matrix multiplication expression to a dense matrix. Due to the explicit
1823    // application of the SFINAE principle this function can only be selected by the compiler in
1824    // case either of the two matrix operands requires an intermediate evaluation and no symmetry
1825    // can be exploited.
1826    */
1827    template< typename MT  // Type of the target dense matrix
1828            , bool SO >    // Storage order of the target sparse matrix
1829    friend inline auto smpAddAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1830       -> EnableIf_t< IsEvaluationRequired_v<MT,MT1,MT2> >
1831    {
1832       BLAZE_FUNCTION_TRACE;
1833 
1834       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1835       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1836 
1837       LT A( rhs.lhs_ );  // Evaluation of the left-hand side dense matrix operand
1838       RT B( rhs.rhs_ );  // Evaluation of the right-hand side sparse matrix operand
1839 
1840       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
1841       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
1842       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
1843       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
1844       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
1845       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
1846 
1847       smpAddAssign( *lhs, A * B );
1848    }
1849    /*! \endcond */
1850    //**********************************************************************************************
1851 
1852    //**Restructuring SMP addition assignment*******************************************************
1853    /*! \cond BLAZE_INTERNAL */
1854    /*!\brief Restructuring SMP addition assignment of a transpose dense matrix-sparse matrix
1855    //        multiplication (\f$ C+=A*B \f$).
1856    // \ingroup dense_matrix
1857    //
1858    // \param lhs The target left-hand side matrix.
1859    // \param rhs The right-hand side multiplication expression to be added.
1860    // \return void
1861    //
1862    // This function implements the symmetry-based restructuring SMP addition assignment of
1863    // a transpose dense matrix-sparse matrix multiplication expression. Due to the explicit
1864    // application of the SFINAE principle this function can only be selected by the compiler
1865    // in case the symmetry of either of the two matrix operands can be exploited.
1866    */
1867    template< typename MT  // Type of the target matrix
1868            , bool SO >    // Storage order of the target matrix
1869    friend inline auto smpAddAssign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1870       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1871    {
1872       BLAZE_FUNCTION_TRACE;
1873 
1874       BLAZE_CONSTRAINT_MUST_NOT_BE_SYMMETRIC_MATRIX_TYPE( MT );
1875 
1876       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1877       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1878 
1879       const ForwardFunctor fwd;
1880 
1881       smpAddAssign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
1882    }
1883    /*! \endcond */
1884    //**********************************************************************************************
1885 
1886    //**SMP addition assignment to sparse matrices**************************************************
1887    // No special implementation for the SMP addition assignment to sparse matrices.
1888    //**********************************************************************************************
1889 
1890    //**SMP subtraction assignment to dense matrices************************************************
1891    /*! \cond BLAZE_INTERNAL */
1892    /*!\brief SMP subtraction assignment of a transpose dense matrix-sparse matrix multiplication
1893    //        to a dense matrix (\f$ C-=A*B \f$).
1894    // \ingroup dense_matrix
1895    //
1896    // \param lhs The target left-hand side dense matrix.
1897    // \param rhs The right-hand side multiplication expression to be subtracted.
1898    // \return void
1899    //
1900    // This function implements the performance optimized subtraction assignment of a transpose
1901    // dense matrix-sparse matrix multiplication expression to a dense matrix. Due to the explicit
1902    // application of the SFINAE principle this function can only be selected by the compiler in
1903    // case either of the two matrix operands requires an intermediate evaluation and no symmetry
1904    // can be exploited.
1905    */
1906    template< typename MT  // Type of the target dense matrix
1907            , bool SO >    // Storage order of the target sparse matrix
1908    friend inline auto smpSubAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1909       -> EnableIf_t< IsEvaluationRequired_v<MT,MT1,MT2> >
1910    {
1911       BLAZE_FUNCTION_TRACE;
1912 
1913       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1914       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1915 
1916       LT A( rhs.lhs_ );  // Evaluation of the left-hand side dense matrix operand
1917       RT B( rhs.rhs_ );  // Evaluation of the right-hand side sparse matrix operand
1918 
1919       BLAZE_INTERNAL_ASSERT( A.rows()    == rhs.lhs_.rows()   , "Invalid number of rows"    );
1920       BLAZE_INTERNAL_ASSERT( A.columns() == rhs.lhs_.columns(), "Invalid number of columns" );
1921       BLAZE_INTERNAL_ASSERT( B.rows()    == rhs.rhs_.rows()   , "Invalid number of rows"    );
1922       BLAZE_INTERNAL_ASSERT( B.columns() == rhs.rhs_.columns(), "Invalid number of columns" );
1923       BLAZE_INTERNAL_ASSERT( A.rows()    == (*lhs).rows()     , "Invalid number of rows"    );
1924       BLAZE_INTERNAL_ASSERT( B.columns() == (*lhs).columns()  , "Invalid number of columns" );
1925 
1926       smpSubAssign( *lhs, A * B );
1927    }
1928    /*! \endcond */
1929    //**********************************************************************************************
1930 
1931    //**Restructuring SMP subtraction assignment****************************************************
1932    /*! \cond BLAZE_INTERNAL */
1933    /*!\brief Restructuring SMP subtraction assignment of a transpose dense matrix-sparse matrix
1934    //        multiplication (\f$ C-=A*B \f$).
1935    // \ingroup dense_matrix
1936    //
1937    // \param lhs The target left-hand side matrix.
1938    // \param rhs The right-hand side multiplication expression to be subtracted.
1939    // \return void
1940    //
1941    // This function implements the symmetry-based restructuring SMP subtraction assignment of
1942    // a transpose dense matrix-sparse matrix multiplication expression. Due to the explicit
1943    // application of the SFINAE principle this function can only be selected by the compiler
1944    // in case the symmetry of either of the two matrix operands can be exploited.
1945    */
1946    template< typename MT  // Type of the target matrix
1947            , bool SO >    // Storage order of the target matrix
1948    friend inline auto smpSubAssign( Matrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1949       -> EnableIf_t< CanExploitSymmetry_v<MT,MT1,MT2> >
1950    {
1951       BLAZE_FUNCTION_TRACE;
1952 
1953       BLAZE_CONSTRAINT_MUST_NOT_BE_SYMMETRIC_MATRIX_TYPE( MT );
1954 
1955       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1956       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1957 
1958       const ForwardFunctor fwd;
1959 
1960       smpSubAssign( *lhs, fwd( rhs.lhs_ * trans( rhs.rhs_ ) ) );
1961    }
1962    /*! \endcond */
1963    //**********************************************************************************************
1964 
1965    //**SMP subtraction assignment to sparse matrices***********************************************
1966    // No special implementation for the SMP subtraction assignment to sparse matrices.
1967    //**********************************************************************************************
1968 
1969    //**SMP Schur product assignment to dense matrices**********************************************
1970    /*! \cond BLAZE_INTERNAL */
1971    /*!\brief SMP Schur product assignment of a transpose dense matrix-sparse matrix multiplication
1972    //        to a dense matrix (\f$ C\circ=A*B \f$).
1973    // \ingroup dense_matrix
1974    //
1975    // \param lhs The target left-hand side dense matrix.
1976    // \param rhs The right-hand side multiplication expression for the Schur product.
1977    // \return void
1978    //
1979    // This function implements the performance optimized Schur product assignment of a transpose
1980    // dense matrix-sparse matrix multiplication expression to a dense matrix.
1981    */
1982    template< typename MT  // Type of the target dense matrix
1983            , bool SO >    // Storage order of the target sparse matrix
smpSchurAssign(DenseMatrix<MT,SO> & lhs,const TDMatSMatMultExpr & rhs)1984    friend inline void smpSchurAssign( DenseMatrix<MT,SO>& lhs, const TDMatSMatMultExpr& rhs )
1985    {
1986       BLAZE_FUNCTION_TRACE;
1987 
1988       BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( ResultType );
1989       BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( ResultType );
1990       BLAZE_CONSTRAINT_MUST_NOT_REQUIRE_EVALUATION( ResultType );
1991 
1992       BLAZE_INTERNAL_ASSERT( (*lhs).rows()    == rhs.rows()   , "Invalid number of rows"    );
1993       BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
1994 
1995       const ResultType tmp( rhs );
1996       smpSchurAssign( *lhs, tmp );
1997    }
1998    /*! \endcond */
1999    //**********************************************************************************************
2000 
2001    //**SMP Schur product assignment to sparse matrices*********************************************
2002    // No special implementation for the SMP Schur product assignment to sparse matrices.
2003    //**********************************************************************************************
2004 
2005    //**SMP multiplication assignment to dense matrices*********************************************
2006    // No special implementation for the SMP multiplication assignment to dense matrices.
2007    //**********************************************************************************************
2008 
2009    //**SMP multiplication assignment to sparse matrices********************************************
2010    // No special implementation for the SMP multiplication assignment to sparse matrices.
2011    //**********************************************************************************************
2012 
2013    //**Compile time checks*************************************************************************
2014    /*! \cond BLAZE_INTERNAL */
2015    BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( MT1 );
2016    BLAZE_CONSTRAINT_MUST_BE_COLUMN_MAJOR_MATRIX_TYPE( MT1 );
2017    BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT2 );
2018    BLAZE_CONSTRAINT_MUST_BE_ROW_MAJOR_MATRIX_TYPE( MT2 );
2019    BLAZE_CONSTRAINT_MUST_NOT_BE_ZERO_TYPE( MT2 );
2020    BLAZE_CONSTRAINT_MUST_FORM_VALID_MATMATMULTEXPR( MT1, MT2 );
2021    /*! \endcond */
2022    //**********************************************************************************************
2023 };
2024 //*************************************************************************************************
2025 
2026 
2027 
2028 
2029 //=================================================================================================
2030 //
2031 //  GLOBAL BINARY ARITHMETIC OPERATORS
2032 //
2033 //=================================================================================================
2034 
2035 //*************************************************************************************************
2036 /*! \cond BLAZE_INTERNAL */
2037 /*!\brief Backend implementation of the multiplication between a column-major dense matrix and a
2038 //        row-major sparse matrix (\f$ A=B*C \f$).
2039 // \ingroup dense_matrix
2040 //
2041 // \param lhs The left-hand side dense matrix for the multiplication.
2042 // \param rhs The right-hand side sparse matrix for the multiplication.
2043 // \return The product of the two matrices.
2044 //
2045 // This function implements a performance optimized treatment of the multiplication between a
2046 // column-major dense matrix and a row-major sparse matrix.
2047 */
2048 template< typename MT1  // Type of the left-hand side dense matrix
2049         , typename MT2  // Type of the right-hand side sparse matrix
2050         , DisableIf_t< ( IsIdentity_v<MT2> &&
2051                          IsSame_v< ElementType_t<MT1>, ElementType_t<MT2> > ) ||
2052                        IsZero_v<MT2> >* = nullptr >
2053 inline const TDMatSMatMultExpr<MT1,MT2,false,false,false,false>
tdmatsmatmult(const DenseMatrix<MT1,true> & lhs,const SparseMatrix<MT2,false> & rhs)2054    tdmatsmatmult( const DenseMatrix<MT1,true>& lhs, const SparseMatrix<MT2,false>& rhs )
2055 {
2056    BLAZE_FUNCTION_TRACE;
2057 
2058    BLAZE_INTERNAL_ASSERT( (*lhs).columns() == (*rhs).rows(), "Invalid matrix sizes" );
2059 
2060    return TDMatSMatMultExpr<MT1,MT2,false,false,false,false>( *lhs, *rhs );
2061 }
2062 /*! \endcond */
2063 //*************************************************************************************************
2064 
2065 
2066 //*************************************************************************************************
2067 /*! \cond BLAZE_INTERNAL */
2068 /*!\brief Backend implementation of the multiplication between a column-major dense matrix and a
2069 //        row-major identity matrix (\f$ A=B*C \f$).
2070 // \ingroup dense_matrix
2071 //
2072 // \param lhs The left-hand side dense matrix for the multiplication.
2073 // \param rhs The right-hand side identity matrix for the multiplication.
2074 // \return Reference to the left-hand side dense matrix.
2075 //
2076 // This function implements a performance optimized treatment of the multiplication between
2077 // a column-major dense matrix and a row-major identity matrix. It returns a reference to the
2078 // left-hand side dense matrix.
2079 */
2080 template< typename MT1  // Type of the left-hand side dense matrix
2081         , typename MT2  // Type of the right-hand side sparse matrix
2082         , EnableIf_t< IsIdentity_v<MT2> &&
2083                       IsSame_v< ElementType_t<MT1>, ElementType_t<MT2> > >* = nullptr >
2084 inline const MT1&
tdmatsmatmult(const DenseMatrix<MT1,true> & lhs,const SparseMatrix<MT2,false> & rhs)2085    tdmatsmatmult( const DenseMatrix<MT1,true>& lhs, const SparseMatrix<MT2,false>& rhs )
2086 {
2087    BLAZE_FUNCTION_TRACE;
2088 
2089    MAYBE_UNUSED( rhs );
2090 
2091    BLAZE_INTERNAL_ASSERT( (*lhs).columns() == (*rhs).rows(), "Invalid matrix sizes" );
2092 
2093    return (*lhs);
2094 }
2095 /*! \endcond */
2096 //*************************************************************************************************
2097 
2098 
2099 //*************************************************************************************************
2100 /*! \cond BLAZE_INTERNAL */
2101 /*!\brief Backend implementation of the multiplication between a column-major dense matrix and a
2102 //        row-major zero matrix (\f$ A=B*C \f$).
2103 // \ingroup dense_matrix
2104 //
2105 // \param lhs The left-hand side dense matrix for the multiplication.
2106 // \param rhs The right-hand side zero matrix for the multiplication.
2107 // \return The resulting zero matrix.
2108 //
2109 // This function implements a performance optimized treatment of the multiplication between a
2110 // column-major dense matrix and a row-major zero matrix. It returns a zero matrix.
2111 */
2112 template< typename MT1  // Type of the left-hand side dense matrix
2113         , typename MT2  // Type of the right-hand side sparse matrix
2114         , EnableIf_t< IsZero_v<MT2> >* = nullptr >
decltype(auto)2115 inline decltype(auto)
2116    tdmatsmatmult( const DenseMatrix<MT1,true>& lhs, const SparseMatrix<MT2,false>& rhs )
2117 {
2118    BLAZE_FUNCTION_TRACE;
2119 
2120    BLAZE_INTERNAL_ASSERT( (*lhs).columns() == (*rhs).rows(), "Invalid matrix sizes" );
2121 
2122    using ReturnType = const MultTrait_t< ResultType_t<MT1>, ResultType_t<MT2> >;
2123 
2124    BLAZE_CONSTRAINT_MUST_BE_ROW_MAJOR_MATRIX_TYPE( ReturnType );
2125    BLAZE_CONSTRAINT_MUST_BE_ZERO_TYPE( ReturnType );
2126 
2127    return ReturnType( (*lhs).rows(), (*rhs).columns() );
2128 }
2129 /*! \endcond */
2130 //*************************************************************************************************
2131 
2132 
2133 //*************************************************************************************************
2134 /*!\brief Multiplication operator for the multiplication of a column-major dense matrix and a
2135 //        row-major sparse matrix (\f$ A=B*C \f$).
2136 // \ingroup dense_matrix
2137 //
2138 // \param lhs The left-hand side dense matrix for the multiplication.
2139 // \param rhs The right-hand side sparse matrix for the multiplication.
2140 // \return The resulting matrix.
2141 // \exception std::invalid_argument Matrix sizes do not match.
2142 //
2143 // This operator represents the multiplication of a column-major dense matrix and a row-major
2144 // sparse matrix:
2145 
2146    \code
2147    using blaze::rowMajor;
2148    using blaze::columnMajor;
2149 
2150    blaze::DynamicMatrix<double,columnMajor> A, C;
2151    blaze::CompressedMatrix<double,rowMajor> B;
2152    // ... Resizing and initialization
2153    C = A * B;
2154    \endcode
2155 
2156 // The operator returns an expression representing a dense matrix of the higher-order element
2157 // type of the two involved matrix element types \a MT1::ElementType and \a MT2::ElementType.
2158 // Both matrix types \a MT1 and \a MT2 as well as the two element types \a MT1::ElementType
2159 // and \a MT2::ElementType have to be supported by the MultTrait class template.\n
2160 // In case the current sizes of the two given matrices don't match, a \a std::invalid_argument
2161 // is thrown.
2162 */
2163 template< typename MT1    // Type of the left-hand side dense matrix
2164         , typename MT2 >  // Type of the right-hand side sparse matrix
decltype(auto)2165 inline decltype(auto)
2166    operator*( const DenseMatrix<MT1,true>& lhs, const SparseMatrix<MT2,false>& rhs )
2167 {
2168    BLAZE_FUNCTION_TRACE;
2169 
2170    if( (*lhs).columns() != (*rhs).rows() ) {
2171       BLAZE_THROW_INVALID_ARGUMENT( "Matrix sizes do not match" );
2172    }
2173 
2174    return tdmatsmatmult( *lhs, *rhs );
2175 }
2176 //*************************************************************************************************
2177 
2178 
2179 
2180 
2181 //=================================================================================================
2182 //
2183 //  GLOBAL FUNCTIONS
2184 //
2185 //=================================================================================================
2186 
2187 //*************************************************************************************************
2188 /*! \cond BLAZE_INTERNAL */
2189 /*!\brief Declares the given non-symmetric matrix multiplication expression as symmetric.
2190 // \ingroup dense_matrix
2191 //
2192 // \param dm The input matrix multiplication expression.
2193 // \return The redeclared dense matrix multiplication expression.
2194 // \exception std::invalid_argument Invalid symmetric matrix specification.
2195 //
2196 // The \a declsym function declares the given non-symmetric matrix multiplication expression
2197 // \a dm as symmetric. The function returns an expression representing the operation. In case
2198 // the given expression does not represent a square matrix, a \a std::invalid_argument exception
2199 // is thrown.\n
2200 // The following example demonstrates the use of the \a declsym function:
2201 
2202    \code
2203    using blaze::rowMajor;
2204    using blaze::columnMajor;
2205 
2206    blaze::DynamicMatrix<double,columnMajor> A, C;
2207    blaze::CompressedMatrix<double,rowMajor> B;
2208    // ... Resizing and initialization
2209    C = declsym( A * B );
2210    \endcode
2211 */
2212 template< typename MT1  // Type of the left-hand side dense matrix
2213         , typename MT2  // Type of the right-hand side dense matrix
2214         , bool SF       // Symmetry flag
2215         , bool HF       // Hermitian flag
2216         , bool LF       // Lower flag
2217         , bool UF >     // Upper flag
decltype(auto)2218 inline decltype(auto) declsym( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>& dm )
2219 {
2220    BLAZE_FUNCTION_TRACE;
2221 
2222    if( !isSquare( dm ) ) {
2223       BLAZE_THROW_INVALID_ARGUMENT( "Invalid symmetric matrix specification" );
2224    }
2225 
2226    using ReturnType = const TDMatSMatMultExpr<MT1,MT2,true,HF,LF,UF>;
2227    return ReturnType( dm.leftOperand(), dm.rightOperand() );
2228 }
2229 /*! \endcond */
2230 //*************************************************************************************************
2231 
2232 
2233 //*************************************************************************************************
2234 /*! \cond BLAZE_INTERNAL */
2235 /*!\brief Declares the given non-Hermitian matrix multiplication expression as Hermitian.
2236 // \ingroup dense_matrix
2237 //
2238 // \param dm The input matrix multiplication expression.
2239 // \return The redeclared dense matrix multiplication expression.
2240 // \exception std::invalid_argument Invalid Hermitian matrix specification.
2241 //
2242 // The \a declherm function declares the given non-Hermitian matrix multiplication expression
2243 // \a dm as Hermitian. The function returns an expression representing the operation. In case
2244 // the given expression does not represent a square matrix, a \a std::invalid_argument exception
2245 // is thrown.\n
2246 // The following example demonstrates the use of the \a declherm function:
2247 
2248    \code
2249    using blaze::rowMajor;
2250    using blaze::columnMajor;
2251 
2252    blaze::DynamicMatrix<double,columnMajor> A, C;
2253    blaze::CompressedMatrix<double,rowMajor> B;
2254    // ... Resizing and initialization
2255    C = declherm( A * B );
2256    \endcode
2257 */
2258 template< typename MT1  // Type of the left-hand side dense matrix
2259         , typename MT2  // Type of the right-hand side dense matrix
2260         , bool SF       // Symmetry flag
2261         , bool HF       // Hermitian flag
2262         , bool LF       // Lower flag
2263         , bool UF >     // Upper flag
decltype(auto)2264 inline decltype(auto) declherm( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>& dm )
2265 {
2266    BLAZE_FUNCTION_TRACE;
2267 
2268    if( !isSquare( dm ) ) {
2269       BLAZE_THROW_INVALID_ARGUMENT( "Invalid Hermitian matrix specification" );
2270    }
2271 
2272    using ReturnType = const TDMatSMatMultExpr<MT1,MT2,SF,true,LF,UF>;
2273    return ReturnType( dm.leftOperand(), dm.rightOperand() );
2274 }
2275 /*! \endcond */
2276 //*************************************************************************************************
2277 
2278 
2279 //*************************************************************************************************
2280 /*! \cond BLAZE_INTERNAL */
2281 /*!\brief Declares the given non-lower matrix multiplication expression as lower.
2282 // \ingroup dense_matrix
2283 //
2284 // \param dm The input matrix multiplication expression.
2285 // \return The redeclared dense matrix multiplication expression.
2286 // \exception std::invalid_argument Invalid lower matrix specification.
2287 //
2288 // The \a decllow function declares the given non-lower matrix multiplication expression
2289 // \a dm as lower. The function returns an expression representing the operation. In case
2290 // the given expression does not represent a square matrix, a \a std::invalid_argument
2291 // exception is thrown.\n
2292 // The following example demonstrates the use of the \a decllow function:
2293 
2294    \code
2295    using blaze::rowMajor;
2296    using blaze::columnMajor;
2297 
2298    blaze::DynamicMatrix<double,columnMajor> A, C;
2299    blaze::CompressedMatrix<double,rowMajor> B;
2300    // ... Resizing and initialization
2301    C = decllow( A * B );
2302    \endcode
2303 */
2304 template< typename MT1  // Type of the left-hand side dense matrix
2305         , typename MT2  // Type of the right-hand side dense matrix
2306         , bool SF       // Symmetry flag
2307         , bool HF       // Hermitian flag
2308         , bool LF       // Lower flag
2309         , bool UF >     // Upper flag
decltype(auto)2310 inline decltype(auto) decllow( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>& dm )
2311 {
2312    BLAZE_FUNCTION_TRACE;
2313 
2314    if( !isSquare( dm ) ) {
2315       BLAZE_THROW_INVALID_ARGUMENT( "Invalid lower matrix specification" );
2316    }
2317 
2318    using ReturnType = const TDMatSMatMultExpr<MT1,MT2,SF,HF,true,UF>;
2319    return ReturnType( dm.leftOperand(), dm.rightOperand() );
2320 }
2321 /*! \endcond */
2322 //*************************************************************************************************
2323 
2324 
2325 //*************************************************************************************************
2326 /*! \cond BLAZE_INTERNAL */
2327 /*!\brief Declares the given non-unilower matrix multiplication expression as unilower.
2328 // \ingroup dense_matrix
2329 //
2330 // \param dm The input matrix multiplication expression.
2331 // \return The redeclared dense matrix multiplication expression.
2332 // \exception std::invalid_argument Invalid unilower matrix specification.
2333 //
2334 // The \a declunilow function declares the given non-unilower matrix multiplication expression
2335 // \a dm as unilower. The function returns an expression representing the operation. In case
2336 // the given expression does not represent a square matrix, a \a std::invalid_argument
2337 // exception is thrown.\n
2338 // The following example demonstrates the use of the \a declunilow function:
2339 
2340    \code
2341    using blaze::rowMajor;
2342    using blaze::columnMajor;
2343 
2344    blaze::DynamicMatrix<double,columnMajor> A, C;
2345    blaze::CompressedMatrix<double,rowMajor> B;
2346    // ... Resizing and initialization
2347    C = declunilow( A * B );
2348    \endcode
2349 */
2350 template< typename MT1  // Type of the left-hand side dense matrix
2351         , typename MT2  // Type of the right-hand side dense matrix
2352         , bool SF       // Symmetry flag
2353         , bool HF       // Hermitian flag
2354         , bool UF >     // Upper flag
decltype(auto)2355 inline decltype(auto) declunilow( const TDMatSMatMultExpr<MT1,MT2,SF,HF,false,UF>& dm )
2356 {
2357    BLAZE_FUNCTION_TRACE;
2358 
2359    if( !isSquare( dm ) ) {
2360       BLAZE_THROW_INVALID_ARGUMENT( "Invalid unilower matrix specification" );
2361    }
2362 
2363    return declunilow( decllow( *dm ) );
2364 }
2365 /*! \endcond */
2366 //*************************************************************************************************
2367 
2368 
2369 //*************************************************************************************************
2370 /*! \cond BLAZE_INTERNAL */
2371 /*!\brief Declares the given non-strictly-lower matrix multiplication expression as strictly lower.
2372 // \ingroup dense_matrix
2373 //
2374 // \param dm The input matrix multiplication expression.
2375 // \return The redeclared dense matrix multiplication expression.
2376 // \exception std::invalid_argument Invalid strlower matrix specification.
2377 //
2378 // The \a declstrlow function declares the given non-strictly-lower matrix multiplication
2379 // expression \a dm as strictly lower. The function returns an expression representing the
2380 // operation. In case the given expression does not represent a square matrix, a
2381 // \a std::invalid_argument exception is thrown.\n
2382 // The following example demonstrates the use of the \a declstrlow function:
2383 
2384    \code
2385    using blaze::rowMajor;
2386    using blaze::columnMajor;
2387 
2388    blaze::DynamicMatrix<double,columnMajor> A, C;
2389    blaze::CompressedMatrix<double,rowMajor> B;
2390    // ... Resizing and initialization
2391    C = declstrlow( A * B );
2392    \endcode
2393 */
2394 template< typename MT1  // Type of the left-hand side dense matrix
2395         , typename MT2  // Type of the right-hand side dense matrix
2396         , bool SF       // Symmetry flag
2397         , bool HF       // Hermitian flag
2398         , bool UF >     // Upper flag
decltype(auto)2399 inline decltype(auto) declstrlow( const TDMatSMatMultExpr<MT1,MT2,SF,HF,false,UF>& dm )
2400 {
2401    BLAZE_FUNCTION_TRACE;
2402 
2403    if( !isSquare( dm ) ) {
2404       BLAZE_THROW_INVALID_ARGUMENT( "Invalid strictly lower matrix specification" );
2405    }
2406 
2407    return declstrlow( decllow( *dm ) );
2408 }
2409 /*! \endcond */
2410 //*************************************************************************************************
2411 
2412 
2413 //*************************************************************************************************
2414 /*! \cond BLAZE_INTERNAL */
2415 /*!\brief Declares the given non-upper matrix multiplication expression as upper.
2416 // \ingroup dense_matrix
2417 //
2418 // \param dm The input matrix multiplication expression.
2419 // \return The redeclared dense matrix multiplication expression.
2420 // \exception std::invalid_argument Invalid upper matrix specification.
2421 //
2422 // The \a declupp function declares the given non-upper matrix multiplication expression
2423 // \a dm as upper. The function returns an expression representing the operation. In case
2424 // the given expression does not represent a square matrix, a \a std::invalid_argument
2425 // exception is thrown.\n
2426 // The following example demonstrates the use of the \a declupp function:
2427 
2428    \code
2429    using blaze::rowMajor;
2430    using blaze::columnMajor;
2431 
2432    blaze::DynamicMatrix<double,columnMajor> A, C;
2433    blaze::CompressedMatrix<double,rowMajor> B;
2434    // ... Resizing and initialization
2435    C = declupp( A * B );
2436    \endcode
2437 */
2438 template< typename MT1  // Type of the left-hand side dense matrix
2439         , typename MT2  // Type of the right-hand side dense matrix
2440         , bool SF       // Symmetry flag
2441         , bool HF       // Hermitian flag
2442         , bool LF       // Lower flag
2443         , bool UF >     // Upper flag
decltype(auto)2444 inline decltype(auto) declupp( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>& dm )
2445 {
2446    BLAZE_FUNCTION_TRACE;
2447 
2448    if( !isSquare( dm ) ) {
2449       BLAZE_THROW_INVALID_ARGUMENT( "Invalid upper matrix specification" );
2450    }
2451 
2452    using ReturnType = const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,true>;
2453    return ReturnType( dm.leftOperand(), dm.rightOperand() );
2454 }
2455 /*! \endcond */
2456 //*************************************************************************************************
2457 
2458 
2459 //*************************************************************************************************
2460 /*! \cond BLAZE_INTERNAL */
2461 /*!\brief Declares the given non-uniupper matrix multiplication expression as uniupper.
2462 // \ingroup dense_matrix
2463 //
2464 // \param dm The input matrix multiplication expression.
2465 // \return The redeclared dense matrix multiplication expression.
2466 // \exception std::invalid_argument Invalid uniupper matrix specification.
2467 //
2468 // The \a decluniupp function declares the given non-uniupper matrix multiplication expression
2469 // \a dm as uniupper. The function returns an expression representing the operation. In case
2470 // the given expression does not represent a square matrix, a \a std::invalid_argument
2471 // exception is thrown.\n
2472 // The following example demonstrates the use of the \a decluniupp function:
2473 
2474    \code
2475    using blaze::rowMajor;
2476    using blaze::columnMajor;
2477 
2478    blaze::DynamicMatrix<double,columnMajor> A, C;
2479    blaze::CompressedMatrix<double,rowMajor> B;
2480    // ... Resizing and initialization
2481    C = decluniupp( A * B );
2482    \endcode
2483 */
2484 template< typename MT1  // Type of the left-hand side dense matrix
2485         , typename MT2  // Type of the right-hand side dense matrix
2486         , bool SF       // Symmetry flag
2487         , bool HF       // Hermitian flag
2488         , bool LF >     // Lower flag
decltype(auto)2489 inline decltype(auto) decluniupp( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,false>& dm )
2490 {
2491    BLAZE_FUNCTION_TRACE;
2492 
2493    if( !isSquare( dm ) ) {
2494       BLAZE_THROW_INVALID_ARGUMENT( "Invalid uniupper matrix specification" );
2495    }
2496 
2497    return decluniupp( declupp( *dm ) );
2498 }
2499 /*! \endcond */
2500 //*************************************************************************************************
2501 
2502 
2503 //*************************************************************************************************
2504 /*! \cond BLAZE_INTERNAL */
2505 /*!\brief Declares the given non-strictly-upper matrix multiplication expression as strictly upper.
2506 // \ingroup dense_matrix
2507 //
2508 // \param dm The input matrix multiplication expression.
2509 // \return The redeclared dense matrix multiplication expression.
2510 // \exception std::invalid_argument Invalid strupper matrix specification.
2511 //
2512 // The \a declstrupp function declares the given non-strictly-upper matrix multiplication
2513 // expression \a dm as strictly upper. The function returns an expression representing the
2514 // operation. In case the given expression does not represent a square matrix, a
2515 // \a std::invalid_argument exception is thrown.\n
2516 // The folupping example demonstrates the use of the \a declstrupp function:
2517 
2518    \code
2519    using blaze::rowMajor;
2520    using blaze::columnMajor;
2521 
2522    blaze::DynamicMatrix<double,columnMajor> A, C;
2523    blaze::CompressedMatrix<double,rowMajor> B;
2524    // ... Resizing and initialization
2525    C = declstrupp( A * B );
2526    \endcode
2527 */
2528 template< typename MT1  // Type of the left-hand side dense matrix
2529         , typename MT2  // Type of the right-hand side dense matrix
2530         , bool SF       // Symmetry flag
2531         , bool HF       // Hermitian flag
2532         , bool LF >     // Lower flag
decltype(auto)2533 inline decltype(auto) declstrupp( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,false>& dm )
2534 {
2535    BLAZE_FUNCTION_TRACE;
2536 
2537    if( !isSquare( dm ) ) {
2538       BLAZE_THROW_INVALID_ARGUMENT( "Invalid strictly upper matrix specification" );
2539    }
2540 
2541    return declstrupp( declupp( *dm ) );
2542 }
2543 /*! \endcond */
2544 //*************************************************************************************************
2545 
2546 
2547 //*************************************************************************************************
2548 /*! \cond BLAZE_INTERNAL */
2549 /*!\brief Declares the given non-diagonal matrix multiplication expression as diagonal.
2550 // \ingroup dense_matrix
2551 //
2552 // \param dm The input matrix multiplication expression.
2553 // \return The redeclared dense matrix multiplication expression.
2554 // \exception std::invalid_argument Invalid diagonal matrix specification.
2555 //
2556 // The \a decldiag function declares the given non-diagonal matrix multiplication expression
2557 // \a dm as diagonal. The function returns an expression representing the operation. In case
2558 // the given expression does not represent a square matrix, a \a std::invalid_argument exception
2559 // is thrown.\n
2560 // The following example demonstrates the use of the \a decldiag function:
2561 
2562    \code
2563    using blaze::rowMajor;
2564    using blaze::columnMajor;
2565 
2566    blaze::DynamicMatrix<double,columnMajor> A, C;
2567    blaze::CompressedMatrix<double,rowMajor> B;
2568    // ... Resizing and initialization
2569    C = decldiag( A * B );
2570    \endcode
2571 */
2572 template< typename MT1  // Type of the left-hand side dense matrix
2573         , typename MT2  // Type of the right-hand side dense matrix
2574         , bool SF       // Symmetry flag
2575         , bool HF       // Hermitian flag
2576         , bool LF       // Lower flag
2577         , bool UF >     // Upper flag
decltype(auto)2578 inline decltype(auto) decldiag( const TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>& dm )
2579 {
2580    BLAZE_FUNCTION_TRACE;
2581 
2582    if( !isSquare( dm ) ) {
2583       BLAZE_THROW_INVALID_ARGUMENT( "Invalid diagonal matrix specification" );
2584    }
2585 
2586    using ReturnType = const TDMatSMatMultExpr<MT1,MT2,SF,HF,true,true>;
2587    return ReturnType( dm.leftOperand(), dm.rightOperand() );
2588 }
2589 /*! \endcond */
2590 //*************************************************************************************************
2591 
2592 
2593 
2594 
2595 //=================================================================================================
2596 //
2597 //  SIZE SPECIALIZATIONS
2598 //
2599 //=================================================================================================
2600 
2601 //*************************************************************************************************
2602 /*! \cond BLAZE_INTERNAL */
2603 template< typename MT1, typename MT2, bool SF, bool HF, bool LF, bool UF >
2604 struct Size< TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>, 0UL >
2605    : public Size<MT1,0UL>
2606 {};
2607 
2608 template< typename MT1, typename MT2, bool SF, bool HF, bool LF, bool UF >
2609 struct Size< TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF>, 1UL >
2610    : public Size<MT2,1UL>
2611 {};
2612 /*! \endcond */
2613 //*************************************************************************************************
2614 
2615 
2616 
2617 
2618 //=================================================================================================
2619 //
2620 //  ISALIGNED SPECIALIZATIONS
2621 //
2622 //=================================================================================================
2623 
2624 //*************************************************************************************************
2625 /*! \cond BLAZE_INTERNAL */
2626 template< typename MT1, typename MT2, bool SF, bool HF, bool LF, bool UF >
2627 struct IsAligned< TDMatSMatMultExpr<MT1,MT2,SF,HF,LF,UF> >
2628    : public IsAligned<MT1>
2629 {};
2630 /*! \endcond */
2631 //*************************************************************************************************
2632 
2633 } // namespace blaze
2634 
2635 #endif
2636