46 #include <ExpressionTemplates/CommutativeTraits.hpp>
47 #include <boost/type_traits.hpp>
62 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
64 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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>
80 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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>
92 const NekMatrix<LhsDataType, StandardMatrixTag>&
lhs,
95 int n = lhs.GetRows();
96 const DataType* mat_ptr = lhs.GetRawPtr();
100 template<
typename LhsDataType,
typename MatrixType>
102 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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>
127 const NekMatrix<LhsDataType, BlockMatrixTag>&
lhs,
129 typename boost::enable_if
131 boost::mpl::not_<
CanGetRawPtr<NekMatrix<LhsDataType, BlockMatrixTag> > >
137 template<
typename DataType,
typename LhsInnerMatrixType>
139 const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>&
lhs,
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>
204 const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>&
lhs,
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);
255 const NekMatrix<NekDouble, StandardMatrixTag>&
lhs,
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);
269 const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>&
lhs,
274 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
276 result[i] *=
lhs.Scale();
280 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
282 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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];
297 const NekMatrix<NekDouble, StandardMatrixTag>&
lhs,
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);
311 const NekMatrix<NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>&
lhs,
316 for(
unsigned int i = 0; i <
lhs.GetColumns(); ++i)
318 result[i] *=
lhs.Scale();
322 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
324 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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
370 boost::mpl::not_<
CanGetRawPtr<NekMatrix<InnerMatrixType, MatrixTag> > >
376 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
378 const NekMatrix<LhsDataType, MatrixType>&
lhs,
381 switch(lhs.GetType())
409 template<
typename DataType,
typename LhsDataType,
typename MatrixType>
411 const NekMatrix<LhsDataType, MatrixType>&
lhs,
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>
426 const NekMatrix<LhsInnerMatrixType, BlockMatrixTag>&
lhs,
441 template<
typename DataType,
typename LhsDataType,
typename MatrixType>