39 template<
typename DataType>
44 m_storagePolicy(
eFULL),
45 m_numberOfSuperDiagonals(
std::numeric_limits<unsigned int>::max()),
46 m_numberOfSubDiagonals(
std::numeric_limits<unsigned int>::max()),
52 template<
typename DataType>
54 unsigned int subDiagonals,
55 unsigned int superDiagonals) :
56 Matrix<DataType>(rows, columns),
59 m_storagePolicy(policy),
60 m_numberOfSuperDiagonals(superDiagonals),
61 m_numberOfSubDiagonals(subDiagonals),
67 template<
typename DataType>
70 unsigned int subDiagonals,
71 unsigned int superDiagonals,
72 unsigned int capacity) :
73 Matrix<DataType>(rows, columns),
76 m_storagePolicy(policy),
77 m_numberOfSuperDiagonals(superDiagonals),
78 m_numberOfSubDiagonals(subDiagonals),
81 unsigned int requiredStorage = this->GetRequiredStorageSize();
82 unsigned int actualSize = std::max(requiredStorage, capacity);
86 template<
typename DataType>
89 unsigned int subDiagonals,
90 unsigned int superDiagonals) :
91 Matrix<DataType>(rows, columns),
94 m_storagePolicy(policy),
95 m_numberOfSuperDiagonals(superDiagonals),
96 m_numberOfSubDiagonals(subDiagonals),
100 std::fill(m_data.begin(), m_data.end(), initValue);
103 template<
typename DataType>
106 unsigned int subDiagonals,
107 unsigned int superDiagonals) :
108 Matrix<DataType>(rows, columns),
110 m_wrapperType(
eCopy),
111 m_storagePolicy(policy),
112 m_numberOfSuperDiagonals(superDiagonals),
113 m_numberOfSubDiagonals(subDiagonals),
116 unsigned int size = GetRequiredStorageSize();
118 std::copy(data, data + size, m_data.begin());
121 template<
typename DataType>
124 unsigned int subDiagonals,
125 unsigned int superDiagonals) :
126 Matrix<DataType>(rows, columns),
128 m_wrapperType(
eCopy),
129 m_storagePolicy(policy),
130 m_numberOfSuperDiagonals(superDiagonals),
131 m_numberOfSubDiagonals(subDiagonals),
138 template<
typename DataType>
141 unsigned int subDiagonals,
142 unsigned int superDiagonals) :
143 Matrix<DataType>(rows, columns),
145 m_wrapperType(wrapperType),
146 m_storagePolicy(policy),
147 m_numberOfSuperDiagonals(superDiagonals),
148 m_numberOfSubDiagonals(subDiagonals),
162 template<
typename DataType>
166 m_wrapperType(rhs.m_wrapperType),
167 m_storagePolicy(rhs.m_storagePolicy),
168 m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
169 m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
172 PerformCopyConstruction(rhs);
176 template<
typename DataType>
179 template<
typename DataType>
192 ResizeDataArrayIfNeeded();
194 unsigned int requiredStorageSize = GetRequiredStorageSize();
201 template<
typename DataType>
204 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
205 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
206 std::string(
" rows"));
207 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
208 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
209 std::string(
" columns"));
211 return this->GetValue(row, column, this->GetTransposeFlag());
214 template<
typename DataType>
217 return this->GetValue(row, column, transpose);
220 template<
typename DataType>
224 this->GetRows(), this->GetColumns(),
225 this->GetNumberOfSubDiagonals(), this->GetNumberOfSuperDiagonals());
228 template<
typename DataType>
231 unsigned int numRows = this->
GetSize()[0];
232 unsigned int numColumns = this->
GetSize()[1];
234 row, col, numRows, numColumns, transpose,
235 m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
238 template<
typename DataType>
241 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
242 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
243 std::string(
" rows"));
244 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
245 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
246 std::string(
" columns"));
248 return GetValue(row, column, this->GetTransposeFlag());
251 template<
typename DataType>
254 static DataType defaultReturnValue;
255 unsigned int index = CalculateIndex(row, column, transpose);
256 if( index != std::numeric_limits<unsigned int>::max() )
258 return m_data[index];
262 return defaultReturnValue;
266 template<
typename DataType>
272 template<
typename DataType>
278 template<
typename DataType>
281 return m_data.data();
284 template<
typename DataType>
287 return begin(this->GetTransposeFlag());
290 template<
typename DataType>
293 if( transpose ==
'N' )
295 return const_iterator(m_data.data(), m_data.data() + m_data.num_elements());
303 template<
typename DataType>
306 return end(this->GetTransposeFlag());
309 template<
typename DataType>
312 if( transpose ==
'N' )
314 return const_iterator(m_data.data(), m_data.data() + m_data.num_elements(),
true);
322 template<
typename DataType>
325 return m_data.num_elements();
328 template<
typename DataType>
331 if( m_numberOfSubDiagonals != std::numeric_limits<unsigned int>::max() )
333 return m_numberOfSubDiagonals;
335 else if( this->GetRows() > 0 )
337 return this->GetRows()-1;
345 template<
typename DataType>
348 if( m_numberOfSuperDiagonals != std::numeric_limits<unsigned int>::max() )
350 return m_numberOfSuperDiagonals;
352 else if( this->GetRows() > 0 )
354 return this->GetRows()-1;
362 template<
typename DataType>
365 return GetNumberOfSubDiagonals() + GetNumberOfSuperDiagonals() + 1;
368 template<
typename DataType>
371 if( this->GetRows() != rhs.
GetRows() ||
379 return std::equal(begin(), end(), rhs.
begin());
383 for(
unsigned int i = 0; i < this->GetRows(); ++i)
385 for(
unsigned int j = 0; j < this->GetColumns(); ++j)
387 if( (*
this)(i,j) != rhs(i,j) )
398 template<
typename DataType>
401 template<
typename DataType>
404 return this->GetRawTransposeFlag();
407 template<
typename DataType>
408 boost::tuples::tuple<unsigned int, unsigned int>
411 return Advance(curRow, curColumn, this->GetTransposeFlag());
414 template<
typename DataType>
415 boost::tuples::tuple<unsigned int, unsigned int>
418 unsigned int numRows = this->GetTransposedRows(transpose);
419 unsigned int numColumns = this->GetTransposedColumns(transpose);
421 switch(m_storagePolicy)
425 numRows, numColumns, curRow, curColumn);
429 numRows, numColumns, curRow, curColumn);
433 numRows, numColumns, curRow, curColumn);
438 numRows, numColumns, curRow, curColumn);
444 numRows, numColumns, curRow, curColumn);
448 numRows, numColumns, curRow, curColumn);
464 return boost::tuples::tuple<unsigned int, unsigned int>(curRow, curColumn);
467 template<
typename DataType>
473 template<
typename DataType>
479 template<
typename DataType>
483 m_wrapperType(wrapperType),
484 m_storagePolicy(rhs.m_storagePolicy),
485 m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
486 m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals)
488 PerformCopyConstruction(rhs);
492 template<
typename DataType>
495 template<
typename DataType>
498 if( m_wrapperType ==
eCopy )
500 unsigned int requiredStorageSize = GetRequiredStorageSize();
501 if( m_data.num_elements() > requiredStorageSize )
504 CopyArrayN(m_data, newArray, requiredStorageSize);
508 else if( m_wrapperType ==
eWrapper )
510 ASSERTL0(
true,
"Can't call RemoveExcessCapacity on a wrapped matrix.");
514 template<
typename DataType>
517 if( m_wrapperType ==
eCopy )
519 if( m_data.num_elements() < requiredStorageSize )
522 std::copy(m_data.data(), m_data.data() + m_data.num_elements(), newData.
data());
526 else if( m_wrapperType ==
eWrapper )
530 ASSERTL0(m_data.num_elements() >= requiredStorageSize,
"Wrapped NekMatrices must have the same dimension in operator=");
534 template<
typename DataType>
537 unsigned int requiredStorageSize = GetRequiredStorageSize();
538 ResizeDataArrayIfNeeded(requiredStorageSize);
541 template<
typename DataType>
555 template<
typename DataType>
561 template<
typename DataType>
567 template<
typename DataType>
573 template<
typename DataType>
581 template<
typename DataType>
584 this->Resize(rows, cols);
591 this->GetData().ChangeSize(this->GetRequiredStorageSize());
592 ASSERTL0(this->GetRequiredStorageSize() <= this->GetData().num_elements(),
"Can't resize matrices if there is not enough capacity.");
596 template<
typename DataType>
599 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
600 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
601 std::string(
" rows"));
602 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
603 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
604 std::string(
" columns"));
606 return (*
this)(row, column, this->GetTransposeFlag());
609 template<
typename DataType>
612 unsigned int index = this->CalculateIndex(row, column, transpose);
613 if( index != std::numeric_limits<unsigned int>::max() )
615 return Proxy(this->GetData()[index]);
624 template<
typename DataType>
627 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
628 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
629 std::string(
" rows"));
630 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
631 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
632 std::string(
" columns"));
633 SetValue(row, column, d, this->GetTransposeFlag());
636 template<
typename DataType>
639 unsigned int index = this->CalculateIndex(row, column, transpose);
640 if( index != std::numeric_limits<unsigned int>::max() )
642 this->GetData()[index] = d;
650 template<
typename DataType>
653 return this->GetData();
656 template<
typename DataType>
659 return this->GetData().data();
662 template<
typename DataType>
665 return begin(this->GetTransposeFlag());
668 template<
typename DataType>
671 if( transpose ==
'N' )
673 return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements());
681 template<
typename DataType>
684 return end(this->GetTransposeFlag());
687 template<
typename DataType>
690 if( transpose ==
'N' )
692 return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements(),
true);
696 return iterator(
this, transpose,
true);
700 template<
typename DataType>
704 int nvals = this->GetColumns();
708 EigenSolve(EigValReal,EigValImag);
710 Vmath::Vmul(nvals,EigValReal,1,EigValReal,1,EigValReal,1);
711 Vmath::Vmul(nvals,EigValImag,1,EigValImag,1,EigValImag,1);
712 Vmath::Vadd(nvals,EigValReal,1,EigValImag,1,EigValReal,1);
719 template<
typename DataType>
724 ASSERTL0(this->GetRows()==this->GetColumns(),
"Only square matrices can be called");
726 switch(this->GetType())
730 this->GetData(), EigValReal,
731 EigValImag, EigVecs);
734 Vmath::Vcopy(this->GetRows(),&(this->GetData())[0],1, &EigValReal[0],1);
752 template<
typename DataType>
755 ASSERTL0(this->GetRows()==this->GetColumns(),
"Only square matrices can be inverted.");
756 ASSERTL0(this->GetTransposeFlag()==
'N',
"Only untransposed matrices may be inverted.");
758 switch(this->GetType())
762 this->GetData(), this->GetTransposeFlag());
789 template<
typename DataType>
792 if( m_tempSpace.capacity() == 0 )
799 template<
typename DataType>
802 std::swap(m_tempSpace, this->GetData());
805 #ifndef NEKTAR_USE_EXPRESSION_TEMPLATES
806 template<
typename DataType>
815 template<
typename DataType>
818 for(
unsigned int i = 0; i < this->GetPtr().num_elements(); ++i)
820 this->GetPtr()[i] *= s;
826 template<
typename DataType>
841 template<
typename DataType>
844 for(
unsigned int i = 0; i < m.
GetRows(); ++i)
846 for(
unsigned int j = 0; j < m.
GetColumns(); ++j)
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
iterator_impl< const DataType, const ThisType > const_iterator
MatrixStorage GetType() const
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
MatrixStorage m_storagePolicy
unsigned int m_numberOfSubDiagonals
unsigned int GetNumberOfSubDiagonals() const
Array< OneD, DataType > m_data
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
unsigned int GetColumns() const
size_type num_elements() const
Returns the array's size.
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
char GetTransposeFlag() const
Matrix< DataType > & operator=(const Matrix< DataType > &rhs)
static unsigned int CalculateIndex(MatrixStorage type, unsigned int row, unsigned int col, unsigned int numRows, unsigned int numColumns, const char transpose= 'N', unsigned int numSubDiags=0, unsigned int numSuperDiags=0)
unsigned int GetNumberOfSuperDiagonals() const
int GetSize(const Array< OneD, const NekDouble > &x)
#define LIB_UTILITIES_EXPORT
void CopyArrayN(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest, unsigned int n)
NekMatrix< NekDouble > Invert(const NekMatrix< NekDouble > &matA)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
unsigned int m_numberOfSuperDiagonals
NekPoint< DataType > operator-(const NekPoint< DataType > &lhs, const NekPoint< DataType > &rhs)
const Array< OneD, const DataType > & GetPtr() const
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
unsigned int GetRows() const
static void EigenSolve(unsigned int n, const Array< OneD, const double > &A, Array< OneD, NekDouble > &EigValReal, Array< OneD, NekDouble > &EigValImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn, char transpose= 'N')
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
static unsigned int GetRequiredStorageSize(MatrixStorage type, unsigned int rows, unsigned int columns, unsigned int subDiags=0, unsigned int superDiags=0)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
void Zero(int n, T *x, const int incx)
Zero vector.
1D Array of constant elements with garbage collection and bounds checking.
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void NegateInPlace(NekVector< DataType > &v)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
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.
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)