1 //=================================================================================================
2 /*!
3 // \file blaze/math/expressions/DMatDMatKronExpr.h
4 // \brief Header file for the dense matrix/dense matrix Kronecker product 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_DMATDMATKRONEXPR_H_
36 #define _BLAZE_MATH_EXPRESSIONS_DMATDMATKRONEXPR_H_
37
38
39 //*************************************************************************************************
40 // Includes
41 //*************************************************************************************************
42
43 #include <utility>
44 #include <blaze/math/Aliases.h>
45 #include <blaze/math/AlignmentFlag.h>
46 #include <blaze/math/constraints/DenseMatrix.h>
47 #include <blaze/math/constraints/MatMatKronExpr.h>
48 #include <blaze/math/constraints/RequiresEvaluation.h>
49 #include <blaze/math/constraints/StorageOrder.h>
50 #include <blaze/math/Exception.h>
51 #include <blaze/math/expressions/Computation.h>
52 #include <blaze/math/expressions/DenseMatrix.h>
53 #include <blaze/math/expressions/Forward.h>
54 #include <blaze/math/expressions/MatMatKronExpr.h>
55 #include <blaze/math/shims/IsDefault.h>
56 #include <blaze/math/shims/Reset.h>
57 #include <blaze/math/shims/Serial.h>
58 #include <blaze/math/traits/KronTrait.h>
59 #include <blaze/math/typetraits/IsExpression.h>
60 #include <blaze/math/typetraits/IsLower.h>
61 #include <blaze/math/typetraits/IsStrictlyLower.h>
62 #include <blaze/math/typetraits/IsStrictlyUpper.h>
63 #include <blaze/math/typetraits/IsTemporary.h>
64 #include <blaze/math/typetraits/IsUpper.h>
65 #include <blaze/util/Assert.h>
66 #include <blaze/util/FunctionTrace.h>
67 #include <blaze/util/mpl/If.h>
68 #include <blaze/util/Types.h>
69
70
71 namespace blaze {
72
73 //=================================================================================================
74 //
75 // CLASS DMATDMATKRONEXPR
76 //
77 //=================================================================================================
78
79 //*************************************************************************************************
80 /*!\brief Expression object for dense matrix-dense matrix Kronecker products.
81 // \ingroup dense_matrix_expression
82 //
83 // The DMatDMatKronExpr class represents the compile time expression for Kronecker products
84 // between dense matrices of any storage order.
85 */
86 template< typename MT1 // Type of the left-hand side dense matrix
87 , typename MT2 // Type of the right-hand side dense matrix
88 , bool SO > // Storage order
89 class DMatDMatKronExpr
90 : public MatMatKronExpr< DenseMatrix< DMatDMatKronExpr<MT1,MT2,SO>, SO > >
91 , private Computation
92 {
93 private:
94 //**Type definitions****************************************************************************
95 using RT1 = ResultType_t<MT1>; //!< Result type of the left-hand side dense matrix expression.
96 using RT2 = ResultType_t<MT2>; //!< Result type of the right-hand side dense matrix expression.
97 using RN1 = ReturnType_t<MT1>; //!< Return type of the left-hand side dense matrix expression.
98 using RN2 = ReturnType_t<MT2>; //!< Return type of the right-hand side dense matrix expression.
99 using CT1 = CompositeType_t<MT1>; //!< Composite type of the left-hand side dense matrix expression.
100 using CT2 = CompositeType_t<MT2>; //!< Composite type of the right-hand side dense matrix expression.
101 using ET1 = ElementType_t<MT1>; //!< Element type of the left-hand side dense matrix expression.
102 using ET2 = ElementType_t<MT2>; //!< Element type of the right-hand side dense matrix expression.
103 //**********************************************************************************************
104
105 //**Return type evaluation**********************************************************************
106 //! Compilation switch for the selection of the subscript operator return type.
107 /*! The \a returnExpr compile time constant expression is a compilation switch for the
108 selection of the \a ReturnType. If either matrix operand returns a temporary vector
109 or matrix, \a returnExpr will be set to \a false and the subscript operator will
110 return it's result by value. Otherwise \a returnExpr will be set to \a true and
111 the subscript operator may return it's result as an expression. */
112 static constexpr bool returnExpr = ( !IsTemporary_v<RN1> && !IsTemporary_v<RN2> );
113
114 //! Expression return type for the subscript operator.
115 using ExprReturnType = decltype( std::declval<RN1>() * std::declval<RN2>() );
116 //**********************************************************************************************
117
118 public:
119 //**Type definitions****************************************************************************
120 //! Type of this DMatDMatKronExpr instance.
121 using This = DMatDMatKronExpr<MT1,MT2,SO>;
122
123 //! Base type of this DMatDMatKronExpr instance.
124 using BaseType = MatMatKronExpr< DenseMatrix<This,SO> >;
125
126 using ResultType = KronTrait_t<RT1,RT2>; //!< Result type for expression template evaluations.
127 using OppositeType = OppositeType_t<ResultType>; //!< Result type with opposite storage order for expression template evaluations.
128 using TransposeType = TransposeType_t<ResultType>; //!< Transpose type for expression template evaluations.
129 using ElementType = ElementType_t<ResultType>; //!< Resulting element type.
130
131 //! Return type for expression template evaluations.
132 using ReturnType = const If_t< returnExpr, ExprReturnType, ElementType >;
133
134 //! Data type for composite expression templates.
135 using CompositeType = const ResultType;
136
137 //! Composite type of the left-hand side dense matrix expression.
138 using LeftOperand = If_t< IsExpression_v<MT1>, const MT1, const MT1& >;
139
140 //! Composite type of the right-hand side dense matrix expression.
141 using RightOperand = If_t< IsExpression_v<MT2>, const MT2, const MT2& >;
142 //**********************************************************************************************
143
144 //**Compilation flags***************************************************************************
145 //! Compilation switch for the expression template evaluation strategy.
146 static constexpr bool simdEnabled = false;
147
148 //! Compilation switch for the expression template assignment strategy.
149 static constexpr bool smpAssignable = false;
150 //**********************************************************************************************
151
152 //**Constructor*********************************************************************************
153 /*!\brief Constructor for the DMatDMatKronExpr class.
154 //
155 // \param lhs The left-hand side operand of the Kronecker product expression.
156 // \param rhs The right-hand side operand of the Kronecker product expression.
157 */
DMatDMatKronExpr(const MT1 & lhs,const MT2 & rhs)158 inline DMatDMatKronExpr( const MT1& lhs, const MT2& rhs ) noexcept
159 : lhs_( lhs ) // Left-hand side dense matrix of the Kronecker product expression
160 , rhs_( rhs ) // Right-hand side dense matrix of the Kronecker product expression
161 {}
162 //**********************************************************************************************
163
164 //**Access operator*****************************************************************************
165 /*!\brief 2D-access to the matrix elements.
166 //
167 // \param i Access index for the row. The index has to be in the range \f$[0..M-1]\f$.
168 // \param j Access index for the column. The index has to be in the range \f$[0..N-1]\f$.
169 // \return The resulting value.
170 */
operator()171 inline ReturnType operator()( size_t i, size_t j ) const {
172 BLAZE_INTERNAL_ASSERT( i < rows() , "Invalid row access index" );
173 BLAZE_INTERNAL_ASSERT( j < columns(), "Invalid column access index" );
174 return lhs_( i/rhs_.rows(), j/rhs_.columns() ) * rhs_( i%rhs_.rows(), j%rhs_.columns() );
175 }
176 //**********************************************************************************************
177
178 //**At function*********************************************************************************
179 /*!\brief Checked access to the matrix elements.
180 //
181 // \param i Access index for the row. The index has to be in the range \f$[0..M-1]\f$.
182 // \param j Access index for the column. The index has to be in the range \f$[0..N-1]\f$.
183 // \return The resulting value.
184 // \exception std::out_of_range Invalid matrix access index.
185 */
at(size_t i,size_t j)186 inline ReturnType at( size_t i, size_t j ) const {
187 if( i >= rows() ) {
188 BLAZE_THROW_OUT_OF_RANGE( "Invalid row access index" );
189 }
190 if( j >= columns() ) {
191 BLAZE_THROW_OUT_OF_RANGE( "Invalid column access index" );
192 }
193 return (*this)(i,j);
194 }
195 //**********************************************************************************************
196
197 //**Rows function*******************************************************************************
198 /*!\brief Returns the current number of rows of the matrix.
199 //
200 // \return The number of rows of the matrix.
201 */
rows()202 inline size_t rows() const noexcept {
203 return lhs_.rows() * rhs_.rows();
204 }
205 //**********************************************************************************************
206
207 //**Columns function****************************************************************************
208 /*!\brief Returns the current number of columns of the matrix.
209 //
210 // \return The number of columns of the matrix.
211 */
columns()212 inline size_t columns() const noexcept {
213 return lhs_.columns() * rhs_.columns();
214 }
215 //**********************************************************************************************
216
217 //**Left operand access*************************************************************************
218 /*!\brief Returns the left-hand side dense matrix operand.
219 //
220 // \return The left-hand side dense matrix operand.
221 */
leftOperand()222 inline LeftOperand leftOperand() const noexcept {
223 return lhs_;
224 }
225 //**********************************************************************************************
226
227 //**Right operand access************************************************************************
228 /*!\brief Returns the right-hand side dense matrix operand.
229 //
230 // \return The right-hand side dense matrix operand.
231 */
rightOperand()232 inline RightOperand rightOperand() const noexcept {
233 return rhs_;
234 }
235 //**********************************************************************************************
236
237 //**********************************************************************************************
238 /*!\brief Returns whether the expression can alias with the given address \a alias.
239 //
240 // \param alias The alias to be checked.
241 // \return \a true in case the expression can alias, \a false otherwise.
242 */
243 template< typename T >
canAlias(const T * alias)244 inline bool canAlias( const T* alias ) const noexcept {
245 return ( lhs_.canAlias( alias ) || rhs_.canAlias( alias ) );
246 }
247 //**********************************************************************************************
248
249 //**********************************************************************************************
250 /*!\brief Returns whether the expression is aliased with the given address \a alias.
251 //
252 // \param alias The alias to be checked.
253 // \return \a true in case an alias effect is detected, \a false otherwise.
254 */
255 template< typename T >
isAliased(const T * alias)256 inline bool isAliased( const T* alias ) const noexcept {
257 return ( lhs_.isAliased( alias ) || rhs_.isAliased( alias ) );
258 }
259 //**********************************************************************************************
260
261 //**********************************************************************************************
262 /*!\brief Returns whether the operands of the expression are properly aligned in memory.
263 //
264 // \return \a true in case the operands are aligned, \a false if not.
265 */
isAligned()266 inline bool isAligned() const noexcept {
267 return rhs_.isAligned();
268 }
269 //**********************************************************************************************
270
271 //**********************************************************************************************
272 /*!\brief Returns whether the expression can be used in SMP assignments.
273 //
274 // \return \a true in case the expression can be used in SMP assignments, \a false if not.
275 */
canSMPAssign()276 inline bool canSMPAssign() const noexcept {
277 return false;
278 }
279 //**********************************************************************************************
280
281 private:
282 //**Member variables****************************************************************************
283 LeftOperand lhs_; //!< Left-hand side dense matrix of the Kronecker product expression.
284 RightOperand rhs_; //!< Right-hand side dense matrix of the Kronecker product expression.
285 //**********************************************************************************************
286
287 //**Assignment to row-major dense matrices******************************************************
288 /*! \cond BLAZE_INTERNAL */
289 /*!\brief Assignment of a dense matrix-dense matrix Kronecker product to a row-major dense
290 // matrix.
291 // \ingroup dense_matrix
292 //
293 // \param lhs The target left-hand side dense matrix.
294 // \param rhs The right-hand side Kronecker product expression to be assigned.
295 // \return void
296 //
297 // This function implements the performance optimized assignment of a dense matrix-dense
298 // matrix Kronecker product expression to a row-major dense matrix.
299 */
300 template< typename MT > // Type of the target dense matrix
assign(DenseMatrix<MT,false> & lhs,const DMatDMatKronExpr & rhs)301 friend inline void assign( DenseMatrix<MT,false>& lhs, const DMatDMatKronExpr& rhs )
302 {
303 BLAZE_FUNCTION_TRACE;
304
305 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
306 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
307
308 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
309 return;
310 }
311
312 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
313 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
314
315 const size_t M( B.rows() );
316 const size_t N( B.columns() );
317
318 if( SO )
319 {
320 for( size_t i=0UL; i<A.rows(); ++i ) {
321 for( size_t j=0UL; j<A.columns(); ++j ) {
322 if( !isDefault<strict>( A(i,j) ) ) {
323 for( size_t l=0UL; l<N; ++l )
324 for( size_t k=0UL; k<M; ++k )
325 (*lhs)(i*M+k,j*N+l) = A(i,j) * B(k,l);
326 }
327 else {
328 for( size_t l=0UL; l<N; ++l )
329 for( size_t k=0UL; k<M; ++k )
330 reset( (*lhs)(i*M+k,j*N+l) );
331 }
332 }
333 }
334 }
335 else
336 {
337 for( size_t i=0UL; i<A.rows(); ++i ) {
338 for( size_t k=0UL; k<M; ++k ) {
339 for( size_t j=0UL; j<A.columns(); ++j ) {
340 if( !isDefault<strict>( A(i,j) ) ) {
341 for( size_t l=0UL; l<N; ++l )
342 (*lhs)(i*M+k,j*N+l) = A(i,j) * B(k,l);
343 }
344 else {
345 for( size_t l=0UL; l<N; ++l )
346 reset( (*lhs)(i*M+k,j*N+l) );
347 }
348 }
349 }
350 }
351 }
352 }
353 /*! \endcond */
354 //**********************************************************************************************
355
356 //**Assignment to column-major dense matrices***************************************************
357 /*! \cond BLAZE_INTERNAL */
358 /*!\brief Assignment of a dense matrix-dense matrix Kronecker product to a column-major dense
359 // matrix.
360 // \ingroup dense_matrix
361 //
362 // \param lhs The target left-hand side dense matrix.
363 // \param rhs The right-hand side Kronecker product expression to be assigned.
364 // \return void
365 //
366 // This function implements the performance optimized assignment of a dense matrix-dense
367 // matrix Kronecker product expression to a column-major dense matrix.
368 */
369 template< typename MT > // Type of the target dense matrix
assign(DenseMatrix<MT,true> & lhs,const DMatDMatKronExpr & rhs)370 friend inline void assign( DenseMatrix<MT,true>& lhs, const DMatDMatKronExpr& rhs )
371 {
372 BLAZE_FUNCTION_TRACE;
373
374 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
375 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
376
377 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
378 return;
379 }
380
381 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
382 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
383
384 const size_t M( B.rows() );
385 const size_t N( B.columns() );
386
387 if( SO )
388 {
389 for( size_t j=0UL; j<A.columns(); ++j ) {
390 for( size_t l=0UL; l<N; ++l ) {
391 for( size_t i=0UL; i<A.rows(); ++i ) {
392 if( !isDefault<strict>( A(i,j) ) ) {
393 for( size_t k=0UL; k<M; ++k )
394 (*lhs)(i*M+k,j*N+l) = A(i,j) * B(k,l);
395 }
396 else {
397 for( size_t k=0UL; k<M; ++k )
398 reset( (*lhs)(i*M+k,j*N+l) );
399 }
400 }
401 }
402 }
403 }
404 else
405 {
406 for( size_t j=0UL; j<A.columns(); ++j ) {
407 for( size_t i=0UL; i<A.rows(); ++i ) {
408 if( !isDefault<strict>( A(i,j) ) ) {
409 for( size_t k=0UL; k<M; ++k )
410 for( size_t l=0UL; l<N; ++l )
411 (*lhs)(i*M+k,j*N+l) = A(i,j) * B(k,l);
412 }
413 else {
414 for( size_t k=0UL; k<M; ++k )
415 for( size_t l=0UL; l<N; ++l )
416 reset( (*lhs)(i*M+k,j*N+l) );
417 }
418 }
419 }
420 }
421 }
422 /*! \endcond */
423 //**********************************************************************************************
424
425 //**Assignment to sparse matrices***************************************************************
426 /*! \cond BLAZE_INTERNAL */
427 /*!\brief Assignment of a dense matrix-dense matrix Kronecker product to a sparse matrix.
428 // \ingroup dense_matrix
429 //
430 // \param lhs The target left-hand side sparse matrix.
431 // \param rhs The right-hand side Kronecker product expression to be assigned.
432 // \return void
433 //
434 // This function implements the performance optimized assignment of a dense matrix-dense
435 // matrix Kronecker product expression to a sparse matrix.
436 */
437 template< typename MT // Type of the target sparse matrix
438 , bool SO2 > // Storage order of the target sparse matrix
assign(SparseMatrix<MT,SO2> & lhs,const DMatDMatKronExpr & rhs)439 friend inline void assign( SparseMatrix<MT,SO2>& lhs, const DMatDMatKronExpr& rhs )
440 {
441 BLAZE_FUNCTION_TRACE;
442
443 using TmpType = If_t< SO == SO2, ResultType, OppositeType >;
444
445 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( ResultType );
446 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( OppositeType );
447 BLAZE_CONSTRAINT_MUST_BE_MATRIX_WITH_STORAGE_ORDER( ResultType, SO );
448 BLAZE_CONSTRAINT_MUST_BE_MATRIX_WITH_STORAGE_ORDER( OppositeType, !SO );
449 BLAZE_CONSTRAINT_MATRICES_MUST_HAVE_SAME_STORAGE_ORDER( MT, TmpType );
450 BLAZE_CONSTRAINT_MUST_NOT_REQUIRE_EVALUATION( TmpType );
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 const TmpType tmp( serial( rhs ) );
456 assign( *lhs, tmp );
457 }
458 /*! \endcond */
459 //**********************************************************************************************
460
461 //**Addition assignment to row-major dense matrices*********************************************
462 /*! \cond BLAZE_INTERNAL */
463 /*!\brief Addition assignment of a dense matrix-dense matrix Kronecker product to a row-major
464 // dense matrix.
465 // \ingroup dense_matrix
466 //
467 // \param lhs The target left-hand side dense matrix.
468 // \param rhs The right-hand side Kronecker product expression to be added.
469 // \return void
470 //
471 // This function implements the performance optimized addition assignment of a dense matrix-
472 // dense matrix Kronecker product expression to a row-major dense matrix.
473 */
474 template< typename MT > // Type of the target dense matrix
addAssign(DenseMatrix<MT,false> & lhs,const DMatDMatKronExpr & rhs)475 friend inline void addAssign( DenseMatrix<MT,false>& lhs, const DMatDMatKronExpr& rhs )
476 {
477 BLAZE_FUNCTION_TRACE;
478
479 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
480 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
481
482 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
483 return;
484 }
485
486 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
487 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
488
489 const size_t M( B.rows() );
490 const size_t N( B.columns() );
491
492 if( SO )
493 {
494 for( size_t i=0UL; i<A.rows(); ++i )
495 for( size_t j=0UL; j<A.columns(); ++j )
496 if( !isDefault<strict>( A(i,j) ) )
497 for( size_t l=0UL; l<N; ++l )
498 {
499 const size_t kbegin( ( IsLower_v<MT2> )
500 ?( ( IsStrictlyLower_v<MT2> ? l+1UL : l ) )
501 :( 0UL ) );
502 const size_t kend( ( IsUpper_v<MT2> )
503 ?( IsStrictlyUpper_v<MT2> ? l : l+1UL )
504 :( M ) );
505 BLAZE_INTERNAL_ASSERT( kbegin <= kend, "Invalid loop indices detected" );
506
507 for( size_t k=kbegin; k<kend; ++k )
508 (*lhs)(i*M+k,j*N+l) += A(i,j) * B(k,l);
509 }
510 }
511 else
512 {
513 for( size_t i=0UL; i<A.rows(); ++i )
514 for( size_t k=0UL; k<M; ++k )
515 {
516 const size_t lbegin( ( IsUpper_v<MT2> )
517 ?( ( IsStrictlyUpper_v<MT2> ? k+1UL : k ) )
518 :( 0UL ) );
519 const size_t lend( ( IsLower_v<MT2> )
520 ?( IsStrictlyLower_v<MT2> ? k : k+1UL )
521 :( N ) );
522 BLAZE_INTERNAL_ASSERT( lbegin <= lend, "Invalid loop indices detected" );
523
524 for( size_t j=0UL; j<A.columns(); ++j )
525 if( !isDefault<strict>( A(i,j) ) )
526 for( size_t l=lbegin; l<lend; ++l )
527 (*lhs)(i*M+k,j*N+l) += A(i,j) * B(k,l);
528 }
529 }
530 }
531 /*! \endcond */
532 //**********************************************************************************************
533
534 //**Addition assignment to column-major dense matrices******************************************
535 /*! \cond BLAZE_INTERNAL */
536 /*!\brief Addition assignment of a dense matrix-dense matrix Kronecker product to a column-major
537 // dense matrix.
538 // \ingroup dense_matrix
539 //
540 // \param lhs The target left-hand side dense matrix.
541 // \param rhs The right-hand side Kronecker product expression to be added.
542 // \return void
543 //
544 // This function implements the performance optimized addition assignment of a dense matrix-
545 // dense matrix Kronecker product expression to a column-major dense matrix.
546 */
547 template< typename MT > // Type of the target dense matrix
addAssign(DenseMatrix<MT,true> & lhs,const DMatDMatKronExpr & rhs)548 friend inline void addAssign( DenseMatrix<MT,true>& lhs, const DMatDMatKronExpr& rhs )
549 {
550 BLAZE_FUNCTION_TRACE;
551
552 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
553 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
554
555 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
556 return;
557 }
558
559 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
560 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
561
562 const size_t M( B.rows() );
563 const size_t N( B.columns() );
564
565 if( SO )
566 {
567 for( size_t j=0UL; j<A.columns(); ++j )
568 for( size_t l=0UL; l<N; ++l )
569 {
570 const size_t kbegin( ( IsLower_v<MT2> )
571 ?( ( IsStrictlyLower_v<MT2> ? l+1UL : l ) )
572 :( 0UL ) );
573 const size_t kend( ( IsUpper_v<MT2> )
574 ?( IsStrictlyUpper_v<MT2> ? l : l+1UL )
575 :( M ) );
576 BLAZE_INTERNAL_ASSERT( kbegin <= kend, "Invalid loop indices detected" );
577
578 for( size_t i=0UL; i<A.rows(); ++i )
579 if( !isDefault<strict>( A(i,j) ) )
580 for( size_t k=kbegin; k<kend; ++k )
581 (*lhs)(i*M+k,j*N+l) += A(i,j) * B(k,l);
582 }
583 }
584 else
585 {
586 for( size_t j=0UL; j<A.columns(); ++j )
587 for( size_t i=0UL; i<A.rows(); ++i )
588 if( !isDefault<strict>( A(i,j) ) )
589 for( size_t k=0UL; k<M; ++k )
590 {
591 const size_t lbegin( ( IsUpper_v<MT2> )
592 ?( ( IsStrictlyUpper_v<MT2> ? k+1UL : k ) )
593 :( 0UL ) );
594 const size_t lend( ( IsLower_v<MT2> )
595 ?( IsStrictlyLower_v<MT2> ? k : k+1UL )
596 :( N ) );
597 BLAZE_INTERNAL_ASSERT( lbegin <= lend, "Invalid loop indices detected" );
598
599 for( size_t l=lbegin; l<lend; ++l )
600 (*lhs)(i*M+k,j*N+l) += A(i,j) * B(k,l);
601 }
602 }
603 }
604 /*! \endcond */
605 //**********************************************************************************************
606
607 //**Addition assignment to sparse matrices******************************************************
608 // No special implementation for the addition assignment to sparse matrices.
609 //**********************************************************************************************
610
611 //**Subtraction assignment to row-major dense matrices******************************************
612 /*! \cond BLAZE_INTERNAL */
613 /*!\brief Subtraction assignment of a dense matrix-dense matrix Kronecker product to a row-major
614 // dense matrix.
615 // \ingroup dense_matrix
616 //
617 // \param lhs The target left-hand side dense matrix.
618 // \param rhs The right-hand side Kronecker product expression to be subtracted.
619 // \return void
620 //
621 // This function implements the performance optimized subtracted assignment of a row-major
622 // dense matrix-dense matrix Kronecker product expression to a dense matrix.
623 */
624 template< typename MT > // Type of the target dense matrix
subAssign(DenseMatrix<MT,false> & lhs,const DMatDMatKronExpr & rhs)625 friend inline void subAssign( DenseMatrix<MT,false>& lhs, const DMatDMatKronExpr& rhs )
626 {
627 BLAZE_FUNCTION_TRACE;
628
629 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
630 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
631
632 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
633 return;
634 }
635
636 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
637 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
638
639 const size_t M( B.rows() );
640 const size_t N( B.columns() );
641
642 if( SO )
643 {
644 for( size_t i=0UL; i<A.rows(); ++i )
645 for( size_t j=0UL; j<A.columns(); ++j )
646 if( !isDefault<strict>( A(i,j) ) )
647 for( size_t l=0UL; l<N; ++l )
648 {
649 const size_t kbegin( ( IsLower_v<MT2> )
650 ?( ( IsStrictlyLower_v<MT2> ? l+1UL : l ) )
651 :( 0UL ) );
652 const size_t kend( ( IsUpper_v<MT2> )
653 ?( IsStrictlyUpper_v<MT2> ? l : l+1UL )
654 :( M ) );
655 BLAZE_INTERNAL_ASSERT( kbegin <= kend, "Invalid loop indices detected" );
656
657 for( size_t k=kbegin; k<kend; ++k )
658 (*lhs)(i*M+k,j*N+l) -= A(i,j) * B(k,l);
659 }
660 }
661 else
662 {
663 for( size_t i=0UL; i<A.rows(); ++i )
664 for( size_t k=0UL; k<M; ++k )
665 {
666 const size_t lbegin( ( IsUpper_v<MT2> )
667 ?( ( IsStrictlyUpper_v<MT2> ? k+1UL : k ) )
668 :( 0UL ) );
669 const size_t lend( ( IsLower_v<MT2> )
670 ?( IsStrictlyLower_v<MT2> ? k : k+1UL )
671 :( N ) );
672 BLAZE_INTERNAL_ASSERT( lbegin <= lend, "Invalid loop indices detected" );
673
674 for( size_t j=0UL; j<A.columns(); ++j )
675 if( !isDefault<strict>( A(i,j) ) )
676 for( size_t l=lbegin; l<lend; ++l )
677 (*lhs)(i*M+k,j*N+l) -= A(i,j) * B(k,l);
678 }
679 }
680 }
681 /*! \endcond */
682 //**********************************************************************************************
683
684 //**Subtraction assignment to column-major dense matrices***************************************
685 /*! \cond BLAZE_INTERNAL */
686 /*!\brief Subtraction assignment of a dense matrix-dense matrix Kronecker product to a
687 // column-major dense matrix.
688 // \ingroup dense_matrix
689 //
690 // \param lhs The target left-hand side dense matrix.
691 // \param rhs The right-hand side Kronecker product expression to be subtracted.
692 // \return void
693 //
694 // This function implements the performance optimized subtracted assignment of a column-major
695 // dense matrix-dense matrix Kronecker product expression to a dense matrix.
696 */
697 template< typename MT > // Type of the target dense matrix
subAssign(DenseMatrix<MT,true> & lhs,const DMatDMatKronExpr & rhs)698 friend inline void subAssign( DenseMatrix<MT,true>& lhs, const DMatDMatKronExpr& rhs )
699 {
700 BLAZE_FUNCTION_TRACE;
701
702 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
703 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
704
705 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
706 return;
707 }
708
709 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
710 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
711
712 const size_t M( B.rows() );
713 const size_t N( B.columns() );
714
715 if( SO )
716 {
717 for( size_t j=0UL; j<A.columns(); ++j )
718 for( size_t l=0UL; l<N; ++l )
719 {
720 const size_t kbegin( ( IsLower_v<MT2> )
721 ?( ( IsStrictlyLower_v<MT2> ? l+1UL : l ) )
722 :( 0UL ) );
723 const size_t kend( ( IsUpper_v<MT2> )
724 ?( IsStrictlyUpper_v<MT2> ? l : l+1UL )
725 :( M ) );
726 BLAZE_INTERNAL_ASSERT( kbegin <= kend, "Invalid loop indices detected" );
727
728 for( size_t i=0UL; i<A.rows(); ++i )
729 if( !isDefault<strict>( A(i,j) ) )
730 for( size_t k=kbegin; k<kend; ++k )
731 (*lhs)(i*M+k,j*N+l) -= A(i,j) * B(k,l);
732 }
733 }
734 else
735 {
736 for( size_t j=0UL; j<A.columns(); ++j )
737 for( size_t i=0UL; i<A.rows(); ++i )
738 if( !isDefault<strict>( A(i,j) ) )
739 for( size_t k=0UL; k<M; ++k )
740 {
741 const size_t lbegin( ( IsUpper_v<MT2> )
742 ?( ( IsStrictlyUpper_v<MT2> ? k+1UL : k ) )
743 :( 0UL ) );
744 const size_t lend( ( IsLower_v<MT2> )
745 ?( IsStrictlyLower_v<MT2> ? k : k+1UL )
746 :( N ) );
747 BLAZE_INTERNAL_ASSERT( lbegin <= lend, "Invalid loop indices detected" );
748
749 for( size_t l=lbegin; l<lend; ++l )
750 (*lhs)(i*M+k,j*N+l) -= A(i,j) * B(k,l);
751 }
752 }
753 }
754 /*! \endcond */
755 //**********************************************************************************************
756
757 //**Subtraction assignment to sparse matrices***************************************************
758 // No special implementation for the subtraction assignment to sparse matrices.
759 //**********************************************************************************************
760
761 //**Schur product assignment to row-major dense matrices****************************************
762 /*! \cond BLAZE_INTERNAL */
763 /*!\brief Schur product assignment of a dense matrix-dense matrix Kronecker product to a
764 // row-major dense matrix.
765 // \ingroup dense_matrix
766 //
767 // \param lhs The target left-hand side dense matrix.
768 // \param rhs The right-hand side Kronecker product expression for the Schur product.
769 // \return void
770 //
771 // This function implements the performance optimized Schur product assignment of a dense
772 // matrix-dense matrix Kronecker product expression to a row-major dense matrix.
773 */
774 template< typename MT > // Type of the target dense matrix
schurAssign(DenseMatrix<MT,false> & lhs,const DMatDMatKronExpr & rhs)775 friend inline void schurAssign( DenseMatrix<MT,false>& lhs, const DMatDMatKronExpr& rhs )
776 {
777 BLAZE_FUNCTION_TRACE;
778
779 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
780 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
781
782 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
783 return;
784 }
785
786 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
787 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
788
789 const size_t M( B.rows() );
790 const size_t N( B.columns() );
791
792 if( SO )
793 {
794 for( size_t i=0UL; i<A.rows(); ++i ) {
795 for( size_t j=0UL; j<A.columns(); ++j ) {
796 if( !isDefault<strict>( A(i,j) ) ) {
797 for( size_t l=0UL; l<N; ++l )
798 for( size_t k=0UL; k<M; ++k )
799 (*lhs)(i*M+k,j*N+l) *= A(i,j) * B(k,l);
800 }
801 else {
802 for( size_t l=0UL; l<N; ++l )
803 for( size_t k=0UL; k<M; ++k )
804 reset( (*lhs)(i*M+k,j*N+l) );
805 }
806 }
807 }
808 }
809 else
810 {
811 for( size_t i=0UL; i<A.rows(); ++i ) {
812 for( size_t k=0UL; k<M; ++k ) {
813 for( size_t j=0UL; j<A.columns(); ++j ) {
814 if( !isDefault<strict>( A(i,j) ) ) {
815 for( size_t l=0UL; l<N; ++l )
816 (*lhs)(i*M+k,j*N+l) *= A(i,j) * B(k,l);
817 }
818 else {
819 for( size_t l=0UL; l<N; ++l )
820 reset( (*lhs)(i*M+k,j*N+l) );
821 }
822 }
823 }
824 }
825 }
826 }
827 /*! \endcond */
828 //**********************************************************************************************
829
830 //**Schur product assignment to column-major dense matrices*************************************
831 /*! \cond BLAZE_INTERNAL */
832 /*!\brief Schur product assignment of a dense matrix-dense matrix Kronecker product to a
833 // column-major dense matrix.
834 // \ingroup dense_matrix
835 //
836 // \param lhs The target left-hand side dense matrix.
837 // \param rhs The right-hand side Kronecker product expression for the Schur product.
838 // \return void
839 //
840 // This function implements the performance optimized Schur product assignment of a dense
841 // matrix-dense matrix Kronecker product expression to a column-major dense matrix.
842 */
843 template< typename MT > // Type of the target dense matrix
schurAssign(DenseMatrix<MT,true> & lhs,const DMatDMatKronExpr & rhs)844 friend inline void schurAssign( DenseMatrix<MT,true>& lhs, const DMatDMatKronExpr& rhs )
845 {
846 BLAZE_FUNCTION_TRACE;
847
848 BLAZE_INTERNAL_ASSERT( (*lhs).rows() == rhs.rows() , "Invalid number of rows" );
849 BLAZE_INTERNAL_ASSERT( (*lhs).columns() == rhs.columns(), "Invalid number of columns" );
850
851 if( rhs.rows() == 0UL || rhs.columns() == 0UL ) {
852 return;
853 }
854
855 CT1 A( serial( rhs.lhs_ ) ); // Evaluation of the left-hand side dense matrix operand
856 CT2 B( serial( rhs.rhs_ ) ); // Evaluation of the right-hand side dense matrix operand
857
858 const size_t M( B.rows() );
859 const size_t N( B.columns() );
860
861 if( SO )
862 {
863 for( size_t j=0UL; j<A.columns(); ++j ) {
864 for( size_t l=0UL; l<N; ++l ) {
865 for( size_t i=0UL; i<A.rows(); ++i ) {
866 if( !isDefault<strict>( A(i,j) ) ) {
867 for( size_t k=0UL; k<M; ++k )
868 (*lhs)(i*M+k,j*N+l) *= A(i,j) * B(k,l);
869 }
870 else {
871 for( size_t k=0UL; k<M; ++k )
872 reset( (*lhs)(i*M+k,j*N+l) );
873 }
874 }
875 }
876 }
877 }
878 else
879 {
880 for( size_t j=0UL; j<A.columns(); ++j ) {
881 for( size_t i=0UL; i<A.rows(); ++i ) {
882 if( !isDefault<strict>( A(i,j) ) ) {
883 for( size_t k=0UL; k<M; ++k )
884 for( size_t l=0UL; l<N; ++l )
885 (*lhs)(i*M+k,j*N+l) *= A(i,j) * B(k,l);
886 }
887 else {
888 for( size_t k=0UL; k<M; ++k )
889 for( size_t l=0UL; l<N; ++l )
890 reset( (*lhs)(i*M+k,j*N+l) );
891 }
892 }
893 }
894 }
895 }
896 /*! \endcond */
897 //**********************************************************************************************
898
899 //**Schur product assignment to sparse matrices*************************************************
900 // No special implementation for the Schur product assignment to sparse matrices.
901 //**********************************************************************************************
902
903 //**Multiplication assignment to dense matrices*************************************************
904 // No special implementation for the multiplication assignment to dense matrices.
905 //**********************************************************************************************
906
907 //**Multiplication assignment to sparse matrices************************************************
908 // No special implementation for the multiplication assignment to sparse matrices.
909 //**********************************************************************************************
910
911 //**Compile time checks*************************************************************************
912 /*! \cond BLAZE_INTERNAL */
913 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( MT1 );
914 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( MT2 );
915 BLAZE_CONSTRAINT_MUST_BE_MATRIX_WITH_STORAGE_ORDER( MT2, SO );
916 BLAZE_CONSTRAINT_MUST_FORM_VALID_MATMATKRONEXPR( MT1, MT2 );
917 /*! \endcond */
918 //**********************************************************************************************
919 };
920 //*************************************************************************************************
921
922
923
924
925 //=================================================================================================
926 //
927 // GLOBAL FUNCTIONS
928 //
929 //=================================================================================================
930
931 //*************************************************************************************************
932 /*!\brief Computes the Kronecker product of two dense matrices (\f$ A=B \otimes C \f$).
933 // \ingroup dense_matrix
934 //
935 // \param lhs The left-hand side dense matrix for the Kronecker product.
936 // \param rhs The right-hand side dense matrix for the Kronecker product.
937 // \return The Kronecker product of the two matrices.
938 //
939 // The kron() function computes the Kronecker product of the two given dense matrices:
940
941 \code
942 blaze::DynamicMatrix<double> A, B, C;
943 // ... Resizing and initialization
944 C = kron( A, B );
945 \endcode
946
947 // The function returns an expression representing a dense matrix of the higher-order element
948 // type of the two involved matrix element types \a MT1::ElementType and \a MT2::ElementType.
949 // Both matrix types \a MT1 and \a MT2 as well as the two element types \a MT1::ElementType
950 // and \a MT2::ElementType have to be supported by the MultTrait class template.
951 */
952 template< typename MT1 // Type of the left-hand side dense matrix
953 , bool SO1 // Storage order of the left-hand side dense matrix
954 , typename MT2 // Type of the right-hand side dense matrix
955 , bool SO2 > // Storage order of the right-hand side dense matrix
decltype(auto)956 inline decltype(auto)
957 kron( const DenseMatrix<MT1,SO1>& lhs, const DenseMatrix<MT2,SO2>& rhs )
958 {
959 BLAZE_FUNCTION_TRACE;
960
961 using ReturnType = const DMatDMatKronExpr<MT1,MT2,SO2>;
962 return ReturnType( *lhs, *rhs );
963 }
964 //*************************************************************************************************
965
966 } // namespace blaze
967
968 #endif
969