Nektar++
MatrixVectorMultiplication.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MatrixVectorMultiplication.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description:
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <type_traits>
36 
37 #include <boost/core/ignore_unused.hpp>
38 
39 #include <LibUtilities/LinearAlgebra/ExplicitInstantiation.h>
40 #include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
41 
42 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
43 
44 #include <LibUtilities/LinearAlgebra/Blas.hpp>
45 #include <LibUtilities/LinearAlgebra/CanGetRawPtr.hpp>
46 
47 #include <LibUtilities/LinearAlgebra/StandardMatrix.hpp>
48 #include <LibUtilities/LinearAlgebra/BlockMatrix.hpp>
49 #include <LibUtilities/LinearAlgebra/ScaledMatrix.hpp>
50 #include <LibUtilities/LinearAlgebra/Blas.hpp>
51 #include <LibUtilities/LinearAlgebra/NekVectorFwd.hpp>
52 #include <LibUtilities/LinearAlgebra/MatrixOperations.hpp>
53 
54 namespace Nektar
55 {
56 
57  template<typename DataType, typename LhsDataType, typename MatrixType>
58  void NekMultiplyUnspecializedMatrixType(DataType* result,
59  const NekMatrix<LhsDataType, MatrixType>& lhs,
60  const DataType* rhs)
61  {
62  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
63  {
64  DataType accum = DataType(0);
65  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
66  {
67  accum += lhs(i,j)*rhs[j];
68  }
69  result[i] = accum;
70  }
71  }
72 
73  template<typename DataType, typename LhsDataType, typename MatrixType>
74  void NekMultiplyDiagonalMatrix(DataType* result,
75  const NekMatrix<LhsDataType, MatrixType>& lhs,
76  const DataType* rhs)
77  {
78  int n = lhs.GetRows();
79  for(unsigned int i = 0; i < n; ++i)
80  {
81  result[i] = lhs(i,i)*rhs[i];
82  }
83  }
84 
85  template<typename DataType, typename LhsDataType>
86  void NekMultiplyDiagonalMatrix(DataType* result,
87  const NekMatrix<LhsDataType, StandardMatrixTag>& lhs,
88  const DataType* rhs)
89  {
90  int n = lhs.GetRows();
91  const DataType* mat_ptr = lhs.GetRawPtr();
92  Vmath::Vmul(n, mat_ptr, 1, rhs, 1, result, 1);
93  }
94 
95  template<typename DataType, typename LhsDataType, typename MatrixType>
96  void NekMultiplyBandedMatrix(DataType* result,
97  const NekMatrix<LhsDataType, MatrixType>& lhs,
98  const DataType* rhs,
99  typename std::enable_if<
100  CanGetRawPtr<NekMatrix<LhsDataType, MatrixType> >::value >::type* p=0)
101  {
102  boost::ignore_unused(p);
103 
104  int m = lhs.GetRows();
105  int n = lhs.GetColumns();
106  int kl = lhs.GetNumberOfSubDiagonals();
107  int ku = lhs.GetNumberOfSuperDiagonals();
108  DataType alpha = lhs.Scale();
109  const DataType* a = lhs.GetRawPtr();
110  int lda = kl + ku + 1;
111  const DataType* x = rhs;
112  int incx = 1;
113  DataType beta = 0.0;
114  DataType* y = result;
115  int incy = 1;
116  Blas::Gbmv('N', m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy);
117 
118  }
119 
120  template<typename DataType, typename LhsDataType>
121  void NekMultiplyBandedMatrix(DataType* result,
122  const NekMatrix<LhsDataType, BlockMatrixTag>& lhs,
123  const DataType* rhs,
124  typename std::enable_if<
125  !CanGetRawPtr<NekMatrix<LhsDataType, BlockMatrixTag> >::value>::type* p=0)
126  {
127  boost::ignore_unused(result, lhs, rhs, p);
128 
129  NEKERROR(ErrorUtil::efatal, "Banded block matrix multiplication not yet implemented");
130  }
131 
132  template<typename DataType, typename LhsInnerMatrixType>
133  void FullBlockMatrixMultiply(NekVector<DataType>& result,
134  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
135  const NekVector<DataType>& rhs)
136  {
137  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
138  unsigned int numberOfBlockColumns = lhs.GetNumberOfBlockColumns();
139  DataType* result_ptr = result.GetRawPtr();
140  const DataType* rhs_ptr = rhs.GetRawPtr();
141 
142  for(unsigned int i = 0; i < result.GetDimension(); ++i)
143  {
144  result_ptr[i] = DataType(0);
145  }
146  Array<OneD, DataType> temp(result.GetDimension());
147  DataType* temp_ptr = temp.get();
148 
149  unsigned int curResultRow = 0;
150  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
151  {
152  unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
153 
154  if( blockRow != 0 )
155  {
156  curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow-1);
157  }
158 
159  if( rowsInBlock == 0 )
160  {
161  continue;
162  }
163 
164  DataType* resultWrapper = result_ptr + curResultRow;
165 
166  unsigned int curWrapperRow = 0;
167  for(unsigned int blockColumn = 0; blockColumn < numberOfBlockColumns; ++blockColumn)
168  {
169  if( blockColumn != 0 )
170  {
171  curWrapperRow += lhs.GetNumberOfColumnsInBlockColumn(blockColumn-1);
172  }
173 
174  //const std::shared_ptr<const LhsInnerMatrixType>& block = lhs.GetBlock(blockRow, blockColumn);
175  const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockColumn);
176  if( !block )
177  {
178  continue;
179  }
180 
181  unsigned int columnsInBlock = lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
182  if( columnsInBlock == 0 )
183  {
184  continue;
185  }
186 
187  const DataType* rhsWrapper = rhs_ptr + curWrapperRow;
188  Multiply(temp_ptr, *block, rhsWrapper);
189  for(unsigned int i = 0; i < rowsInBlock; ++i)
190  {
191  resultWrapper[i] += temp_ptr[i];
192  }
193  }
194  }
195  }
196 
197  template<typename DataType, typename LhsInnerMatrixType>
198  void DiagonalBlockMatrixMultiply(NekVector<DataType>& result,
199  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
200  const NekVector<DataType>& rhs)
201  {
202  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
203  DataType* result_ptr = result.GetRawPtr();
204  const DataType* rhs_ptr = rhs.GetRawPtr();
205 
206  Array<OneD, unsigned int> rowSizes;
207  Array<OneD, unsigned int> colSizes;
208  lhs.GetBlockSizes(rowSizes, colSizes);
209 
210  unsigned int curResultRow = 0;
211  unsigned int curWrapperRow = 0;
212  unsigned int rowsInBlock = 0;
213  unsigned int columnsInBlock = 0;
214  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
215  {
216  curResultRow += rowsInBlock;
217  curWrapperRow += columnsInBlock;
218  if ( blockRow == 0)
219  {
220  rowsInBlock = rowSizes[blockRow] + 1;
221  columnsInBlock = colSizes[blockRow] + 1;
222  }
223  else
224  {
225  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
226  columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
227  }
228 
229  if( rowsInBlock == 0)
230  {
231  continue;
232  }
233  if( columnsInBlock == 0)
234  {
235  std::fill(result.begin()+curResultRow,
236  result.begin()+curResultRow + rowsInBlock, 0.0);
237  continue;
238  }
239 
240  const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockRow);
241  if( !block )
242  {
243  continue;
244  }
245 
246  DataType* resultWrapper = result_ptr + curResultRow;
247  const DataType* rhsWrapper = rhs_ptr + curWrapperRow;
248  Multiply(resultWrapper, *block, rhsWrapper);
249  //resultWrapper = (*block)*rhsWrapper;
250  }
251  curResultRow += rowsInBlock;
252  if (curResultRow < result.GetRows())
253  {
254  std::fill(result.begin()+curResultRow, result.end(), 0.0);
255  }
256  }
257 
258  void DiagonalBlockFullScalMatrixMultiply(NekVector<double>& result,
259  const NekMatrix<NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>, BlockMatrixTag>& lhs,
260  const NekVector<double>& rhs)
261  {
262  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
263  double* result_ptr = result.GetRawPtr();
264  const double* rhs_ptr = rhs.GetRawPtr();
265 
266  Array<OneD, unsigned int> rowSizes;
267  Array<OneD, unsigned int> colSizes;
268  lhs.GetBlockSizes(rowSizes, colSizes);
269 
270  unsigned int curResultRow = 0;
271  unsigned int curWrapperRow = 0;
272  unsigned int rowsInBlock = 0;
273  unsigned int columnsInBlock = 0;
274  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
275  {
276  curResultRow += rowsInBlock;
277  curWrapperRow += columnsInBlock;
278  if ( blockRow == 0)
279  {
280  rowsInBlock = rowSizes[blockRow] + 1;
281  columnsInBlock = colSizes[blockRow] + 1;
282  }
283  else
284  {
285  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
286  columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
287  }
288 
289  if( rowsInBlock == 0)
290  {
291  continue;
292  }
293  if( columnsInBlock == 0)
294  {
295  std::fill(result.begin()+curResultRow,
296  result.begin()+curResultRow + rowsInBlock, 0.0);
297  continue;
298  }
299 
300  const DNekScalMat* block = lhs.GetBlockPtr(blockRow, blockRow);
301  if( !block )
302  {
303  continue;
304  }
305 
306  double* resultWrapper = result_ptr + curResultRow;
307  const double* rhsWrapper = rhs_ptr + curWrapperRow;
308 
309  // Multiply
310  const unsigned int* size = block->GetSize();
311  Blas::Gemv('N', size[0], size[1], block->Scale(),
312  block->GetRawPtr(), size[0], rhsWrapper, 1,
313  0.0, resultWrapper, 1);
314  }
315  curResultRow += rowsInBlock;
316  if (curResultRow < result.GetRows())
317  {
318  std::fill(result.begin()+curResultRow, result.end(), 0.0);
319  }
320  }
321 
322  void DiagonalBlockFullScalMatrixMultiply(NekVector<NekSingle>& result,
323  const NekMatrix<NekMatrix<NekMatrix<NekSingle, StandardMatrixTag>, ScaledMatrixTag>, BlockMatrixTag>& lhs,
324  const NekVector<NekSingle>& rhs)
325  {
326  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
327  NekSingle* result_ptr = result.GetRawPtr();
328  const NekSingle* rhs_ptr = rhs.GetRawPtr();
329 
330  Array<OneD, unsigned int> rowSizes;
331  Array<OneD, unsigned int> colSizes;
332  lhs.GetBlockSizes(rowSizes, colSizes);
333 
334  unsigned int curResultRow = 0;
335  unsigned int curWrapperRow = 0;
336  unsigned int rowsInBlock = 0;
337  unsigned int columnsInBlock = 0;
338  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
339  {
340  curResultRow += rowsInBlock;
341  curWrapperRow += columnsInBlock;
342  if ( blockRow == 0)
343  {
344  rowsInBlock = rowSizes[blockRow] + 1;
345  columnsInBlock = colSizes[blockRow] + 1;
346  }
347  else
348  {
349  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
350  columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
351  }
352 
353  if( rowsInBlock == 0)
354  {
355  continue;
356  }
357  if( columnsInBlock == 0)
358  {
359  std::fill(result.begin()+curResultRow,
360  result.begin()+curResultRow + rowsInBlock, 0.0);
361  continue;
362  }
363 
364  const SNekScalMat* block = lhs.GetBlockPtr(blockRow, blockRow);
365  if( !block )
366  {
367  continue;
368  }
369 
370  NekSingle* resultWrapper = result_ptr + curResultRow;
371  const NekSingle* rhsWrapper = rhs_ptr + curWrapperRow;
372 
373  // Multiply
374  const unsigned int* size = block->GetSize();
375  Blas::Gemv('N', size[0], size[1], block->Scale(),
376  block->GetRawPtr(), size[0], rhsWrapper, 1,
377  0.0, resultWrapper, 1);
378  }
379  curResultRow += rowsInBlock;
380  if (curResultRow < result.GetRows())
381  {
382  std::fill(result.begin()+curResultRow, result.end(), 0.0);
383  }
384  }
385 
386  template<typename DataType>
387  void NekMultiplyLowerTriangularMatrix(DataType* result,
388  const NekMatrix<DataType, StandardMatrixTag>& lhs,
389  const DataType* rhs)
390  {
391  int vectorSize = lhs.GetColumns();
392  std::copy(rhs, rhs+vectorSize, result);
393  int n = lhs.GetRows();
394  const DataType* a = lhs.GetRawPtr();
395  DataType* x = result;
396  int incx = 1;
397 
398  Blas::Tpmv('L', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
399  }
400 
401  template<typename DataType>
402  void NekMultiplyLowerTriangularMatrix(DataType* result,
403  const NekMatrix<NekMatrix<DataType, StandardMatrixTag>, ScaledMatrixTag>& lhs,
404  const DataType* rhs)
405  {
406  NekMultiplyLowerTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
407 
408  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
409  {
410  result[i] *= lhs.Scale();
411  }
412  }
413 
414  template<typename DataType, typename LhsDataType, typename MatrixType>
415  void NekMultiplyLowerTriangularMatrix(DataType* result,
416  const NekMatrix<LhsDataType, MatrixType>& lhs,
417  const DataType* rhs)
418  {
419  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
420  {
421  DataType accum = DataType(0);
422  for(unsigned int j = 0; j <= i; ++j)
423  {
424  accum += lhs(i,j)*rhs[j];
425  }
426  result[i] = accum;
427  }
428  }
429 
430  template<typename DataType>
431  void NekMultiplyUpperTriangularMatrix(DataType* result,
432  const NekMatrix<DataType, StandardMatrixTag>& lhs,
433  const DataType* rhs)
434  {
435  int vectorSize = lhs.GetColumns();
436  std::copy(rhs, rhs+vectorSize, result);
437  int n = lhs.GetRows();
438  const DataType* a = lhs.GetRawPtr();
439  DataType* x = result;
440  int incx = 1;
441 
442  Blas::Tpmv('U', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
443  }
444 
445  template<typename DataType>
446  void NekMultiplyUpperTriangularMatrix(DataType* result,
447  const NekMatrix<NekMatrix<DataType, StandardMatrixTag>, ScaledMatrixTag>& lhs,
448  const DataType* rhs)
449  {
450  NekMultiplyUpperTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
451 
452  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
453  {
454  result[i] *= lhs.Scale();
455  }
456  }
457 
458  template<typename DataType, typename LhsDataType, typename MatrixType>
459  void NekMultiplyUpperTriangularMatrix(DataType* result,
460  const NekMatrix<LhsDataType, MatrixType>& lhs,
461  const DataType* rhs)
462  {
463  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
464  {
465  DataType accum = DataType(0);
466  for(unsigned int j = i; j < lhs.GetColumns(); ++j)
467  {
468  accum += lhs(i,j)*rhs[j];
469  }
470  result[i] = accum;
471  }
472  }
473 
474  template<typename DataType, typename InnerMatrixType, typename MatrixTag>
475  void NekMultiplySymmetricMatrix(DataType* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const DataType* rhs,
476  typename std::enable_if<CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> >::value >::type* p=0)
477  {
478  boost::ignore_unused(p);
479 
480  const unsigned int* size = lhs.GetSize();
481 
482  DataType alpha = lhs.Scale();
483  const DataType* a = lhs.GetRawPtr();
484  const DataType* x = rhs;
485  int incx = 1;
486  DataType beta = 0.0;
487  DataType* y = result;
488  int incy = 1;
489 
490  Blas::Spmv('U', size[0], alpha, a, x, incx, beta, y, incy);
491  }
492 
493  template<typename DataType, typename InnerMatrixType, typename MatrixTag>
494  void NekMultiplySymmetricMatrix(DataType* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const DataType* rhs,
495  typename std::enable_if<!CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> >::value >::type* p = 0)
496  {
497  boost::ignore_unused(p);
498 
499  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
500  }
501 
502  template<typename DataType, typename InnerMatrixType, typename MatrixTag>
503  void NekMultiplyFullMatrix(DataType* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const DataType* rhs,
504  typename std::enable_if<CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> >::value >::type* p=0)
505  {
506  boost::ignore_unused(p);
507 
508  const unsigned int* size = lhs.GetSize();
509 
510  char t = lhs.GetTransposeFlag();
511 
512  DataType alpha = lhs.Scale();
513  const DataType* a = lhs.GetRawPtr();
514  int lda = size[0];
515  const DataType* x = rhs;
516  int incx = 1;
517  DataType beta = 0.0;
518  DataType* y = result;
519  int incy = 1;
520 
521  Blas::Gemv(t, size[0], size[1], alpha, a, lda, x, incx, beta, y, incy);
522  }
523 
524  template<typename DataType, typename InnerMatrixType, typename MatrixTag>
525  void NekMultiplyFullMatrix(DataType* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const DataType* rhs,
526  typename std::enable_if<!CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> >::value>::type* p = 0)
527  {
528  boost::ignore_unused(p);
529 
530  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
531  }
532 
533  template<typename DataType, typename LhsDataType, typename MatrixType>
534  void Multiply(DataType* result,
535  const NekMatrix<LhsDataType, MatrixType>& lhs,
536  const DataType* rhs)
537  {
538  switch(lhs.GetType())
539  {
540  case eFULL:
541  NekMultiplyFullMatrix(result, lhs, rhs);
542  break;
543  case eDIAGONAL:
544  NekMultiplyDiagonalMatrix(result, lhs, rhs);
545  break;
546  case eUPPER_TRIANGULAR:
547  NekMultiplyUpperTriangularMatrix(result, lhs, rhs);
548  break;
549  case eLOWER_TRIANGULAR:
550  NekMultiplyLowerTriangularMatrix(result, lhs, rhs);
551  break;
552  case eSYMMETRIC:
553  NekMultiplySymmetricMatrix(result, lhs, rhs);
554  break;
555  case eBANDED:
556  NekMultiplyBandedMatrix(result, lhs, rhs);
557  break;
558  case eSYMMETRIC_BANDED:
559  case eUPPER_TRIANGULAR_BANDED:
560  case eLOWER_TRIANGULAR_BANDED:
561  default:
562  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
563  }
564  }
565 
566  template<typename DataType, typename LhsDataType, typename MatrixType>
567  void Multiply(NekVector<DataType>& result,
568  const NekMatrix<LhsDataType, MatrixType>& lhs,
569  const NekVector<DataType>& rhs)
570  {
571 
572  ASSERTL1(lhs.GetColumns() == rhs.GetRows(), std::string("A left side matrix with column count ") +
573  std::to_string(lhs.GetColumns()) +
574  std::string(" and a right side vector with row count ") +
575  std::to_string(rhs.GetRows()) + std::string(" can't be multiplied."));
576  Multiply(result.GetRawPtr(), lhs, rhs.GetRawPtr());
577  }
578 
579  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)))
580  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES_SINGLE, (1, (void)), (1, (NekVector<NekSingle>&)), (1,(const NekVector<NekSingle>&)))
581 
582  template<typename DataType, typename LhsInnerMatrixType>
583  void Multiply(NekVector<DataType>& result,
584  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
585  const NekVector<DataType>& rhs)
586  {
587  if( lhs.GetStorageType() == eDIAGONAL )
588  {
589  DiagonalBlockMatrixMultiply(result, lhs, rhs);
590  }
591  else
592  {
593  FullBlockMatrixMultiply(result, lhs, rhs);
594  }
595  }
596 
597  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_BLOCK_MATRIX_TYPES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)))
598  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_BLOCK_MATRIX_TYPES_SINGLE, (1, (void)), (1, (NekVector<NekSingle>&)), (1,(const NekVector<NekSingle>&)))
599 
600  template<typename DataType, typename LhsDataType, typename MatrixType>
601  NekVector<DataType>
602  Multiply(const NekMatrix<LhsDataType, MatrixType>& lhs,
603  const NekVector<DataType>& rhs)
604  {
605  NekVector<DataType> result(lhs.GetRows(), DataType(0));
606  Multiply(result, lhs, rhs);
607  return result;
608  }
609 
610  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (NekVector<NekDouble>)), (0, ()), (1,(const NekVector<NekDouble>&)))
611  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES_SINGLE, (1, (NekVector<NekSingle>)), (0, ()), (1,(const NekVector<NekSingle>&)))
612 
613 }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
NEKTAR_STANDARD_AND_SCALED_MATRICES
#define NEKTAR_STANDARD_AND_SCALED_MATRICES
Definition: ExplicitInstantiation.h:58
NEKTAR_ALL_MATRIX_TYPES
#define NEKTAR_ALL_MATRIX_TYPES
Definition: ExplicitInstantiation.h:56
NEKTAR_BLOCK_MATRIX_TYPES
#define NEKTAR_BLOCK_MATRIX_TYPES
Definition: ExplicitInstantiation.h:57
Nektar::ConstMatrix::GetRows
unsigned int GetRows() const
Definition: MatrixBase.cpp:58
Nektar::ConstMatrix::GetTransposeFlag
char GetTransposeFlag() const
Definition: MatrixBase.cpp:113
Nektar::ConstMatrix::GetColumns
unsigned int GetColumns() const
Definition: MatrixBase.cpp:77
Nektar::NekVector::GetDimension
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:209
Nektar::NekVector::GetRows
unsigned int GetRows() const
Definition: NekVector.cpp:215
Nektar::NekVector::begin
iterator begin()
Definition: NekVector.cpp:239
Nektar::NekVector::end
iterator end()
Definition: NekVector.cpp:242
Nektar::NekVector::GetRawPtr
DataType * GetRawPtr()
Definition: NekVector.cpp:221
Blas::Gemv
static void Gemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:247
Blas::Gbmv
static void Gbmv(const char &trans, const int &m, const int &n, const int &kl, const int &ku, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
Definition: Blas.hpp:273
Blas::Tpmv
static void Tpmv(const char &uplo, const char &trans, const char &diag, const int &n, const double *ap, double *x, const int &incx)
Definition: Blas.hpp:305
Blas::Spmv
static void Spmv(const char &uplo, const int &n, const double &alpha, const double *a, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A is symmetric packed.
Definition: Blas.hpp:327
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
CellMLToNektar.pycml.copy
def copy(self)
Definition: pycml.py:2663
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX
SNekMat NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(AddEqualNegatedLhs, NEKTAR_ALL_MATRIX_TYPES,(1,(void)),(1,(DNekMat &)),(0,())) NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(AddEqualNegatedLhs
Nektar::NekMultiplyFullMatrix
void NekMultiplyFullMatrix(DataType *result, const NekMatrix< InnerMatrixType, MatrixTag > &lhs, const DataType *rhs, typename std::enable_if< CanGetRawPtr< NekMatrix< InnerMatrixType, MatrixTag > >::value >::type *p=0)
Definition: MatrixVectorMultiplication.cpp:503
Nektar::DiagonalBlockMatrixMultiply
void DiagonalBlockMatrixMultiply(NekVector< DataType > &result, const NekMatrix< LhsInnerMatrixType, BlockMatrixTag > &lhs, const NekVector< DataType > &rhs)
Definition: MatrixVectorMultiplication.cpp:198
Nektar::Multiply
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
Definition: MatrixOperations.cpp:41
Nektar::NekMultiplyLowerTriangularMatrix
void NekMultiplyLowerTriangularMatrix(DataType *result, const NekMatrix< DataType, StandardMatrixTag > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:387
Nektar::NEKTAR_ALL_MATRIX_TYPES_SINGLE
NEKTAR_ALL_MATRIX_TYPES_SINGLE
Definition: MatrixOperations.cpp:59
Nektar::NekMultiplyDiagonalMatrix
void NekMultiplyDiagonalMatrix(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:74
Nektar::NekMultiplySymmetricMatrix
void NekMultiplySymmetricMatrix(DataType *result, const NekMatrix< InnerMatrixType, MatrixTag > &lhs, const DataType *rhs, typename std::enable_if< CanGetRawPtr< NekMatrix< InnerMatrixType, MatrixTag > >::value >::type *p=0)
Definition: MatrixVectorMultiplication.cpp:475
Nektar::DiagonalBlockFullScalMatrixMultiply
void DiagonalBlockFullScalMatrixMultiply(NekVector< double > &result, const NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > &lhs, const NekVector< double > &rhs)
Definition: MatrixVectorMultiplication.cpp:258
Nektar::NekMultiplyBandedMatrix
void NekMultiplyBandedMatrix(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs, typename std::enable_if< CanGetRawPtr< NekMatrix< LhsDataType, MatrixType > >::value >::type *p=0)
Definition: MatrixVectorMultiplication.cpp:96
Nektar::NekMultiplyUnspecializedMatrixType
void NekMultiplyUnspecializedMatrixType(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:58
Nektar::eLOWER_TRIANGULAR_BANDED
@ eLOWER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:53
Nektar::eUPPER_TRIANGULAR_BANDED
@ eUPPER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:52
Nektar::FullBlockMatrixMultiply
void FullBlockMatrixMultiply(NekVector< DataType > &result, const NekMatrix< LhsInnerMatrixType, BlockMatrixTag > &lhs, const NekVector< DataType > &rhs)
Definition: MatrixVectorMultiplication.cpp:133
Nektar::NekMultiplyUpperTriangularMatrix
void NekMultiplyUpperTriangularMatrix(DataType *result, const NekMatrix< DataType, StandardMatrixTag > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:431
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
Generated on Wed Nov 24 2021 10:33:27 for Nektar++ by doxygen 1.9.1