Nektar++
MatrixVectorMultiplication.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MatrixVectorMultiplication.cpp
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/Blas.hpp>
48 #include <LibUtilities/LinearAlgebra/BlockMatrix.hpp>
49 #include <LibUtilities/LinearAlgebra/MatrixOperations.hpp>
50 #include <LibUtilities/LinearAlgebra/NekVectorFwd.hpp>
51 #include <LibUtilities/LinearAlgebra/ScaledMatrix.hpp>
52 #include <LibUtilities/LinearAlgebra/StandardMatrix.hpp>
53 
54 namespace Nektar
55 {
56 
57 template <typename DataType, typename LhsDataType, typename MatrixType>
58 void NekMultiplyUnspecializedMatrixType(
59  DataType *result, 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(
87  DataType *result, 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(
97  DataType *result, 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 template <typename DataType, typename LhsDataType>
120 void NekMultiplyBandedMatrix(
121  DataType *result, const NekMatrix<LhsDataType, BlockMatrixTag> &lhs,
122  const DataType *rhs,
123  typename std::enable_if<
124  !CanGetRawPtr<NekMatrix<LhsDataType, BlockMatrixTag>>::value>::type *p =
125  0)
126 {
127  boost::ignore_unused(result, lhs, rhs, p);
128 
129  NEKERROR(ErrorUtil::efatal,
130  "Banded block matrix multiplication not yet implemented");
131 }
132 
133 template <typename DataType, typename LhsInnerMatrixType>
134 void FullBlockMatrixMultiply(
135  NekVector<DataType> &result,
136  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag> &lhs,
137  const NekVector<DataType> &rhs)
138 {
139  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
140  unsigned int numberOfBlockColumns = lhs.GetNumberOfBlockColumns();
141  DataType *result_ptr = result.GetRawPtr();
142  const DataType *rhs_ptr = rhs.GetRawPtr();
143 
144  for (unsigned int i = 0; i < result.GetDimension(); ++i)
145  {
146  result_ptr[i] = DataType(0);
147  }
148  Array<OneD, DataType> temp(result.GetDimension());
149  DataType *temp_ptr = temp.get();
150 
151  unsigned int curResultRow = 0;
152  for (unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
153  {
154  unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
155 
156  if (blockRow != 0)
157  {
158  curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow - 1);
159  }
160 
161  if (rowsInBlock == 0)
162  {
163  continue;
164  }
165 
166  DataType *resultWrapper = result_ptr + curResultRow;
167 
168  unsigned int curWrapperRow = 0;
169  for (unsigned int blockColumn = 0; blockColumn < numberOfBlockColumns;
170  ++blockColumn)
171  {
172  if (blockColumn != 0)
173  {
174  curWrapperRow +=
175  lhs.GetNumberOfColumnsInBlockColumn(blockColumn - 1);
176  }
177 
178  // const std::shared_ptr<const LhsInnerMatrixType>& block =
179  // lhs.GetBlock(blockRow, blockColumn);
180  const LhsInnerMatrixType *block =
181  lhs.GetBlockPtr(blockRow, blockColumn);
182  if (!block)
183  {
184  continue;
185  }
186 
187  unsigned int columnsInBlock =
188  lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
189  if (columnsInBlock == 0)
190  {
191  continue;
192  }
193 
194  const DataType *rhsWrapper = rhs_ptr + curWrapperRow;
195  Multiply(temp_ptr, *block, rhsWrapper);
196  for (unsigned int i = 0; i < rowsInBlock; ++i)
197  {
198  resultWrapper[i] += temp_ptr[i];
199  }
200  }
201  }
202 }
203 
204 template <typename DataType, typename LhsInnerMatrixType>
205 void DiagonalBlockMatrixMultiply(
206  NekVector<DataType> &result,
207  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag> &lhs,
208  const NekVector<DataType> &rhs)
209 {
210  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
211  DataType *result_ptr = result.GetRawPtr();
212  const DataType *rhs_ptr = rhs.GetRawPtr();
213 
214  Array<OneD, unsigned int> rowSizes;
215  Array<OneD, unsigned int> colSizes;
216  lhs.GetBlockSizes(rowSizes, colSizes);
217 
218  unsigned int curResultRow = 0;
219  unsigned int curWrapperRow = 0;
220  unsigned int rowsInBlock = 0;
221  unsigned int columnsInBlock = 0;
222  for (unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
223  {
224  curResultRow += rowsInBlock;
225  curWrapperRow += columnsInBlock;
226  if (blockRow == 0)
227  {
228  rowsInBlock = rowSizes[blockRow] + 1;
229  columnsInBlock = colSizes[blockRow] + 1;
230  }
231  else
232  {
233  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow - 1];
234  columnsInBlock = colSizes[blockRow] - colSizes[blockRow - 1];
235  }
236 
237  if (rowsInBlock == 0)
238  {
239  continue;
240  }
241  if (columnsInBlock == 0)
242  {
243  std::fill(result.begin() + curResultRow,
244  result.begin() + curResultRow + rowsInBlock, 0.0);
245  continue;
246  }
247 
248  const LhsInnerMatrixType *block = lhs.GetBlockPtr(blockRow, blockRow);
249  if (!block)
250  {
251  continue;
252  }
253 
254  DataType *resultWrapper = result_ptr + curResultRow;
255  const DataType *rhsWrapper = rhs_ptr + curWrapperRow;
256  Multiply(resultWrapper, *block, rhsWrapper);
257  // resultWrapper = (*block)*rhsWrapper;
258  }
259  curResultRow += rowsInBlock;
260  if (curResultRow < result.GetRows())
261  {
262  std::fill(result.begin() + curResultRow, result.end(), 0.0);
263  }
264 }
265 
266 void DiagonalBlockFullScalMatrixMultiply(
267  NekVector<double> &result,
268  const NekMatrix<
269  NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>,
270  BlockMatrixTag> &lhs,
271  const NekVector<double> &rhs)
272 {
273  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
274  double *result_ptr = result.GetRawPtr();
275  const double *rhs_ptr = rhs.GetRawPtr();
276 
277  Array<OneD, unsigned int> rowSizes;
278  Array<OneD, unsigned int> colSizes;
279  lhs.GetBlockSizes(rowSizes, colSizes);
280 
281  unsigned int curResultRow = 0;
282  unsigned int curWrapperRow = 0;
283  unsigned int rowsInBlock = 0;
284  unsigned int columnsInBlock = 0;
285  for (unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
286  {
287  curResultRow += rowsInBlock;
288  curWrapperRow += columnsInBlock;
289  if (blockRow == 0)
290  {
291  rowsInBlock = rowSizes[blockRow] + 1;
292  columnsInBlock = colSizes[blockRow] + 1;
293  }
294  else
295  {
296  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow - 1];
297  columnsInBlock = colSizes[blockRow] - colSizes[blockRow - 1];
298  }
299 
300  if (rowsInBlock == 0)
301  {
302  continue;
303  }
304  if (columnsInBlock == 0)
305  {
306  std::fill(result.begin() + curResultRow,
307  result.begin() + curResultRow + rowsInBlock, 0.0);
308  continue;
309  }
310 
311  const DNekScalMat *block = lhs.GetBlockPtr(blockRow, blockRow);
312  if (!block)
313  {
314  continue;
315  }
316 
317  double *resultWrapper = result_ptr + curResultRow;
318  const double *rhsWrapper = rhs_ptr + curWrapperRow;
319 
320  // Multiply
321  const unsigned int *size = block->GetSize();
322  Blas::Gemv('N', size[0], size[1], block->Scale(), block->GetRawPtr(),
323  size[0], rhsWrapper, 1, 0.0, resultWrapper, 1);
324  }
325  curResultRow += rowsInBlock;
326  if (curResultRow < result.GetRows())
327  {
328  std::fill(result.begin() + curResultRow, result.end(), 0.0);
329  }
330 }
331 
332 void DiagonalBlockFullScalMatrixMultiply(
333  NekVector<NekSingle> &result,
334  const NekMatrix<
335  NekMatrix<NekMatrix<NekSingle, StandardMatrixTag>, ScaledMatrixTag>,
336  BlockMatrixTag> &lhs,
337  const NekVector<NekSingle> &rhs)
338 {
339  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
340  NekSingle *result_ptr = result.GetRawPtr();
341  const NekSingle *rhs_ptr = rhs.GetRawPtr();
342 
343  Array<OneD, unsigned int> rowSizes;
344  Array<OneD, unsigned int> colSizes;
345  lhs.GetBlockSizes(rowSizes, colSizes);
346 
347  unsigned int curResultRow = 0;
348  unsigned int curWrapperRow = 0;
349  unsigned int rowsInBlock = 0;
350  unsigned int columnsInBlock = 0;
351  for (unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
352  {
353  curResultRow += rowsInBlock;
354  curWrapperRow += columnsInBlock;
355  if (blockRow == 0)
356  {
357  rowsInBlock = rowSizes[blockRow] + 1;
358  columnsInBlock = colSizes[blockRow] + 1;
359  }
360  else
361  {
362  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow - 1];
363  columnsInBlock = colSizes[blockRow] - colSizes[blockRow - 1];
364  }
365 
366  if (rowsInBlock == 0)
367  {
368  continue;
369  }
370  if (columnsInBlock == 0)
371  {
372  std::fill(result.begin() + curResultRow,
373  result.begin() + curResultRow + rowsInBlock, 0.0);
374  continue;
375  }
376 
377  const SNekScalMat *block = lhs.GetBlockPtr(blockRow, blockRow);
378  if (!block)
379  {
380  continue;
381  }
382 
383  NekSingle *resultWrapper = result_ptr + curResultRow;
384  const NekSingle *rhsWrapper = rhs_ptr + curWrapperRow;
385 
386  // Multiply
387  const unsigned int *size = block->GetSize();
388  Blas::Gemv('N', size[0], size[1], block->Scale(), block->GetRawPtr(),
389  size[0], rhsWrapper, 1, 0.0, resultWrapper, 1);
390  }
391  curResultRow += rowsInBlock;
392  if (curResultRow < result.GetRows())
393  {
394  std::fill(result.begin() + curResultRow, result.end(), 0.0);
395  }
396 }
397 
398 template <typename DataType>
399 void NekMultiplyLowerTriangularMatrix(
400  DataType *result, const NekMatrix<DataType, StandardMatrixTag> &lhs,
401  const DataType *rhs)
402 {
403  int vectorSize = lhs.GetColumns();
404  std::copy(rhs, rhs + vectorSize, result);
405  int n = lhs.GetRows();
406  const DataType *a = lhs.GetRawPtr();
407  DataType *x = result;
408  int incx = 1;
409 
410  Blas::Tpmv('L', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
411 }
412 
413 template <typename DataType>
414 void NekMultiplyLowerTriangularMatrix(
415  DataType *result,
416  const NekMatrix<NekMatrix<DataType, StandardMatrixTag>, ScaledMatrixTag>
417  &lhs,
418  const DataType *rhs)
419 {
420  NekMultiplyLowerTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
421 
422  for (unsigned int i = 0; i < lhs.GetColumns(); ++i)
423  {
424  result[i] *= lhs.Scale();
425  }
426 }
427 
428 template <typename DataType, typename LhsDataType, typename MatrixType>
429 void NekMultiplyLowerTriangularMatrix(
430  DataType *result, const NekMatrix<LhsDataType, MatrixType> &lhs,
431  const DataType *rhs)
432 {
433  for (unsigned int i = 0; i < lhs.GetRows(); ++i)
434  {
435  DataType accum = DataType(0);
436  for (unsigned int j = 0; j <= i; ++j)
437  {
438  accum += lhs(i, j) * rhs[j];
439  }
440  result[i] = accum;
441  }
442 }
443 
444 template <typename DataType>
445 void NekMultiplyUpperTriangularMatrix(
446  DataType *result, const NekMatrix<DataType, StandardMatrixTag> &lhs,
447  const DataType *rhs)
448 {
449  int vectorSize = lhs.GetColumns();
450  std::copy(rhs, rhs + vectorSize, result);
451  int n = lhs.GetRows();
452  const DataType *a = lhs.GetRawPtr();
453  DataType *x = result;
454  int incx = 1;
455 
456  Blas::Tpmv('U', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
457 }
458 
459 template <typename DataType>
460 void NekMultiplyUpperTriangularMatrix(
461  DataType *result,
462  const NekMatrix<NekMatrix<DataType, StandardMatrixTag>, ScaledMatrixTag>
463  &lhs,
464  const DataType *rhs)
465 {
466  NekMultiplyUpperTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
467 
468  for (unsigned int i = 0; i < lhs.GetColumns(); ++i)
469  {
470  result[i] *= lhs.Scale();
471  }
472 }
473 
474 template <typename DataType, typename LhsDataType, typename MatrixType>
475 void NekMultiplyUpperTriangularMatrix(
476  DataType *result, const NekMatrix<LhsDataType, MatrixType> &lhs,
477  const DataType *rhs)
478 {
479  for (unsigned int i = 0; i < lhs.GetRows(); ++i)
480  {
481  DataType accum = DataType(0);
482  for (unsigned int j = i; j < lhs.GetColumns(); ++j)
483  {
484  accum += lhs(i, j) * rhs[j];
485  }
486  result[i] = accum;
487  }
488 }
489 
490 template <typename DataType, typename InnerMatrixType, typename MatrixTag>
491 void NekMultiplySymmetricMatrix(
492  DataType *result, const NekMatrix<InnerMatrixType, MatrixTag> &lhs,
493  const DataType *rhs,
494  typename std::enable_if<
495  CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag>>::value>::type *p =
496  0)
497 {
498  boost::ignore_unused(p);
499 
500  const unsigned int *size = lhs.GetSize();
501 
502  DataType alpha = lhs.Scale();
503  const DataType *a = lhs.GetRawPtr();
504  const DataType *x = rhs;
505  int incx = 1;
506  DataType beta = 0.0;
507  DataType *y = result;
508  int incy = 1;
509 
510  Blas::Spmv('U', size[0], alpha, a, x, incx, beta, y, incy);
511 }
512 
513 template <typename DataType, typename InnerMatrixType, typename MatrixTag>
514 void NekMultiplySymmetricMatrix(
515  DataType *result, const NekMatrix<InnerMatrixType, MatrixTag> &lhs,
516  const DataType *rhs,
517  typename std::enable_if<
518  !CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag>>::value>::type *p =
519  0)
520 {
521  boost::ignore_unused(p);
522 
523  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
524 }
525 
526 template <typename DataType, typename InnerMatrixType, typename MatrixTag>
527 void NekMultiplyFullMatrix(
528  DataType *result, const NekMatrix<InnerMatrixType, MatrixTag> &lhs,
529  const DataType *rhs,
530  typename std::enable_if<
531  CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag>>::value>::type *p =
532  0)
533 {
534  boost::ignore_unused(p);
535 
536  const unsigned int *size = lhs.GetSize();
537 
538  char t = lhs.GetTransposeFlag();
539 
540  DataType alpha = lhs.Scale();
541  const DataType *a = lhs.GetRawPtr();
542  int lda = size[0];
543  const DataType *x = rhs;
544  int incx = 1;
545  DataType beta = 0.0;
546  DataType *y = result;
547  int incy = 1;
548 
549  Blas::Gemv(t, size[0], size[1], alpha, a, lda, x, incx, beta, y, incy);
550 }
551 
552 template <typename DataType, typename InnerMatrixType, typename MatrixTag>
553 void NekMultiplyFullMatrix(
554  DataType *result, const NekMatrix<InnerMatrixType, MatrixTag> &lhs,
555  const DataType *rhs,
556  typename std::enable_if<
557  !CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag>>::value>::type *p =
558  0)
559 {
560  boost::ignore_unused(p);
561 
562  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
563 }
564 
565 template <typename DataType, typename LhsDataType, typename MatrixType>
566 void Multiply(DataType *result, const NekMatrix<LhsDataType, MatrixType> &lhs,
567  const DataType *rhs)
568 {
569  switch (lhs.GetType())
570  {
571  case eFULL:
572  NekMultiplyFullMatrix(result, lhs, rhs);
573  break;
574  case eDIAGONAL:
575  NekMultiplyDiagonalMatrix(result, lhs, rhs);
576  break;
577  case eUPPER_TRIANGULAR:
578  NekMultiplyUpperTriangularMatrix(result, lhs, rhs);
579  break;
580  case eLOWER_TRIANGULAR:
581  NekMultiplyLowerTriangularMatrix(result, lhs, rhs);
582  break;
583  case eSYMMETRIC:
584  NekMultiplySymmetricMatrix(result, lhs, rhs);
585  break;
586  case eBANDED:
587  NekMultiplyBandedMatrix(result, lhs, rhs);
588  break;
589  case eSYMMETRIC_BANDED:
590  case eUPPER_TRIANGULAR_BANDED:
591  case eLOWER_TRIANGULAR_BANDED:
592  default:
593  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
594  }
595 }
596 
597 template <typename DataType, typename LhsDataType, typename MatrixType>
598 void Multiply(NekVector<DataType> &result,
599  const NekMatrix<LhsDataType, MatrixType> &lhs,
600  const NekVector<DataType> &rhs)
601 {
602 
603  ASSERTL1(lhs.GetColumns() == rhs.GetRows(),
604  std::string("A left side matrix with column count ") +
605  std::to_string(lhs.GetColumns()) +
606  std::string(" and a right side vector with row count ") +
607  std::to_string(rhs.GetRows()) +
608  std::string(" can't be multiplied."));
609  Multiply(result.GetRawPtr(), lhs, rhs.GetRawPtr());
610 }
611 
612 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
613  Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES, (1, (void)),
614  (1, (NekVector<NekDouble> &)), (1, (const NekVector<NekDouble> &)))
615 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
616  Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES_SINGLE, (1, (void)),
617  (1, (NekVector<NekSingle> &)), (1, (const NekVector<NekSingle> &)))
618 
619 template <typename DataType, typename LhsInnerMatrixType>
620 void Multiply(NekVector<DataType> &result,
621  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag> &lhs,
622  const NekVector<DataType> &rhs)
623 {
624  if (lhs.GetStorageType() == eDIAGONAL)
625  {
626  DiagonalBlockMatrixMultiply(result, lhs, rhs);
627  }
628  else
629  {
630  FullBlockMatrixMultiply(result, lhs, rhs);
631  }
632 }
633 
634 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
635  Multiply, NEKTAR_BLOCK_MATRIX_TYPES, (1, (void)),
636  (1, (NekVector<NekDouble> &)), (1, (const NekVector<NekDouble> &)))
637 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
638  Multiply, NEKTAR_BLOCK_MATRIX_TYPES_SINGLE, (1, (void)),
639  (1, (NekVector<NekSingle> &)), (1, (const NekVector<NekSingle> &)))
640 
641 template <typename DataType, typename LhsDataType, typename MatrixType>
642 NekVector<DataType> Multiply(const NekMatrix<LhsDataType, MatrixType> &lhs,
643  const NekVector<DataType> &rhs)
644 {
645  NekVector<DataType> result(lhs.GetRows(), DataType(0));
646  Multiply(result, lhs, rhs);
647  return result;
648 }
649 
650 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
651  Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (NekVector<NekDouble>)), (0, ()),
652  (1, (const NekVector<NekDouble> &)))
653 NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(
654  Multiply, NEKTAR_ALL_MATRIX_TYPES_SINGLE, (1, (NekVector<NekSingle>)),
655  (0, ()), (1, (const NekVector<NekSingle> &)))
656 
657 } // namespace Nektar
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:249
NEKTAR_STANDARD_AND_SCALED_MATRICES
#define NEKTAR_STANDARD_AND_SCALED_MATRICES
Definition: ExplicitInstantiation.h:65
NEKTAR_ALL_MATRIX_TYPES
#define NEKTAR_ALL_MATRIX_TYPES
Definition: ExplicitInstantiation.h:58
NEKTAR_BLOCK_MATRIX_TYPES
#define NEKTAR_BLOCK_MATRIX_TYPES
Definition: ExplicitInstantiation.h:62
Nektar::ConstMatrix::GetRows
unsigned int GetRows() const
Definition: MatrixBase.cpp:65
Nektar::ConstMatrix::GetTransposeFlag
char GetTransposeFlag() const
Definition: MatrixBase.cpp:122
Nektar::ConstMatrix::GetColumns
unsigned int GetColumns() const
Definition: MatrixBase.cpp:84
Nektar::NekVector::GetDimension
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:201
Nektar::NekVector::GetRows
unsigned int GetRows() const
Definition: NekVector.cpp:206
Nektar::NekVector::begin
iterator begin()
Definition: NekVector.cpp:235
Nektar::NekVector::end
iterator end()
Definition: NekVector.cpp:241
Nektar::NekVector::GetRawPtr
DataType * GetRawPtr()
Definition: NekVector.cpp:211
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:228
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:254
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:281
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:304
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
CellMLToNektar.pycml.copy
def copy(self)
Definition: pycml.py:2663
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
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::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::DiagonalBlockMatrixMultiply
void DiagonalBlockMatrixMultiply(NekVector< DataType > &result, const NekMatrix< LhsInnerMatrixType, BlockMatrixTag > &lhs, const NekVector< DataType > &rhs)
Definition: MatrixVectorMultiplication.cpp:205
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:399
Nektar::NEKTAR_ALL_MATRIX_TYPES_SINGLE
NEKTAR_ALL_MATRIX_TYPES_SINGLE
Definition: MatrixOperations.cpp:61
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:491
Nektar::NekMultiplyDiagonalMatrix
void NekMultiplyDiagonalMatrix(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:74
Nektar::DiagonalBlockFullScalMatrixMultiply
void DiagonalBlockFullScalMatrixMultiply(NekVector< double > &result, const NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > &lhs, const NekVector< double > &rhs)
Definition: MatrixVectorMultiplication.cpp:266
Nektar::NekMultiplyUnspecializedMatrixType
void NekMultiplyUnspecializedMatrixType(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:58
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:527
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:134
Nektar::NekMultiplyUpperTriangularMatrix
void NekMultiplyUpperTriangularMatrix(DataType *result, const NekMatrix< DataType, StandardMatrixTag > &lhs, const DataType *rhs)
Definition: MatrixVectorMultiplication.cpp:445
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:209
Generated on Sat Mar 25 2023 12:35:51 for Nektar++ by doxygen 1.9.1