Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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,
64  const NekMatrix<LhsDataType, MatrixType>& lhs,
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,
80  const NekMatrix<LhsDataType, MatrixType>& lhs,
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,
92  const NekMatrix<LhsDataType, StandardMatrixTag>& lhs,
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>
102  const NekMatrix<LhsDataType, MatrixType>& lhs,
103  const NekDouble* rhs,
104  typename boost::enable_if
105  <
106  CanGetRawPtr<NekMatrix<LhsDataType, MatrixType> >
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,
127  const NekMatrix<LhsDataType, BlockMatrixTag>& lhs,
128  const DataType* rhs,
129  typename boost::enable_if
130  <
131  boost::mpl::not_<CanGetRawPtr<NekMatrix<LhsDataType, BlockMatrixTag> > >
132  >::type* p=0)
133  {
134  NEKERROR(ErrorUtil::efatal, "Banded block matrix multiplication not yet implemented");
135  }
136 
137  template<typename DataType, typename LhsInnerMatrixType>
139  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
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>
204  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
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  std::fill(result.begin(), result.end(), 0.0);
212 
213  unsigned int curResultRow = 0;
214  unsigned int curWrapperRow = 0;
215  for(unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
216  {
217 
218  if( blockRow != 0 )
219  {
220  curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow-1);
221  }
222 
223  unsigned int blockColumn = blockRow;
224  if( blockColumn != 0 )
225  {
226  curWrapperRow += lhs.GetNumberOfColumnsInBlockColumn(blockColumn-1);
227  }
228 
229  unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
230  if( rowsInBlock == 0 )
231  {
232  continue;
233  }
234 
235  unsigned int columnsInBlock = lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
236  if( columnsInBlock == 0 )
237  {
238  continue;
239  }
240 
241  const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockColumn);
242  if( !block )
243  {
244  continue;
245  }
246 
247  double* resultWrapper = result_ptr + curResultRow;
248  const double* rhsWrapper = rhs_ptr + curWrapperRow;
249  Multiply(resultWrapper, *block, rhsWrapper);
250  //resultWrapper = (*block)*rhsWrapper;
251  }
252  }
253 
255  const NekMatrix<NekDouble, StandardMatrixTag>& lhs,
256  const NekDouble* rhs)
257  {
258  int vectorSize = lhs.GetColumns();
259  std::copy(rhs, rhs+vectorSize, result);
260  int n = lhs.GetRows();
261  const double* a = lhs.GetRawPtr();
262  double* x = result;
263  int incx = 1;
264 
265  Blas::Dtpmv('L', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
266  }
267 
269  const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>& lhs,
270  const NekDouble* rhs)
271  {
272  NekMultiplyLowerTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
273 
274  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
275  {
276  result[i] *= lhs.Scale();
277  }
278  }
279 
280  template<typename DataType, typename LhsDataType, typename MatrixType>
281  void NekMultiplyLowerTriangularMatrix(DataType* result,
282  const NekMatrix<LhsDataType, MatrixType>& lhs,
283  const DataType* rhs)
284  {
285  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
286  {
287  DataType accum = DataType(0);
288  for(unsigned int j = 0; j <= i; ++j)
289  {
290  accum += lhs(i,j)*rhs[j];
291  }
292  result[i] = accum;
293  }
294  }
295 
297  const NekMatrix<NekDouble, StandardMatrixTag>& lhs,
298  const NekDouble* rhs)
299  {
300  int vectorSize = lhs.GetColumns();
301  std::copy(rhs, rhs+vectorSize, result);
302  int n = lhs.GetRows();
303  const double* a = lhs.GetRawPtr();
304  double* x = result;
305  int incx = 1;
306 
307  Blas::Dtpmv('U', lhs.GetTransposeFlag(), 'N', n, a, x, incx);
308  }
309 
311  const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>& lhs,
312  const NekDouble* rhs)
313  {
314  NekMultiplyUpperTriangularMatrix(result, *lhs.GetOwnedMatrix(), rhs);
315 
316  for(unsigned int i = 0; i < lhs.GetColumns(); ++i)
317  {
318  result[i] *= lhs.Scale();
319  }
320  }
321 
322  template<typename DataType, typename LhsDataType, typename MatrixType>
323  void NekMultiplyUpperTriangularMatrix(DataType* result,
324  const NekMatrix<LhsDataType, MatrixType>& lhs,
325  const DataType* rhs)
326  {
327  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
328  {
329  DataType accum = DataType(0);
330  for(unsigned int j = i; j < lhs.GetColumns(); ++j)
331  {
332  accum += lhs(i,j)*rhs[j];
333  }
334  result[i] = accum;
335  }
336  }
337 
338  template<typename InnerMatrixType, typename MatrixTag>
339  void NekMultiplyFullMatrix(NekDouble* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const NekDouble* rhs,
340  typename boost::enable_if
341  <
342  CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> >
343  >::type* p=0)
344  {
345  int m = lhs.GetRows();
346  int n = lhs.GetColumns();
347 
348  char t = lhs.GetTransposeFlag();
349  if( t == 'T' )
350  {
351  std::swap(m, n);
352  }
353 
354  double alpha = lhs.Scale();
355  const double* a = lhs.GetRawPtr();
356  int lda = m;
357  const double* x = rhs;
358  int incx = 1;
359  double beta = 0.0;
360  double* y = result;
361  int incy = 1;
362 
363  Blas::Dgemv(t, m, n, alpha, a, lda, x, incx, beta, y, incy);
364  }
365 
366  template<typename InnerMatrixType, typename MatrixTag>
367  void NekMultiplyFullMatrix(NekDouble* result, const NekMatrix<InnerMatrixType, MatrixTag>& lhs, const NekDouble* rhs,
368  typename boost::enable_if
369  <
370  boost::mpl::not_<CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> > >
371  >::type* p = 0)
372  {
373  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
374  }
375 
376  template<typename DataType, typename LhsDataType, typename MatrixType>
377  void Multiply(DataType* result,
378  const NekMatrix<LhsDataType, MatrixType>& lhs,
379  const DataType* rhs)
380  {
381  switch(lhs.GetType())
382  {
383  case eFULL:
384  NekMultiplyFullMatrix(result, lhs, rhs);
385  break;
386  case eDIAGONAL:
387  NekMultiplyDiagonalMatrix(result, lhs, rhs);
388  break;
389  case eUPPER_TRIANGULAR:
390  NekMultiplyUpperTriangularMatrix(result, lhs, rhs);
391  break;
392  case eLOWER_TRIANGULAR:
393  NekMultiplyLowerTriangularMatrix(result, lhs, rhs);
394  break;
395  case eSYMMETRIC:
396  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
397  break;
398  case eBANDED:
399  NekMultiplyBandedMatrix(result, lhs, rhs);
400  break;
401  case eSYMMETRIC_BANDED:
404  default:
405  NekMultiplyUnspecializedMatrixType(result, lhs, rhs);
406  }
407  }
408 
409  template<typename DataType, typename LhsDataType, typename MatrixType>
411  const NekMatrix<LhsDataType, MatrixType>& lhs,
412  const NekVector<DataType>& rhs)
413  {
414 
415  ASSERTL1(lhs.GetColumns() == rhs.GetRows(), std::string("A left side matrix with column count ") +
416  boost::lexical_cast<std::string>(lhs.GetColumns()) +
417  std::string(" and a right side vector with row count ") +
418  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be multiplied."));
419  Multiply(result.GetRawPtr(), lhs, rhs.GetRawPtr());
420  }
421 
422  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_STANDARD_AND_SCALED_MATRICES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)));
423 
424  template<typename DataType, typename LhsInnerMatrixType>
426  const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>& lhs,
427  const NekVector<DataType>& rhs)
428  {
429  if( lhs.GetStorageType() == eDIAGONAL )
430  {
431  DiagonalBlockMatrixMultiply(result, lhs, rhs);
432  }
433  else
434  {
435  FullBlockMatrixMultiply(result, lhs, rhs);
436  }
437  }
438 
439  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_BLOCK_MATRIX_TYPES, (1, (void)), (1, (NekVector<NekDouble>&)), (1,(const NekVector<NekDouble>&)));
440 
441  template<typename DataType, typename LhsDataType, typename MatrixType>
442  NekVector<DataType>
443  Multiply(const NekMatrix<LhsDataType, MatrixType>& lhs,
444  const NekVector<DataType>& rhs)
445  {
446  NekVector<DataType> result(lhs.GetRows(), DataType(0));
447  Multiply(result, lhs, rhs);
448  return result;
449  }
450 
451  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (NekVector<NekDouble>)), (0, ()), (1,(const NekVector<NekDouble>&)));
452 
453 
454 }