36 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STANDARD_MATRIX_HPP
37 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STANDARD_MATRIX_HPP
44 #include <ExpressionTemplates/ExpressionTemplates.hpp>
52 #include <ExpressionTemplates/Node.hpp>
53 #include <ExpressionTemplates/ExpressionEvaluator.hpp>
62 template<
typename DataType>
74 template<
typename T,
typename MatrixType>
83 typedef std::forward_iterator_tag
type;
87 struct TagType<const Z>
89 typedef std::input_iterator_tag
type;
94 typedef typename boost::call_traits<value_type>::reference
reference;
96 typedef typename boost::remove_reference<value_type>::type*
pointer;
102 m_curRow(
std::numeric_limits<unsigned int>::max()),
103 m_curColumn(
std::numeric_limits<unsigned int>::max()),
105 m_curIndex(
std::numeric_limits<unsigned int>::max()),
121 m_transpose(transpose)
125 m_curRow = std::numeric_limits<unsigned int>::max();
126 m_curColumn = std::numeric_limits<unsigned int>::max();
127 m_curIndex = std::numeric_limits<unsigned int>::max();
134 m_curRow(rhs.m_curRow),
135 m_curColumn(rhs.m_curColumn),
136 m_matrix(rhs.m_matrix),
137 m_curIndex(rhs.m_curIndex),
138 m_transpose(rhs.m_transpose)
142 iterator_impl<T, MatrixType>&
operator=(
const iterator_impl<T, MatrixType>& rhs)
146 m_curRow = rhs.m_curRow;
147 m_curColumn = rhs.m_curColumn;
148 m_matrix = rhs.m_matrix;
149 m_curIndex = rhs.m_curIndex;
150 m_transpose = rhs.m_transpose;
158 ASSERTL1(m_data < m_end,
"Attempt to dereference matrix iterator after its end.");
163 return m_matrix->GetPtr()[m_curIndex];
171 ASSERTL1(m_data < m_end,
"Attempt to dereference matrix iterator after its end.");
176 return m_matrix->GetPtr()[m_curIndex];
189 boost::tie(m_curRow, m_curColumn) =
190 m_matrix->Advance(m_curRow, m_curColumn, m_transpose);
191 if( m_curRow == std::numeric_limits<unsigned int>::max() )
193 m_curIndex = m_curRow;
197 m_curIndex = m_matrix->CalculateIndex(m_curRow, m_curColumn, m_transpose);
206 iterator_impl<T, MatrixType> result = *
this;
213 return m_data == rhs.m_data &&
214 m_end == rhs.m_end &&
215 m_curRow == rhs.m_curRow &&
216 m_curColumn == rhs.m_curColumn &&
217 m_matrix == rhs.m_matrix &&
218 m_curIndex == rhs.m_curIndex &&
219 m_transpose == rhs.m_transpose;
224 return !(*
this == rhs);
248 unsigned int subDiagonals = std::numeric_limits<unsigned int>::max(),
249 unsigned int superDiagonals = std::numeric_limits<unsigned int>::max());
253 unsigned int subDiagonals,
254 unsigned int superDiagonals,
255 unsigned int capacity);
266 unsigned int subDiagonals = std::numeric_limits<unsigned int>::max(),
267 unsigned int superDiagonals = std::numeric_limits<unsigned int>::max());
278 unsigned int subDiagonals = std::numeric_limits<unsigned int>::max(),
279 unsigned int superDiagonals = std::numeric_limits<unsigned int>::max());
290 unsigned int subDiagonals = std::numeric_limits<unsigned int>::max(),
291 unsigned int superDiagonals = std::numeric_limits<unsigned int>::max());
304 unsigned int subDiagonals = std::numeric_limits<unsigned int>::max(),
305 unsigned int superDiagonals = std::numeric_limits<unsigned int>::max());
309 #ifdef NEKTAR_USE_EXPRESSION_TEMPLATES
310 template<
typename L,
typename Op,
typename R>
311 NekMatrix(
const expt::Node<L, Op, R>& rhs) :
313 m_wrapperType(
eCopy),
314 m_storagePolicy(
eFULL),
315 m_numberOfSuperDiagonals(
std::numeric_limits<unsigned int>::max()),
316 m_numberOfSubDiagonals(
std::numeric_limits<unsigned int>::max()),
319 boost::tuple<unsigned int, unsigned int, unsigned int> sizes =
320 MatrixSize<expt::Node<L, Op, R>,
typename expt::Node<L, Op, R>::Indices, 0>::GetRequiredSize(rhs.GetData());
322 unsigned int rows = sizes.get<0>();
323 unsigned int columns = sizes.get<1>();
324 unsigned int bufferCapacity = sizes.get<2>();
327 this->Resize(rows, columns);
329 expt::ExpressionEvaluator::EvaluateWithoutAliasingCheck(rhs, *
this);
330 this->RemoveExcessCapacity();
333 #endif //NEKTAR_USE_EXPRESSION_TEMPLATES
340 template<
typename InnerMatrixType>
343 BaseType::operator=(rhs);
344 m_storagePolicy = rhs.GetType();
345 m_numberOfSubDiagonals = rhs.GetNumberOfSubDiagonals();
346 m_numberOfSuperDiagonals = rhs.GetNumberOfSuperDiagonals();
348 ResizeDataArrayIfNeeded();
350 unsigned int requiredStorageSize = GetRequiredStorageSize();
351 DataType scale = rhs.Scale();
353 DataType* lhs_array = m_data.data();
354 const DataType* rhs_array = rhs.GetRawPtr();
356 for(
unsigned int i = 0; i < requiredStorageSize; ++i)
358 lhs_array[i] = scale*rhs_array[i];
366 #ifdef NEKTAR_USE_EXPRESSION_TEMPLATES
367 template<
typename L,
typename Op,
typename R>
368 ThisType& operator=(
const expt::Node<L, Op, R>& rhs)
370 if( this->GetWrapperType() ==
eWrapper ||
371 expt::ExpressionEvaluator::ContainsAlias(rhs, *
this) )
388 boost::tuple<unsigned int, unsigned int, unsigned int> sizes =
389 MatrixSize<expt::Node<L, Op, R>,
typename expt::Node<L, Op, R>::Indices, 0>::GetRequiredSize(rhs.GetData());
390 unsigned int rows = sizes.get<0>();
391 unsigned int columns = sizes.get<1>();
392 unsigned int bufferSize = sizes.get<2>();
394 if( this->GetRows() != rows ||
395 this->GetColumns() != columns )
397 this->Resize(rows, columns);
402 this->ResizeDataArrayIfNeeded(bufferSize);
404 this->SetTransposeFlag(
'N');
406 expt::ExpressionEvaluator::Evaluate(rhs, *
this);
407 this->RemoveExcessCapacity();
412 #endif //NEKTAR_USE_EXPRESSION_TEMPLATES
420 Proxy() : m_value(defaultReturnValue) {}
421 explicit Proxy(DataType& value) : m_value(value) {}
422 Proxy(
const Proxy& rhs) : m_value(rhs.m_value) {}
425 m_value = rhs.m_value;
431 operator DataType&() {
return m_value; }
432 operator const DataType&()
const {
return m_value; }
435 if( &m_value != &defaultReturnValue )
451 LIB_UTILITIES_EXPORT ConstGetValueType operator()(
unsigned int row,
unsigned int column,
char transpose)
const;
455 LIB_UTILITIES_EXPORT unsigned int CalculateIndex(
unsigned int row,
unsigned int col,
const char transpose)
const;
458 LIB_UTILITIES_EXPORT typename boost::call_traits<DataType>::const_reference GetValue(
unsigned int row,
unsigned int column)
const;
465 LIB_UTILITIES_EXPORT ConstGetValueType GetValue(
unsigned int row,
unsigned int column,
char transpose)
const;
477 typedef iterator_impl<DataType, ThisType>
iterator;
519 Advance(
unsigned int curRow,
unsigned int curColumn)
const;
522 Advance(
unsigned int curRow,
unsigned int curColumn,
char transpose)
const;
526 LIB_UTILITIES_EXPORT static boost::shared_ptr<ThisType> CreateWrapper(
const boost::shared_ptr<ThisType>& rhs);
535 LIB_UTILITIES_EXPORT void SetValue(
unsigned int row,
unsigned int column,
typename boost::call_traits<DataType>::const_reference d);
537 LIB_UTILITIES_EXPORT void SetValue(
unsigned int row,
unsigned int column,
typename boost::call_traits<DataType>::const_reference d,
char transpose);
559 #ifdef NEKTAR_USE_EXPRESSION_TEMPLATES
560 expt::Node<expt::Node<NekMatrix<DataType, StandardMatrixTag> >, expt::NegateOp,
void >
operator-()
const
562 expt::Node<NekMatrix<DataType, StandardMatrixTag> > leafNode(*
this);
563 return expt::Node<expt::Node<NekMatrix<DataType, StandardMatrixTag> >, expt::NegateOp,
void >(leafNode);
584 LIB_UTILITIES_EXPORT virtual typename boost::call_traits<DataType>::value_type v_GetValue(
unsigned int row,
unsigned int column)
const;
591 LIB_UTILITIES_EXPORT virtual void v_SetValue(
unsigned int row,
unsigned int column,
typename boost::call_traits<DataType>::const_reference d);
606 template<
typename DataType>
609 template<
typename DataType>
613 template<
typename DataType>
617 #ifdef NEKTAR_USE_EXPRESSION_TEMPLATES
620 template<
typename DataType>
621 struct IsAlias<
Nektar::NekMatrix<DataType, Nektar::StandardMatrixTag>,
Nektar::NekMatrix<DataType, Nektar::StandardMatrixTag> >
625 return lhs.GetPtr().Overlaps(rhs.GetPtr());
629 template<
typename DataType,
typename NodeType,
typename Indices,
unsigned int StartIndex>
630 struct CreateFromTree<
Nektar::NekMatrix<DataType, Nektar::StandardMatrixTag>, NodeType, Indices, StartIndex>
632 template<
typename ArgVectorType>
635 boost::tuple<unsigned int, unsigned int, unsigned int> sizes =
636 Nektar::MatrixSize<NodeType, Indices, StartIndex>::GetRequiredSize(tree);
638 unsigned int rows = sizes.get<0>();
639 unsigned int columns = sizes.get<1>();
640 unsigned int bufferSize = sizes.get<2>();
644 std::numeric_limits<unsigned int>::max(), bufferSize);
651 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STANDARD_MATRIX_HPP
iterator_impl(pointer d, pointer e, bool isEnd=false)
iterator_impl< const DataType, const ThisType > const_iterator
unsigned int difference_type
static Array< OneD, NekDouble > NullNekDouble1DArray
boost::remove_reference< value_type >::type * pointer
TagType< T >::type iterator_category
MatrixStorage m_storagePolicy
bool operator!=(const iterator_impl< T, MatrixType > &rhs)
unsigned int m_numberOfSubDiagonals
Matrix< DataType > BaseType
bool operator==(const Array< OneD, NekDouble > &lhs, const Array< OneD, NekDouble > &rhs)
std::input_iterator_tag type
PointerWrapper m_wrapperType
Array< OneD, DataType > m_data
boost::call_traits< value_type >::const_reference const_reference
static DataType defaultReturnValue
boost::call_traits< value_type >::reference reference
NekMatrix< DataType, StandardMatrixTag > ThisType
ThisType & operator=(const NekMatrix< InnerMatrixType, ScaledMatrixTag > &rhs)
iterator_impl< T, MatrixType > operator++(int)
increment operator.
#define LIB_UTILITIES_EXPORT
NekMatrix< NekDouble > Invert(const NekMatrix< NekDouble > &matA)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
const DataType & operator*() const
void operator=(const DataType &newValue)
unsigned int m_numberOfSuperDiagonals
const_reference operator*() const
std::forward_iterator_tag type
NekPoint< DataType > operator-(const NekPoint< DataType > &lhs, const NekPoint< DataType > &rhs)
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
const DataType & ConstGetValueType
iterator_impl(MatrixType *m, char transpose, bool isEnd=false)
Array< OneD, DataType > m_tempSpace
iterator_impl< DataType, ThisType > iterator
iterator_impl< T, MatrixType > & operator++()
Prefix increment operator.
1D Array of constant elements with garbage collection and bounds checking.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void NegateInPlace(NekVector< DataType > &v)
iterator_impl(const iterator_impl< T, MatrixType > &rhs)
iterator_impl< T, MatrixType > & operator=(const iterator_impl< T, MatrixType > &rhs)
Proxy & operator=(const Proxy &rhs)
bool operator==(const iterator_impl< T, MatrixType > &rhs)