Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DgemmOverride.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // For more information, please see: http://www.nektar.info
4 //
5 // The MIT License
6 //
7 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
8 // Department of Aeronautics, Imperial College London (UK), and Scientific
9 // Computing and Imaging Institute, University of Utah (USA).
10 //
11 // License for the specific language governing rights and limitations under
12 // Permission is hereby granted, free of charge, to any person obtaining a
13 // copy of this software and associated documentation files (the "Software"),
14 // to deal in the Software without restriction, including without limitation
15 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 // and/or sell copies of the Software, and to permit persons to whom the
17 // Software is furnished to do so, subject to the following conditions:
18 //
19 // The above copyright notice and this permission notice shall be included
20 // in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 // DEALINGS IN THE SOFTWARE.
29 //
30 // Description: Defines the global functions needed for matrix operations.
31 //
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_DGEMM_OVERRIDE_HPP
35 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_DGEMM_OVERRIDE_HPP
36 
37 #ifdef NEKTAR_USE_EXPRESSION_TEMPLATES
38 
39 #include <ExpressionTemplates/ExpressionTemplates.hpp>
40 #include <boost/utility/enable_if.hpp>
41 #include <boost/type_traits.hpp>
42 #include <LibUtilities/LinearAlgebra/CanGetRawPtr.hpp>
43 #include <LibUtilities/LinearAlgebra/Blas.hpp>
44 #include <LibUtilities/LinearAlgebra/BlasOverrideUtil.hpp>
45 
46 namespace Nektar
47 {
48  template<typename ADataType, typename BDataType,
49  typename AMatrixType, typename BMatrixType>
50  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
51  double alpha, const NekMatrix<ADataType, AMatrixType>& A, const NekMatrix<BDataType, BMatrixType>& B)
52  {
53  if( A.GetType() != eFULL || B.GetType() != eFULL )
54  {
55  Multiply(result, A, B);
56  MultiplyEqual(result, alpha);
57  return;
58  }
59 
60  unsigned int M = A.GetRows();
61  unsigned int N = B.GetColumns();
62  unsigned int K = A.GetColumns();
63 
64  unsigned int LDA = M;
65  if( A.GetTransposeFlag() == 'T' )
66  {
67  LDA = K;
68  }
69 
70  unsigned int LDB = K;
71  if( B.GetTransposeFlag() == 'T' )
72  {
73  LDB = N;
74  }
75 
76  Blas::Dgemm(A.GetTransposeFlag(), B.GetTransposeFlag(), M, N, K,
77  A.Scale()*B.Scale()*alpha, A.GetRawPtr(), LDA,
78  B.GetRawPtr(), LDB, 0.0,
79  result.GetRawPtr(), result.GetRows());
80  }
81 
82  template<typename ADataType, typename BDataType,
83  typename AMatrixType, typename BMatrixType>
84  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
85  const NekMatrix<ADataType, AMatrixType>& A, double alpha, const NekMatrix<BDataType, BMatrixType>& B)
86  {
87  Dgemm(result, alpha, A, B);
88  }
89 
90  template<typename ADataType, typename BDataType,
91  typename AMatrixType, typename BMatrixType>
92  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
93  const NekMatrix<ADataType, AMatrixType>& A, const NekMatrix<BDataType, BMatrixType>& B, double alpha)
94  {
95  Dgemm(result, alpha, A, B);
96  }
97 
98  template<typename ADataType, typename BDataType, typename CDataType,
99  typename AMatrixType, typename BMatrixType, typename CMatrixType>
100  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
101  double alpha, const NekMatrix<ADataType, AMatrixType>& A, const NekMatrix<BDataType, BMatrixType>& B,
102  double beta, const NekMatrix<CDataType, CMatrixType>& C)
103  {
104  if( A.GetType() != eFULL || B.GetType() != eFULL || C.GetType() != eFULL )
105  {
106  Multiply(result, A, B);
107  MultiplyEqual(result, alpha);
108  NekMatrix<double> temp = beta*C;
109  AddEqual(result, temp);
110  return;
111  }
112 
113  if( C.GetTransposeFlag() == 'T' )
114  {
115  // In this case, the data array in C is not ordered correctly,
116  // so we need to copy it into result to get the right ordering.
117  result.SetSize(C.GetRows(), C.GetColumns());
118 
119  for(unsigned int i = 0; i < C.GetRows(); ++i)
120  {
121  for(unsigned int j = 0; j < C.GetColumns(); ++j)
122  {
123  result(i,j) = C(i,j);
124  }
125  }
126  }
127  else
128  {
129  result = C;
130  }
131 
132  unsigned int M = A.GetRows();
133  unsigned int N = B.GetColumns();
134  unsigned int K = A.GetColumns();
135 
136  unsigned int LDA = M;
137  if( A.GetTransposeFlag() == 'T' )
138  {
139  LDA = K;
140  }
141 
142  unsigned int LDB = K;
143  if( B.GetTransposeFlag() == 'T' )
144  {
145  LDB = N;
146  }
147 
148  Blas::Dgemm(A.GetTransposeFlag(), B.GetTransposeFlag(), M, N, K,
149  A.Scale()*B.Scale()*alpha, A.GetRawPtr(), LDA,
150  B.GetRawPtr(), LDB, beta*result.Scale(),
151  result.GetRawPtr(), result.GetRows());
152  }
153 
154  template<typename ADataType, typename BDataType, typename CDataType,
155  typename AMatrixType, typename BMatrixType, typename CMatrixType>
156  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
157  const NekMatrix<ADataType, AMatrixType>& A, double alpha, const NekMatrix<BDataType, BMatrixType>& B,
158  double beta, const NekMatrix<CDataType, CMatrixType>& C)
159  {
160  Dgemm(result, alpha, A, B, beta, C);
161  }
162 
163  template<typename ADataType, typename BDataType, typename CDataType,
164  typename AMatrixType, typename BMatrixType, typename CMatrixType>
165  void Dgemm(NekMatrix<double, StandardMatrixTag>& result,
166  const NekMatrix<ADataType, AMatrixType>& A, const NekMatrix<BDataType, BMatrixType>& B, double alpha,
167  double beta, const NekMatrix<CDataType, CMatrixType>& C)
168  {
169  Dgemm(result, alpha, A, B, beta, C);
170  }
171 }
172 
173 namespace expt
174 {
175  namespace impl
176  {
177  // The DgemmNodeEvaluator creates an evaluator for the requested node, or the
178  // child of a negate node. The negation operation will occur as part of the dgemm
179  // call, so we do not need to do it up front.
180  template<typename NodeType, typename IndicesType, unsigned int index>
181  struct DgemmNodeEvaluator
182  {
183  typedef EvaluateNodeWithTemporaryIfNeeded<NodeType, IndicesType, index> Evaluator;
184  typedef typename Evaluator::ResultType Type;
185  static double GetScale() { return 1.0; }
186  };
187 
188  template<typename LhsType, typename IndicesType, unsigned int index>
189  struct DgemmNodeEvaluator<Node<LhsType, NegateOp>, IndicesType, index>
190  {
191  typedef EvaluateNodeWithTemporaryIfNeeded<LhsType, IndicesType, index> Evaluator;
192  typedef typename Evaluator::ResultType Type;
193  static double GetScale() { return -1.0; }
194  };
195  }
196 
197  namespace impl
198  {
199  // To handle alpha*A*B - beta*C, this class provides the scale depending on the operator.
200  template<typename OpType>
201  struct BetaScale
202  {
203  static double GetScale() { return 1.0; }
204  };
205 
206  template<>
207  struct BetaScale<SubtractOp>
208  {
209  static double GetScale() { return -1.0; }
210  };
211 
212  template<typename T>
213  struct IsDouble : public boost::false_type{};
214 
215  template<>
216  struct IsDouble<double> : public boost::true_type {};
217 
218  }
219 
220  namespace impl
221  {
222  template<typename NodeType, typename IndicesType, unsigned int index, typename enabled=void>
223  struct AlphaABParameterAccessImpl : public boost::false_type {};
224 
225  template<typename A1, typename A2, typename A3, typename IndicesType, unsigned int index>
226  struct AlphaABParameterAccessImpl< Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, IndicesType, index,
227  typename boost::enable_if
228  <
229  Test3ArgumentAssociativeNode<Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, IsDouble, MultiplyOp,
230  Nektar::CanGetRawPtr, MultiplyOp, Nektar::CanGetRawPtr>
231  //boost::mpl::and_
232  //<
233  // boost::is_same<typename A1::ResultType, double>,
234  // CanGetRawPtr<typename A2::ResultType>,
235  // CanGetRawPtr<typename A3::ResultType>
236  //>
237  >::type> : public boost::true_type
238  {
239  typedef DgemmNodeEvaluator<A1, IndicesType, index> AlphaWrappedEvaluator;
240  typedef typename AlphaWrappedEvaluator::Evaluator AlphaEvaluator;
241 
242  static const unsigned int A2Index = index + A1::TotalCount;
243  typedef DgemmNodeEvaluator<A2, IndicesType, A2Index> AWrappedEvaluator;
244  typedef typename AWrappedEvaluator::Evaluator AEvaluator;
245 
246  static const unsigned int A3Index = A2Index + A2::TotalCount;
247  typedef DgemmNodeEvaluator<A3, IndicesType, A3Index> BWrappedEvaluator;
248  typedef typename BWrappedEvaluator::Evaluator BEvaluator;
249 
250  template<typename ArgumentVectorType>
251  static double GetAlpha(const ArgumentVectorType& args)
252  {
253  return AlphaEvaluator::Evaluate(args) * AWrappedEvaluator::GetScale() * BWrappedEvaluator::GetScale();
254  }
255  };
256 
257  template<typename A1, typename A2, typename A3, typename IndicesType, unsigned int index>
258  struct AlphaABParameterAccessImpl< Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, IndicesType, index,
259  typename boost::enable_if
260  <
261  Test3ArgumentAssociativeNode<Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, Nektar::CanGetRawPtr, MultiplyOp,
262  IsDouble, MultiplyOp, Nektar::CanGetRawPtr>
263  //boost::mpl::and_
264  //<
265  // CanGetRawPtr<typename A1::ResultType>,
266  // boost::is_same<typename A2::ResultType, double>,
267  // CanGetRawPtr<typename A3::ResultType>
268  //>
269  >::type> : public boost::true_type
270  {
271  typedef DgemmNodeEvaluator<A1, IndicesType, index> AWrappedEvaluator;
272  typedef typename AWrappedEvaluator::Evaluator AEvaluator;
273 
274  static const unsigned int A2Index = index + A1::TotalCount;
275  typedef DgemmNodeEvaluator<A2, IndicesType, A2Index> AlphaWrappedEvaluator;
276  typedef typename AlphaWrappedEvaluator::Evaluator AlphaEvaluator;
277 
278  static const unsigned int A3Index = A2Index + A2::TotalCount;
279  typedef DgemmNodeEvaluator<A3, IndicesType, A3Index> BWrappedEvaluator;
280  typedef typename BWrappedEvaluator::Evaluator BEvaluator;
281 
282  template<typename ArgumentVectorType>
283  static double GetAlpha(const ArgumentVectorType& args)
284  {
285  return AlphaEvaluator::Evaluate(args) * AWrappedEvaluator::GetScale() * BWrappedEvaluator::GetScale();
286  }
287  };
288 
289  template<typename A1, typename A2, typename A3, typename IndicesType, unsigned int index>
290  struct AlphaABParameterAccessImpl< Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, IndicesType, index,
291  typename boost::enable_if
292  <
293  Test3ArgumentAssociativeNode<Node<Node<A1, MultiplyOp, A2>, MultiplyOp, A3>, Nektar::CanGetRawPtr, MultiplyOp,
294  Nektar::CanGetRawPtr, MultiplyOp, IsDouble>
295  //boost::mpl::and_
296  //<
297  // CanGetRawPtr<typename A1::ResultType>,
298  // CanGetRawPtr<typename A2::ResultType>,
299  // boost::is_same<typename A3::ResultType, double>
300  //>
301  >::type> : public boost::true_type
302  {
303  typedef DgemmNodeEvaluator<A1, IndicesType, index> AWrappedEvaluator;
304  typedef typename AWrappedEvaluator::Evaluator AEvaluator;
305 
306  static const unsigned int A2Index = index + A1::TotalCount;
307  typedef DgemmNodeEvaluator<A2, IndicesType, A2Index> BWrappedEvaluator;
308  typedef typename BWrappedEvaluator::Evaluator BEvaluator;
309 
310  static const unsigned int A3Index = A2Index + A2::TotalCount;
311  typedef DgemmNodeEvaluator<A3, IndicesType, A3Index> AlphaWrappedEvaluator;
312  typedef typename AlphaWrappedEvaluator::Evaluator AlphaEvaluator;
313 
314  template<typename ArgumentVectorType>
315  static double GetAlpha(const ArgumentVectorType& args)
316  {
317  return AlphaEvaluator::Evaluate(args) * AWrappedEvaluator::GetScale() * BWrappedEvaluator::GetScale();
318  }
319  };
320 
321  // Plain M*M
322  template<typename A1, typename A2, typename IndicesType, unsigned int index>
323  struct AlphaABParameterAccessImpl< Node<A1, MultiplyOp, A2>, IndicesType, index,
324  typename boost::enable_if
325  <
326  boost::mpl::and_
327  <
328  TestBinaryNode<Node<A1, MultiplyOp, A2>, Nektar::CanGetRawPtr, MultiplyOp, Nektar::CanGetRawPtr>,
329  boost::mpl::not_<Test3ArgumentAssociativeNode<Node<A1, MultiplyOp, A2>, Nektar::CanGetRawPtr, MultiplyOp, Nektar::CanGetRawPtr, MultiplyOp, IsDouble> >,
330  boost::mpl::not_<Test3ArgumentAssociativeNode<Node<A1, MultiplyOp, A2>, Nektar::CanGetRawPtr, MultiplyOp, IsDouble, MultiplyOp, Nektar::CanGetRawPtr> >,
331  boost::mpl::not_<Test3ArgumentAssociativeNode<Node<A1, MultiplyOp, A2>, IsDouble, MultiplyOp, Nektar::CanGetRawPtr, MultiplyOp, Nektar::CanGetRawPtr> >
332  >
333  //boost::mpl::and_
334  //<
335  // CanGetRawPtr<typename A1::ResultType>,
336  // CanGetRawPtr<typename A2::ResultType>,
337  // boost::is_same<typename A3::ResultType, double>
338  //>
339  >::type> : public boost::true_type
340  {
341  typedef DgemmNodeEvaluator<A1, IndicesType, index> AWrappedEvaluator;
342  typedef typename AWrappedEvaluator::Evaluator AEvaluator;
343 
344  static const unsigned int A2Index = index + A1::TotalCount;
345  typedef DgemmNodeEvaluator<A2, IndicesType, A2Index> BWrappedEvaluator;
346  typedef typename BWrappedEvaluator::Evaluator BEvaluator;
347 
348  template<typename ArgumentVectorType>
349  static double GetAlpha(const ArgumentVectorType& args)
350  {
351  return 1.0 * AWrappedEvaluator::GetScale() * BWrappedEvaluator::GetScale();
352  }
353  };
354 
355  // The separation here allows us to separate the negatin from the AlphaABNode.
356  template<typename NodeType, typename IndicesType, unsigned int index>
357  struct AlphaABParameterAccess : public AlphaABParameterAccessImpl<NodeType, IndicesType, index> {};
358 
359  template<typename T, typename IndicesType, unsigned int index>
360  struct AlphaABParameterAccess<Node<T, NegateOp, void>, IndicesType, index> : public AlphaABParameterAccessImpl<T, IndicesType, index>
361  {
362  typedef Node<T, NegateOp, void> NodeType;
363 
364  template<typename ArgumentVectorType>
365  static double GetAlpha(const ArgumentVectorType& args)
366  {
367  return -AlphaABParameterAccessImpl<T, IndicesType, index>::GetAlpha(args);
368  }
369  };
370 
371 
372  template<typename NodeType, typename IndicesType, unsigned int index, typename enabled=void>
373  struct BetaCParameterAccessImpl : public boost::false_type {};
374 
375  template<typename T, typename IndicesType, unsigned int index>
376  struct BetaCParameterAccessImpl< Node<T, void, void>, IndicesType, index> : public boost::true_type
377  {
378  typedef Node<T, void, void> NodeType;
379  typedef DgemmNodeEvaluator<NodeType, IndicesType, index> CWrappedEvaluator;
380  typedef typename CWrappedEvaluator::Evaluator CEvaluator;
381 
382  template<typename ArgumentVectorType>
383  static double GetBeta(const ArgumentVectorType& args)
384  {
385  return 1.0;
386  }
387  };
388 
389  template<typename L, typename R, typename IndicesType, unsigned int index>
390  struct BetaCParameterAccessImpl< Node<L, expt::MultiplyOp, R>, IndicesType, index,
391  typename boost::enable_if
392  <
393  boost::mpl::and_
394  <
395  Nektar::CanGetRawPtr<typename L::ResultType>,
396  boost::is_same<typename R::ResultType, double>
397  >
398  >::type> : public boost::true_type
399  {
400  typedef DgemmNodeEvaluator<L, IndicesType, index> CWrappedEvaluator;
401  typedef typename CWrappedEvaluator::Evaluator CEvaluator;
402 
403  static const unsigned int nextIndex = index + L::TotalCount;
404  typedef DgemmNodeEvaluator<R, IndicesType, nextIndex> BetaWrappedEvaluator;
405  typedef typename BetaWrappedEvaluator::Evaluator BetaEvaluator;
406 
407  template<typename ArgumentVectorType>
408  static double GetBeta(const ArgumentVectorType& args)
409  {
410  return BetaEvaluator::Evaluate(args) * CWrappedEvaluator::GetScale();
411  }
412  };
413 
414  template<typename L, typename R, typename IndicesType, unsigned int index>
415  struct BetaCParameterAccessImpl< Node<L, expt::MultiplyOp, R>, IndicesType, index,
416  typename boost::enable_if
417  <
418  boost::mpl::and_
419  <
420  Nektar::CanGetRawPtr<typename R::ResultType>,
421  boost::is_same<typename L::ResultType, double>
422  >
423  >::type> : public boost::true_type
424  {
425  typedef DgemmNodeEvaluator<L, IndicesType, index> BetaWrappedEvaluator;
426  typedef typename BetaWrappedEvaluator::Evaluator BetaEvaluator;
427 
428  static const unsigned int nextIndex = index + L::TotalCount;
429  typedef DgemmNodeEvaluator<R, IndicesType, nextIndex> CWrappedEvaluator;
430  typedef typename CWrappedEvaluator::Evaluator CEvaluator;
431 
432  template<typename ArgumentVectorType>
433  static double GetBeta(const ArgumentVectorType& args)
434  {
435  return BetaEvaluator::Evaluate(args) * CWrappedEvaluator::GetScale();
436  }
437  };
438 
439  template<typename NodeType, typename IndicesType, unsigned int index>
440  struct BetaCParameterAccess : public BetaCParameterAccessImpl<NodeType, IndicesType, index> {};
441 
442  template<typename T, typename IndicesType, unsigned int index>
443  struct BetaCParameterAccess<Node<T, NegateOp, void>, IndicesType, index> : public BetaCParameterAccessImpl<T, IndicesType, index>
444  {
445  template<typename ArgumentVectorType>
446  static double GetBeta(const ArgumentVectorType& args)
447  {
448  return -BetaCParameterAccessImpl<T, IndicesType, index>::GetBeta(args);
449  }
450  };
451 
452  }
453 
454  // Three cases to handle.
455  // 1. alpha*A*B +/- beta*C
456  // 2. beta*C +/- alpha*A*B - Implement through commutative transform.
457  // 3. alpha*A*B +/- beta*C*D - This has to be considered separately because it matches both 1 and 2.
458  // 4. alpha*A*B
459 
460  ///////////////////////////////////////////////////////////
461  // Case 1: alpha*A*B +/- beta*C
462  // Case 3 as well.
463  ///////////////////////////////////////////////////////////
464  template<typename L, typename OpType, typename R, typename IndicesType, unsigned int index>
465  struct BinaryBinaryEvaluateNodeOverride<L, OpType, R, IndicesType, index,
466  typename boost::enable_if
467  <
468  boost::mpl::and_
469  <
470  IsAdditiveOperator<OpType>,
471  impl::AlphaABParameterAccess<L, IndicesType, index>,
472  impl::BetaCParameterAccess<R, IndicesType, index + L::TotalCount>
473  //IsAlphaABNode<L>,
474  //IsBetaCNode<R>
475  >
476  >::type> : public boost::true_type
477  {
478  typedef Node<L, OpType, R> NodeType;
479  typedef impl::AlphaABParameterAccess<L, IndicesType, index> LhsAccess;
480  static const unsigned int rhsIndex = index + L::TotalCount;
481  typedef impl::BetaCParameterAccess<R, IndicesType, rhsIndex> RhsAccess;
482  typedef typename LhsAccess::AEvaluator AEvaluator;
483  typedef typename LhsAccess::BEvaluator BEvaluator;
484  typedef typename RhsAccess::CEvaluator CEvaluator;
485 
486  template<typename ResultType, typename ArgumentVectorType>
487  static void Evaluate(ResultType& accumulator, const ArgumentVectorType& args)
488  {
489  typename AEvaluator::ResultType a = AEvaluator::Evaluate(args);
490  typename BEvaluator::ResultType b = BEvaluator::Evaluate(args);
491  typename CEvaluator::ResultType c = CEvaluator::Evaluate(args);
492 
493  double alpha = LhsAccess::GetAlpha(args);
494  double beta = RhsAccess::GetBeta(args);
495 
496  beta *= impl::BetaScale<OpType>::GetScale();
497  Nektar::Dgemm(accumulator, alpha, a, b, beta, c);
498  }
499  };
500 
501  ////////////////////////////////////////////////////////////
502  // Case 2 - beta*C +/- alpha*A*B
503  ////////////////////////////////////////////////////////////
504  template<typename L, typename OpType, typename R, typename IndicesType, unsigned int index>
505  struct BinaryBinaryEvaluateNodeOverride<L, OpType, R, IndicesType, index,
506  typename boost::enable_if
507  <
508  boost::mpl::and_
509  <
510  IsAdditiveOperator<OpType>,
511  impl::AlphaABParameterAccess<R, IndicesType, index + L::TotalCount>,
512  impl::BetaCParameterAccess<L, IndicesType, index>,
513  boost::mpl::not_<impl::AlphaABParameterAccess<L, IndicesType, index> >
514  //IsAlphaABNode<R>,
515  //IsBetaCNode<L>,
516  //boost::mpl::not_<IsAlphaABNode<L> >
517  >
518  >::type> : public boost::true_type
519  {
520  typedef Node<L, OpType, R> NodeType;
521 
522  typedef impl::BetaCParameterAccess<L, IndicesType, index> LhsAccess;
523  static const unsigned int rhsIndex = index + L::TotalCount;
524  typedef impl::AlphaABParameterAccess<R, IndicesType, rhsIndex> RhsAccess;
525 
526  typedef typename RhsAccess::AEvaluator AEvaluator;
527  typedef typename RhsAccess::BEvaluator BEvaluator;
528  typedef typename LhsAccess::CEvaluator CEvaluator;
529 
530  template<typename ResultType, typename ArgumentVectorType>
531  static void Evaluate(ResultType& accumulator, const ArgumentVectorType& args)
532  {
533  typename AEvaluator::ResultType a = AEvaluator::Evaluate(args);
534  typename BEvaluator::ResultType b = BEvaluator::Evaluate(args);
535  typename CEvaluator::ResultType c = CEvaluator::Evaluate(args);
536 
537  double alpha = RhsAccess::GetAlpha(args);
538  double beta = LhsAccess::GetBeta(args);
539 
540  alpha *= impl::BetaScale<OpType>::GetScale();
541  Nektar::Dgemm(accumulator, alpha, a, b, beta, c);
542  }
543  };
544 
545  ////////////////////////////////////////////////////////////
546  // Case 4 - alpha*A*B
547  ////////////////////////////////////////////////////////////
548  template<typename L, typename OpType, typename R, typename IndicesType, unsigned int index>
549  struct BinaryBinaryEvaluateNodeOverride<L, OpType, R, IndicesType, index,
550  typename boost::enable_if
551  <
552  //IsAlphaABNode<Node<L, OpType, R> >
553  impl::AlphaABParameterAccess<Node<L, OpType, R>, IndicesType, index>
554  >::type> : public boost::true_type
555  {
556  typedef Node<L, OpType, R> NodeType;
557  typedef impl::AlphaABParameterAccess<NodeType, IndicesType, index> LhsAccess;
558  typedef typename LhsAccess::AEvaluator AEvaluator;
559  typedef typename LhsAccess::BEvaluator BEvaluator;
560 
561  template<typename ResultType, typename ArgumentVectorType>
562  static void Evaluate(ResultType& accumulator, const ArgumentVectorType& args)
563  {
564  typename AEvaluator::ResultType a = AEvaluator::Evaluate(args);
565  typename BEvaluator::ResultType b = BEvaluator::Evaluate(args);
566 
567  double alpha = LhsAccess::GetAlpha(args);
568  Nektar::Dgemm(accumulator, alpha, a, b);
569  }
570  };
571 
572 
573 }
574 
575 #endif
576 #endif
Nektar::MultiplyEqual
void MultiplyEqual(NekMatrix< LhsDataType, StandardMatrixTag > &lhs, typename boost::call_traits< LhsDataType >::const_reference rhs)
Nektar::AddEqual
RhsMatrixType void AddEqual(NekMatrix< DataType, StandardMatrixTag > &result, const NekMatrix< RhsDataType, RhsMatrixType > &rhs)
Definition: MatrixOperations.cpp:169
Nektar::Multiply
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
Definition: MatrixOperations.cpp:42
CanGetRawPtr.hpp
L
NekDouble L
Definition: CompressibleBL.cpp:79
Blas.hpp
BlasOverrideUtil.hpp
Generated on Fri Dec 2 2016 13:10:46 for Nektar++ by   doxygen 1.8.8