Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description:
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
39 
42 
45 
46 #include <ExpressionTemplates/CommutativeTraits.hpp>
47 #include <boost/type_traits.hpp>
48 
56 
57 
58 
59 namespace Nektar
60 {
61 
62  template<typename DataType, typename LhsDataType, typename MatrixType>
63  void NekMultiplyUnspecializedMatrixType(DataType* result,
65  const DataType* rhs)
66  {
67  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
68  {
69  DataType accum = DataType(0);
70  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
71  {
72  accum += lhs(i,j)*rhs[j];
73  }
74  result[i] = accum;
75  }
76  }
77 
78  template<typename DataType, typename LhsDataType, typename MatrixType>
79  void NekMultiplyDiagonalMatrix(DataType* result,
81  const DataType* rhs)
82  {
83  int n = lhs.GetRows();
84  for(unsigned int i = 0; i < n; ++i)
85  {
86  result[i] = lhs(i,i)*rhs[i];
87  }
88  }
89 
90  template<typename DataType, typename LhsDataType>
91  void NekMultiplyDiagonalMatrix(DataType* result,
93  const DataType* rhs)
94  {
95  int n = lhs.GetRows();
96  const DataType* mat_ptr = lhs.GetRawPtr();
97  Vmath::Vmul(n, mat_ptr, 1, rhs, 1, result, 1);
98  }
99 
100  template<typename LhsDataType, typename MatrixType>
103  const NekDouble* rhs,
104  typename boost::enable_if
105  <
107  >::type* p=0)
108  {
109  int m = lhs.GetRows();
110  int n = lhs.GetColumns();
111  int kl = lhs.GetNumberOfSubDiagonals();
112  int ku = lhs.GetNumberOfSuperDiagonals();
113  double alpha = lhs.Scale();
114  const double* a = lhs.GetRawPtr();
115  int lda = kl + ku + 1;
116  const double* x = rhs;
117  int incx = 1;
118  double beta = 0.0;
119  double* y = result;
120  int incy = 1;
121  Blas::Dgbmv('N', m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy);
122 
123  }
124 
125  template<typename DataType, typename LhsDataType>
126  void NekMultiplyBandedMatrix(DataType* result,
128  const DataType* rhs,
129  typename boost::enable_if
130  <
132  >::type* p=0)
133  {
134  NEKERROR(ErrorUtil::efatal, "Banded block matrix multiplication not yet implemented");
135  }
136 
137  template<typename DataType, typename LhsInnerMatrixType>
140  const NekVector<DataType>& rhs)
141  {
142  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
143  unsigned int numberOfBlockColumns = lhs.GetNumberOfBlockColumns();
144  DataType* result_ptr = result.GetRawPtr();
145  const DataType* rhs_ptr = rhs.GetRawPtr();
146 
147  for(unsigned int i = 0; i < result.GetDimension(); ++i)
148  {
149  result_ptr[i] = DataType(0);
150  }
151  Array<OneD, DataType> temp(result.GetDimension());
152  DataType* temp_ptr = temp.get();
153 
154  unsigned int curResultRow = 0;
155  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
156  {
157  unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
158 
159  if( blockRow != 0 )
160  {
161  curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow-1);
162  }
163 
164  if( rowsInBlock == 0 )
165  {
166  continue;
167  }
168 
169  DataType* resultWrapper = result_ptr + curResultRow;
170 
171  unsigned int curWrapperRow = 0;
172  for(unsigned int blockColumn = 0; blockColumn < numberOfBlockColumns; ++blockColumn)
173  {
174  if( blockColumn != 0 )
175  {
176  curWrapperRow += lhs.GetNumberOfColumnsInBlockColumn(blockColumn-1);
177  }
178 
179  //const boost::shared_ptr<const LhsInnerMatrixType>& block = lhs.GetBlock(blockRow, blockColumn);
180  const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockColumn);
181  if( !block )
182  {
183  continue;
184  }
185 
186  unsigned int columnsInBlock = lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
187  if( columnsInBlock == 0 )
188  {
189  continue;
190  }
191 
192  const DataType* rhsWrapper = rhs_ptr + curWrapperRow;
193  Multiply(temp_ptr, *block, rhsWrapper);
194  for(unsigned int i = 0; i < rowsInBlock; ++i)
195  {
196  resultWrapper[i] += temp_ptr[i];
197  }
198  }
199  }
200  }
201 
202  template<typename LhsInnerMatrixType>
205  const NekVector<double>& rhs)
206  {
207  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
208  double* result_ptr = result.GetRawPtr();
209  const double* rhs_ptr = rhs.GetRawPtr();
210 
211  Array<OneD, unsigned int> rowSizes;
212  Array<OneD, unsigned int> colSizes;
213  lhs.GetBlockSizes(rowSizes, colSizes);
214 
215  unsigned int curResultRow = 0;
216  unsigned int curWrapperRow = 0;
217  unsigned int rowsInBlock = 0;
218  unsigned int columnsInBlock = 0;
219  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
220  {
221  curResultRow += rowsInBlock;
222  curWrapperRow += columnsInBlock;
223  if ( blockRow == 0)
224  {
225  rowsInBlock = rowSizes[blockRow] + 1;
226  columnsInBlock = colSizes[blockRow] + 1;
227  }
228  else
229  {
230  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
231  columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
232  }
233 
234  if( rowsInBlock == 0)
235  {
236  continue;
237  }
238  if( columnsInBlock == 0)
239  {
240  std::fill(result.begin()+curResultRow,
241  result.begin()+curResultRow + rowsInBlock, 0.0);
242  continue;
243  }
244 
245  const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockRow);
246  if( !block )
247  {
248  continue;
249  }
250 
251  double* resultWrapper = result_ptr + curResultRow;
252  const double* rhsWrapper = rhs_ptr + curWrapperRow;
253  Multiply(resultWrapper, *block, rhsWrapper);
254  //resultWrapper = (*block)*rhsWrapper;
255  }
256  curResultRow += rowsInBlock;
257  if (curResultRow < result.GetRows())
258  {
259  std::fill(result.begin()+curResultRow, result.end(), 0.0);
260  }
261  }
262 
264  const NekMatrix<NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>, BlockMatrixTag>& lhs,
265  const NekVector<double>& rhs)
266  {
267  unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
268  double* result_ptr = result.GetRawPtr();
269  const double* rhs_ptr = rhs.GetRawPtr();
270 
271  Array<OneD, unsigned int> rowSizes;
272  Array<OneD, unsigned int> colSizes;
273  lhs.GetBlockSizes(rowSizes, colSizes);
274 
275  unsigned int curResultRow = 0;
276  unsigned int curWrapperRow = 0;
277  unsigned int rowsInBlock = 0;
278  unsigned int columnsInBlock = 0;
279  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
280  {
281  curResultRow += rowsInBlock;
282  curWrapperRow += columnsInBlock;
283  if ( blockRow == 0)
284  {
285  rowsInBlock = rowSizes[blockRow] + 1;
286  columnsInBlock = colSizes[blockRow] + 1;
287  }
288  else
289  {
290  rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
291  columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
292  }
293 
294  if( rowsInBlock == 0)
295  {
296  continue;
297  }
298  if( columnsInBlock == 0)
299  {
300  std::fill(result.begin()+curResultRow,
301  result.begin()+curResultRow + rowsInBlock, 0.0);
302  continue;
303  }
304 
305  const DNekScalMat* block = lhs.GetBlockPtr(blockRow, blockRow);
306  if( !block )
307  {
308  continue;
309  }
310 
311  double* resultWrapper = result_ptr + curResultRow;
312  const double* rhsWrapper = rhs_ptr + curWrapperRow;
313 
314  // Multiply
315  const unsigned int* size = block->GetSize();
316  Blas::Dgemv('N', size[0], size[1], block->Scale(),
317  block->GetRawPtr(), size[0], rhsWrapper, 1,
318  0.0, resultWrapper, 1);
319  }
320  curResultRow += rowsInBlock;
321  if (curResultRow < result.GetRows())
322  {
323  std::fill(result.begin()+curResultRow, result.end(), 0.0);
324  }
325  }
326 
329  const NekDouble* rhs)
330  {
331  int vectorSize = lhs.GetColumns();
332  std::copy(rhs, rhs+vectorSize, result);
333  int n = lhs.GetRows();
334  const double* a = lhs.GetRawPtr();
335  double* x = result;
336  int incx = 1;
337 
338  Blas::Dtpmv('L', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
339  }
340 
342  const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>& lhs,
343  const NekDouble* rhs)
344  {
345  NekMultiplyLowerTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
346 
347  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
348  {
349  result[i] *= lhs.Scale();
350  }
351  }
352 
353  template<typename DataType, typename LhsDataType, typename MatrixType>
354  void NekMultiplyLowerTriangularMatrix(DataType* result,
356  const DataType* rhs)
357  {
358  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
359  {
360  DataType accum = DataType(0);
361  for(unsigned int j = 0; j <= i; ++j)
362  {
363  accum += lhs(i,j)*rhs[j];
364  }
365  result[i] = accum;
366  }
367  }
368 
371  const NekDouble* rhs)
372  {
373  int vectorSize = lhs.GetColumns();
374  std::copy(rhs, rhs+vectorSize, result);
375  int n = lhs.GetRows();
376  const double* a = lhs.GetRawPtr();
377  double* x = result;
378  int incx = 1;
379 
380  Blas::Dtpmv('U', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
381  }
382 
384  const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>& lhs,
385  const NekDouble* rhs)
386  {
387  NekMultiplyUpperTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
388 
389  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
390  {
391  result[i] *= lhs.Scale();
392  }
393  }
394 
395  template<typename DataType, typename LhsDataType, typename MatrixType>
396  void NekMultiplyUpperTriangularMatrix(DataType* result,
398  const DataType* rhs)
399  {
400  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
401  {
402  DataType accum = DataType(0);
403  for(unsigned int j = i; j < lhs.GetColumns(); ++j)
404  {
405  accum += lhs(i,j)*rhs[j];
406  }
407  result[i] = accum;
408  }
409  }
410 
411  template<typename InnerMatrixType, typename MatrixTag>
413  typename boost::enable_if
414  <
416  >::type* p=0)
417  {
418  const unsigned int* size = lhs.GetSize();
419 
420  double alpha = lhs.Scale();
421  const double* a = lhs.GetRawPtr();
422  const double* x = rhs;
423  int incx = 1;
424  double beta = 0.0;
425  double* y = result;
426  int incy = 1;
427 
428  Blas::Dspmv('U', size[0], alpha, a, x, incx, beta, y, incy);
429  }
430 
431  template<typename InnerMatrixType, typename MatrixTag>
433  typename boost::enable_if
434  <
436  >::type* p = 0)
437  {
438  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
439  }
440 
441  template<typename InnerMatrixType, typename MatrixTag>
443  typename boost::enable_if
444  <
446  >::type* p=0)
447  {
448  const unsigned int* size = lhs.GetSize();
449 
450  char t = lhs.GetTransposeFlag();
451 
452  double alpha = lhs.Scale();
453  const double* a = lhs.GetRawPtr();
454  int lda = size[0];
455  const double* x = rhs;
456  int incx = 1;
457  double beta = 0.0;
458  double* y = result;
459  int incy = 1;
460 
461  Blas::Dgemv(t, size[0], size[1], alpha, a, lda, x, incx, beta, y, incy);
462  }
463 
464  template<typename InnerMatrixType, typename MatrixTag>
466  typename boost::enable_if
467  <
469  >::type* p = 0)
470  {
471  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
472  }
473 
474  template<typename DataType, typename LhsDataType, typename MatrixType>
475  void Multiply(DataType* result,
477  const DataType* rhs)
478  {
479  switch(lhs.GetType())
480  {
481  case eFULL:
482  NekMultiplyFullMatrix(result, lhs, rhs);
483  break;
484  case eDIAGONAL:
485  NekMultiplyDiagonalMatrix(result, lhs, rhs);
486  break;
487  case eUPPER_TRIANGULAR:
488  NekMultiplyUpperTriangularMatrix(result, lhs, rhs);
489  break;
490  case eLOWER_TRIANGULAR:
491  NekMultiplyLowerTriangularMatrix(result, lhs, rhs);
492  break;
493  case eSYMMETRIC:
494  NekMultiplySymmetricMatrix(result, lhs, rhs);
495  break;
496  case eBANDED:
497  NekMultiplyBandedMatrix(result, lhs, rhs);
498  break;
499  case eSYMMETRIC_BANDED:
502  default:
503  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
504  }
505  }
506 
507  template<typename DataType, typename LhsDataType, typename MatrixType>
510  const NekVector<DataType>& rhs)
511  {
512 
513  ASSERTL1(lhs.GetColumns() == rhs.GetRows(), std::string("A left side matrix with column count ") +
514  boost::lexical_cast<std::string>(lhs.GetColumns()) +
515  std::string(" and a right side vector with row count ") +
516  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be multiplied."));
517  Multiply(result.GetRawPtr(), lhs, rhs.GetRawPtr());
518  }
519 
520  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)));
521 
522  template<typename DataType, typename LhsInnerMatrixType>
525  const NekVector<DataType>& rhs)
526  {
527  if( lhs.GetStorageType() == eDIAGONAL )
528  {
529  DiagonalBlockMatrixMultiply(result, lhs, rhs);
530  }
531  else
532  {
533  FullBlockMatrixMultiply(result, lhs, rhs);
534  }
535  }
536 
537  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_BLOCK_MATRIX_TYPES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)));
538 
539  template<typename DataType, typename LhsDataType, typename MatrixType>
540  NekVector<DataType>
542  const NekVector<DataType>& rhs)
543  {
544  NekVector<DataType> result(lhs.GetRows(), DataType(0));
545  Multiply(result, lhs, rhs);
546  return result;
547  }
548 
549  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (NekVector<NekDouble>)), (0, ()), (1,(const NekVector<NekDouble>&)));
550 
551 
552 }
void NekMultiplyUpperTriangularMatrix(NekDouble *result, const NekMatrix< NekDouble, StandardMatrixTag > &lhs, const NekDouble *rhs)
void FullBlockMatrixMultiply(NekVector< DataType > &result, const NekMatrix< LhsInnerMatrixType, BlockMatrixTag > &lhs, const NekVector< DataType > &rhs)
iterator begin()
Definition: NekVector.cpp:242
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
void DiagonalBlockFullScalMatrixMultiply(NekVector< double > &result, const NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > &lhs, const NekVector< double > &rhs)
#define NEKTAR_STANDARD_AND_SCALED_MATRICES
void NekMultiplyUnspecializedMatrixType(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
void NekMultiplyFullMatrix(NekDouble *result, const NekMatrix< InnerMatrixType, MatrixTag > &lhs, const NekDouble *rhs, typename boost::enable_if< CanGetRawPtr< NekMatrix< InnerMatrixType, MatrixTag > > >::type *p=0)
StandardMatrixTag & lhs
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
DataType * GetRawPtr()
Definition: NekVector.cpp:224
void NekMultiplySymmetricMatrix(NekDouble *result, const NekMatrix< InnerMatrixType, MatrixTag > &lhs, const NekDouble *rhs, typename boost::enable_if< CanGetRawPtr< NekMatrix< InnerMatrixType, MatrixTag > > >::type *p=0)
double NekDouble
void NekMultiplyLowerTriangularMatrix(NekDouble *result, const NekMatrix< NekDouble, StandardMatrixTag > &lhs, const NekDouble *rhs)
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:212
unsigned int GetRows() const
Definition: NekVector.cpp:218
void NekMultiplyBandedMatrix(NekDouble *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const NekDouble *rhs, typename boost::enable_if< CanGetRawPtr< NekMatrix< LhsDataType, MatrixType > > >::type *p=0)
iterator end()
Definition: NekVector.cpp:245
NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES,(1,(void)),(1,(DNekMat &)),(1,(const NekDouble &))) template< typename DataType
void NekMultiplyDiagonalMatrix(DataType *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const DataType *rhs)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void DiagonalBlockMatrixMultiply(NekVector< double > &result, const NekMatrix< LhsInnerMatrixType, BlockMatrixTag > &lhs, const NekVector< double > &rhs)
#define NEKTAR_BLOCK_MATRIX_TYPES
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:183