39 template<
typename DataType>
44 m_numberOfSuperDiagonals(
std::numeric_limits<unsigned int>::max()),
45 m_numberOfSubDiagonals(
std::numeric_limits<unsigned int>::max()),
51 template<
typename DataType>
53 unsigned int subDiagonals,
54 unsigned int superDiagonals) :
55 Matrix<DataType>(rows, columns, policy),
58 m_numberOfSuperDiagonals(superDiagonals),
59 m_numberOfSubDiagonals(subDiagonals),
65 template<
typename DataType>
68 unsigned int subDiagonals,
69 unsigned int superDiagonals,
70 unsigned int capacity) :
71 Matrix<DataType>(rows, columns, policy),
74 m_numberOfSuperDiagonals(superDiagonals),
75 m_numberOfSubDiagonals(subDiagonals),
78 unsigned int requiredStorage = this->GetRequiredStorageSize();
79 unsigned int actualSize = std::max(requiredStorage, capacity);
83 template<
typename DataType>
86 unsigned int subDiagonals,
87 unsigned int superDiagonals) :
88 Matrix<DataType>(rows, columns, policy),
91 m_numberOfSuperDiagonals(superDiagonals),
92 m_numberOfSubDiagonals(subDiagonals),
96 std::fill(m_data.begin(), m_data.end(), initValue);
99 template<
typename DataType>
102 unsigned int subDiagonals,
103 unsigned int superDiagonals) :
104 Matrix<DataType>(rows, columns, policy),
106 m_wrapperType(
eCopy),
107 m_numberOfSuperDiagonals(superDiagonals),
108 m_numberOfSubDiagonals(subDiagonals),
111 unsigned int size = GetRequiredStorageSize();
113 std::copy(data, data + size, m_data.begin());
116 template<
typename DataType>
119 unsigned int subDiagonals,
120 unsigned int superDiagonals) :
121 Matrix<DataType>(rows, columns, policy),
123 m_wrapperType(
eCopy),
124 m_numberOfSuperDiagonals(superDiagonals),
125 m_numberOfSubDiagonals(subDiagonals),
132 template<
typename DataType>
135 unsigned int subDiagonals,
136 unsigned int superDiagonals) :
137 Matrix<DataType>(rows, columns, policy),
139 m_wrapperType(wrapperType),
140 m_numberOfSuperDiagonals(superDiagonals),
141 m_numberOfSubDiagonals(subDiagonals),
155 template<
typename DataType>
159 m_wrapperType(rhs.m_wrapperType),
160 m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
161 m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
164 PerformCopyConstruction(rhs);
167 template<
typename DataType>
179 ResizeDataArrayIfNeeded();
181 unsigned int requiredStorageSize = GetRequiredStorageSize();
188 template<
typename DataType>
191 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
192 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
193 std::string(
" rows"));
194 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
195 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
196 std::string(
" columns"));
198 return this->GetValue(row, column, this->GetTransposeFlag());
201 template<
typename DataType>
204 return this->GetValue(row, column, transpose);
207 template<
typename DataType>
211 this->GetRows(), this->GetColumns(),
212 this->GetNumberOfSubDiagonals(), this->GetNumberOfSuperDiagonals());
215 template<
typename DataType>
218 unsigned int numRows = this->
GetSize()[0];
219 unsigned int numColumns = this->
GetSize()[1];
221 row, col, numRows, numColumns, transpose,
222 m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
225 template<
typename DataType>
228 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
229 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
230 std::string(
" rows"));
231 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
232 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
233 std::string(
" columns"));
235 return GetValue(row, column, this->GetTransposeFlag());
238 template<
typename DataType>
241 static DataType defaultReturnValue;
242 unsigned int index = CalculateIndex(row, column, transpose);
243 if( index != std::numeric_limits<unsigned int>::max() )
245 return m_data[index];
249 return defaultReturnValue;
253 template<
typename DataType>
259 template<
typename DataType>
265 template<
typename DataType>
268 return m_data.data();
271 template<
typename DataType>
274 return begin(this->GetTransposeFlag());
277 template<
typename DataType>
280 if( transpose ==
'N' )
282 return const_iterator(m_data.data(), m_data.data() + m_data.num_elements());
290 template<
typename DataType>
293 return end(this->GetTransposeFlag());
296 template<
typename DataType>
299 if( transpose ==
'N' )
301 return const_iterator(m_data.data(), m_data.data() + m_data.num_elements(),
true);
309 template<
typename DataType>
312 return m_data.num_elements();
315 template<
typename DataType>
318 if( m_numberOfSubDiagonals != std::numeric_limits<unsigned int>::max() )
320 return m_numberOfSubDiagonals;
322 else if( this->GetRows() > 0 )
324 return this->GetRows()-1;
332 template<
typename DataType>
335 if( m_numberOfSuperDiagonals != std::numeric_limits<unsigned int>::max() )
337 return m_numberOfSuperDiagonals;
339 else if( this->GetRows() > 0 )
341 return this->GetRows()-1;
349 template<
typename DataType>
352 return GetNumberOfSubDiagonals() + GetNumberOfSuperDiagonals() + 1;
355 template<
typename DataType>
358 if( this->GetRows() != rhs.
GetRows() ||
366 return std::equal(begin(), end(), rhs.
begin());
370 for(
unsigned int i = 0; i < this->GetRows(); ++i)
372 for(
unsigned int j = 0; j < this->GetColumns(); ++j)
374 if( (*
this)(i,j) != rhs(i,j) )
385 template<
typename DataType>
388 template<
typename DataType>
389 boost::tuples::tuple<unsigned int, unsigned int>
392 return Advance(curRow, curColumn, this->GetTransposeFlag());
395 template<
typename DataType>
396 boost::tuples::tuple<unsigned int, unsigned int>
399 unsigned int numRows = this->GetTransposedRows(transpose);
400 unsigned int numColumns = this->GetTransposedColumns(transpose);
402 switch(this->GetStorageType())
406 numRows, numColumns, curRow, curColumn);
410 numRows, numColumns, curRow, curColumn);
414 numRows, numColumns, curRow, curColumn);
419 numRows, numColumns, curRow, curColumn);
425 numRows, numColumns, curRow, curColumn);
429 numRows, numColumns, curRow, curColumn);
445 return boost::tuples::tuple<unsigned int, unsigned int>(curRow, curColumn);
448 template<
typename DataType>
454 template<
typename DataType>
460 template<
typename DataType>
464 m_wrapperType(wrapperType),
465 m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
466 m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals)
468 PerformCopyConstruction(rhs);
472 template<
typename DataType>
475 template<
typename DataType>
478 if( m_wrapperType ==
eCopy )
480 unsigned int requiredStorageSize = GetRequiredStorageSize();
481 if( m_data.num_elements() > requiredStorageSize )
484 CopyArrayN(m_data, newArray, requiredStorageSize);
488 else if( m_wrapperType ==
eWrapper )
490 ASSERTL0(
true,
"Can't call RemoveExcessCapacity on a wrapped matrix.");
494 template<
typename DataType>
497 if( m_wrapperType ==
eCopy )
499 if( m_data.num_elements() < requiredStorageSize )
502 std::copy(m_data.data(), m_data.data() + m_data.num_elements(), newData.
data());
506 else if( m_wrapperType ==
eWrapper )
510 ASSERTL0(m_data.num_elements() >= requiredStorageSize,
"Wrapped NekMatrices must have the same dimension in operator=");
514 template<
typename DataType>
517 unsigned int requiredStorageSize = GetRequiredStorageSize();
518 ResizeDataArrayIfNeeded(requiredStorageSize);
521 template<
typename DataType>
535 template<
typename DataType>
541 template<
typename DataType>
547 template<
typename DataType>
555 template<
typename DataType>
558 this->Resize(rows, cols);
565 this->GetData().ChangeSize(this->GetRequiredStorageSize());
566 ASSERTL0(this->GetRequiredStorageSize() <= this->GetData().num_elements(),
"Can't resize matrices if there is not enough capacity.");
570 template<
typename DataType>
573 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
574 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
575 std::string(
" rows"));
576 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
577 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
578 std::string(
" columns"));
580 return (*
this)(row, column, this->GetTransposeFlag());
583 template<
typename DataType>
586 unsigned int index = this->CalculateIndex(row, column, transpose);
587 if( index != std::numeric_limits<unsigned int>::max() )
589 return Proxy(this->GetData()[index]);
598 template<
typename DataType>
601 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
602 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
603 std::string(
" rows"));
604 ASSERTL2(column < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
605 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
606 std::string(
" columns"));
607 SetValue(row, column, d, this->GetTransposeFlag());
610 template<
typename DataType>
613 unsigned int index = this->CalculateIndex(row, column, transpose);
614 if( index != std::numeric_limits<unsigned int>::max() )
616 this->GetData()[index] = d;
624 template<
typename DataType>
627 return this->GetData();
630 template<
typename DataType>
633 return this->GetData().data();
636 template<
typename DataType>
639 return begin(this->GetTransposeFlag());
642 template<
typename DataType>
645 if( transpose ==
'N' )
647 return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements());
655 template<
typename DataType>
658 return end(this->GetTransposeFlag());
661 template<
typename DataType>
664 if( transpose ==
'N' )
666 return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements(),
true);
670 return iterator(
this, transpose,
true);
674 template<
typename DataType>
678 int nvals = this->GetColumns();
682 EigenSolve(EigValReal,EigValImag);
684 Vmath::Vmul(nvals,EigValReal,1,EigValReal,1,EigValReal,1);
685 Vmath::Vmul(nvals,EigValImag,1,EigValImag,1,EigValImag,1);
686 Vmath::Vadd(nvals,EigValReal,1,EigValImag,1,EigValReal,1);
693 template<
typename DataType>
698 ASSERTL0(this->GetRows()==this->GetColumns(),
"Only square matrices can be called");
700 switch(this->GetType())
704 this->GetData(), EigValReal,
705 EigValImag, EigVecs);
708 Vmath::Vcopy(this->GetRows(),&(this->GetData())[0],1, &EigValReal[0],1);
726 template<
typename DataType>
729 ASSERTL0(this->GetRows()==this->GetColumns(),
"Only square matrices can be inverted.");
730 ASSERTL0(this->GetTransposeFlag()==
'N',
"Only untransposed matrices may be inverted.");
732 switch(this->GetType())
736 this->GetData(), this->GetTransposeFlag());
766 template<
typename DataType>
769 if( m_tempSpace.capacity() == 0 )
776 template<
typename DataType>
779 std::swap(m_tempSpace, this->GetData());
782 #ifndef NEKTAR_USE_EXPRESSION_TEMPLATES
783 template<
typename DataType>
792 template<
typename DataType>
795 for(
unsigned int i = 0; i < this->GetPtr().num_elements(); ++i)
797 this->GetPtr()[i] *= s;
803 template<
typename DataType>
818 template<
typename DataType>
821 for(
unsigned int i = 0; i < m.
GetRows(); ++i)
823 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
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.
unsigned int m_numberOfSubDiagonals
char GetTransposeFlag() const
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)
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
#define LIB_UTILITIES_EXPORT
void CopyArrayN(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest, unsigned int n)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
int GetSize(const ConstArray< OneD, NekDouble > &x)
unsigned int m_numberOfSuperDiagonals
NekMatrix< NekDouble > Invert(const NekMatrix< NekDouble > &matA)
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
MatrixStorage GetType() 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 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)