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