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();
211 std::fill(result.
begin(), result.
end(), 0.0);
213 unsigned int curResultRow = 0;
214 unsigned int curWrapperRow = 0;
215 for(
unsigned int blockRow = 0; blockRow < numberOfBlockRows; ++blockRow)
220 curResultRow += lhs.GetNumberOfRowsInBlockRow(blockRow-1);
223 unsigned int blockColumn = blockRow;
224 if( blockColumn != 0 )
226 curWrapperRow += lhs.GetNumberOfColumnsInBlockColumn(blockColumn-1);
229 unsigned int rowsInBlock = lhs.GetNumberOfRowsInBlockRow(blockRow);
230 if( rowsInBlock == 0 )
235 unsigned int columnsInBlock = lhs.GetNumberOfColumnsInBlockColumn(blockColumn);
236 if( columnsInBlock == 0 )
241 const LhsInnerMatrixType* block = lhs.GetBlockPtr(blockRow, blockColumn);
247 double* resultWrapper = result_ptr + curResultRow;
248 const double* rhsWrapper = rhs_ptr + curWrapperRow;
249 Multiply(resultWrapper, *block, rhsWrapper);
258 int vectorSize = lhs.GetColumns();
259 std::copy(rhs, rhs+vectorSize, result);
260 int n = lhs.GetRows();
261 const double* a = lhs.GetRawPtr();
265 Blas::Dtpmv(
'L', lhs.GetTransposeFlag(),
'N', n, a, x, incx);
274 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
276 result[i] *=
lhs.Scale();
280 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
285 for(
unsigned int i = 0; i < lhs.GetRows(); ++i)
287 DataType accum = DataType(0);
288 for(
unsigned int j = 0; j <= i; ++j)
290 accum +=
lhs(i,j)*rhs[j];
300 int vectorSize = lhs.GetColumns();
301 std::copy(rhs, rhs+vectorSize, result);
302 int n = lhs.GetRows();
303 const double* a = lhs.GetRawPtr();
307 Blas::Dtpmv(
'U', lhs.GetTransposeFlag(),
'N', n, a, x, incx);
316 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
318 result[i] *=
lhs.Scale();
322 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
327 for(
unsigned int i = 0; i < lhs.GetRows(); ++i)
329 DataType accum = DataType(0);
330 for(
unsigned int j = i; j < lhs.GetColumns(); ++j)
332 accum +=
lhs(i,j)*rhs[j];
338 template<
typename InnerMatrixType,
typename MatrixTag>
340 typename boost::enable_if
345 int m = lhs.GetRows();
346 int n = lhs.GetColumns();
348 char t = lhs.GetTransposeFlag();
354 double alpha = lhs.Scale();
355 const double* a = lhs.GetRawPtr();
357 const double* x = rhs;
363 Blas::Dgemv(t, m, n, alpha, a, lda, x, incx, beta, y, incy);
366 template<
typename InnerMatrixType,
typename MatrixTag>
368 typename boost::enable_if
376 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
381 switch(lhs.GetType())
409 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
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."));
424 template<
typename DataType,
typename LhsInnerMatrixType>
441 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...
#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 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.