46 #include <ExpressionTemplates/CommutativeTraits.hpp>
47 #include <boost/type_traits.hpp>
62 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
67 for(
unsigned int i = 0; i < lhs.GetRows(); ++i)
69 DataType accum = DataType(0);
70 for(
unsigned int j = 0; j < lhs.GetColumns(); ++j)
72 accum +=
lhs(i,j)*rhs[j];
78 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
83 int n = lhs.GetRows();
84 for(
unsigned int i = 0; i < n; ++i)
86 result[i] =
lhs(i,i)*rhs[i];
90 template<
typename DataType,
typename LhsDataType>
95 int n = lhs.GetRows();
96 const DataType* mat_ptr = lhs.GetRawPtr();
100 template<
typename LhsDataType,
typename MatrixType>
104 typename boost::enable_if
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;
121 Blas::Dgbmv(
'N', m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy);
125 template<
typename DataType,
typename LhsDataType>
129 typename boost::enable_if
137 template<
typename DataType,
typename LhsInnerMatrixType>
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();
149 result_ptr[i] = DataType(0);
152 DataType* temp_ptr = temp.
get();
154 unsigned int curResultRow = 0;
155 for(
unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
157 unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
161 curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow-1);
164 if( rowsInBlock == 0 )
169 DataType* resultWrapper = result_ptr + curResultRow;
171 unsigned int curWrapperRow = 0;
172 for(
unsigned int blockColumn = 0; blockColumn < numberOfBlockColumns; ++blockColumn)
174 if( blockColumn != 0 )
176 curWrapperRow += lhs.GetNumberOfColumnsInBlockColumn(blockColumn-1);
180 const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockColumn);
186 unsigned int columnsInBlock = lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
187 if( columnsInBlock == 0 )
192 const DataType* rhsWrapper = rhs_ptr + curWrapperRow;
193 Multiply(temp_ptr, *block, rhsWrapper);
194 for(
unsigned int i = 0; i < rowsInBlock; ++i)
196 resultWrapper[i] += temp_ptr[i];
202 template<
typename LhsInnerMatrixType>
207 unsigned int numberOfBlockRows = lhs.GetNumberOfBlockRows();
213 lhs.GetBlockSizes(rowSizes, colSizes);
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)
221 curResultRow += rowsInBlock;
222 curWrapperRow += columnsInBlock;
225 rowsInBlock = rowSizes[blockRow] + 1;
226 columnsInBlock = colSizes[blockRow] + 1;
230 rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
231 columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
234 if( rowsInBlock == 0)
238 if( columnsInBlock == 0)
240 std::fill(result.
begin()+curResultRow,
241 result.
begin()+curResultRow + rowsInBlock, 0.0);
245 const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockRow);
251 double* resultWrapper = result_ptr + curResultRow;
252 const double* rhsWrapper = rhs_ptr + curWrapperRow;
253 Multiply(resultWrapper, *block, rhsWrapper);
256 curResultRow += rowsInBlock;
257 if (curResultRow < result.
GetRows())
259 std::fill(result.
begin()+curResultRow, result.
end(), 0.0);
267 unsigned int numberOfBlockRows =
lhs.GetNumberOfBlockRows();
273 lhs.GetBlockSizes(rowSizes, colSizes);
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)
281 curResultRow += rowsInBlock;
282 curWrapperRow += columnsInBlock;
285 rowsInBlock = rowSizes[blockRow] + 1;
286 columnsInBlock = colSizes[blockRow] + 1;
290 rowsInBlock = rowSizes[blockRow] - rowSizes[blockRow-1];
291 columnsInBlock = colSizes[blockRow] - colSizes[blockRow-1];
294 if( rowsInBlock == 0)
298 if( columnsInBlock == 0)
300 std::fill(result.
begin()+curResultRow,
301 result.
begin()+curResultRow + rowsInBlock, 0.0);
311 double* resultWrapper = result_ptr + curResultRow;
312 const double* rhsWrapper = rhs_ptr + curWrapperRow;
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);
320 curResultRow += rowsInBlock;
321 if (curResultRow < result.
GetRows())
323 std::fill(result.
begin()+curResultRow, result.
end(), 0.0);
331 int vectorSize = lhs.GetColumns();
333 int n = lhs.GetRows();
334 const double* a = lhs.GetRawPtr();
338 Blas::Dtpmv(
'L', lhs.GetTransposeFlag(),
'N', n, a, x, incx);
347 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
349 result[i] *=
lhs.Scale();
353 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
358 for(
unsigned int i = 0; i < lhs.GetRows(); ++i)
360 DataType accum = DataType(0);
361 for(
unsigned int j = 0; j <= i; ++j)
363 accum +=
lhs(i,j)*rhs[j];
373 int vectorSize = lhs.GetColumns();
375 int n = lhs.GetRows();
376 const double* a = lhs.GetRawPtr();
380 Blas::Dtpmv(
'U', lhs.GetTransposeFlag(),
'N', n, a, x, incx);
389 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
391 result[i] *=
lhs.Scale();
395 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
400 for(
unsigned int i = 0; i < lhs.GetRows(); ++i)
402 DataType accum = DataType(0);
403 for(
unsigned int j = i; j < lhs.GetColumns(); ++j)
405 accum +=
lhs(i,j)*rhs[j];
411 template<
typename InnerMatrixType,
typename MatrixTag>
413 typename boost::enable_if
418 const unsigned int* size = lhs.GetSize();
420 double alpha = lhs.Scale();
421 const double* a = lhs.GetRawPtr();
422 const double* x = rhs;
428 Blas::Dspmv(
'U', size[0], alpha, a, x, incx, beta, y, incy);
431 template<
typename InnerMatrixType,
typename MatrixTag>
433 typename boost::enable_if
441 template<
typename InnerMatrixType,
typename MatrixTag>
443 typename boost::enable_if
448 const unsigned int* size = lhs.GetSize();
450 char t = lhs.GetTransposeFlag();
452 double alpha = lhs.Scale();
453 const double* a = lhs.GetRawPtr();
455 const double* x = rhs;
461 Blas::Dgemv(t, size[0], size[1], alpha, a, lda, x, incx, beta, y, incy);
464 template<
typename InnerMatrixType,
typename MatrixTag>
466 typename boost::enable_if
474 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
479 switch(lhs.GetType())
507 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
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."));
522 template<
typename DataType,
typename LhsInnerMatrixType>
539 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
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)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
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)
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
void NekMultiplySymmetricMatrix(NekDouble *result, const NekMatrix< InnerMatrixType, MatrixTag > &lhs, const NekDouble *rhs, typename boost::enable_if< CanGetRawPtr< NekMatrix< InnerMatrixType, MatrixTag > > >::type *p=0)
void NekMultiplyLowerTriangularMatrix(NekDouble *result, const NekMatrix< NekDouble, StandardMatrixTag > &lhs, const NekDouble *rhs)
unsigned int GetDimension() const
Returns the number of dimensions for the point.
unsigned int GetRows() const
void NekMultiplyBandedMatrix(NekDouble *result, const NekMatrix< LhsDataType, MatrixType > &lhs, const NekDouble *rhs, typename boost::enable_if< CanGetRawPtr< NekMatrix< LhsDataType, MatrixType > > >::type *p=0)
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...
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.