41 template<
typename DataType,
typename InnerMatrixType>
42 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
MatrixStorage type) :
48 m_numberOfBlockRows(0),
49 m_numberOfBlockColumns(0),
54 template<
typename DataType,
typename InnerMatrixType>
55 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
unsigned int numberOfBlockRows,
unsigned int numberOfBlockColumns,
56 unsigned int rowsPerBlock,
unsigned int columnsPerBlock,
58 BaseType(numberOfBlockRows*rowsPerBlock, numberOfBlockColumns*columnsPerBlock),
60 m_rowSizes(numberOfBlockRows),
61 m_columnSizes(numberOfBlockColumns),
63 m_numberOfBlockRows(numberOfBlockRows),
64 m_numberOfBlockColumns(numberOfBlockColumns),
67 m_storageSize = GetRequiredStorageSize();
68 m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
69 for(
unsigned int i = 1; i <= numberOfBlockRows; ++i)
71 m_rowSizes[i-1] = i*rowsPerBlock-1;
74 for(
unsigned int i = 1; i <= numberOfBlockColumns; ++i)
76 m_columnSizes[i-1] = i*columnsPerBlock-1;
80 template<
typename DataType,
typename InnerMatrixType>
81 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
unsigned int numberOfBlockRows,
unsigned int numberOfBlockColumns,
82 const unsigned int* rowsPerBlock,
const unsigned int* columnsPerBlock,
84 BaseType(std::accumulate(rowsPerBlock, rowsPerBlock + numberOfBlockRows, 0),
85 std::accumulate(columnsPerBlock, columnsPerBlock + numberOfBlockColumns, 0)),
87 m_rowSizes(numberOfBlockRows),
88 m_columnSizes(numberOfBlockColumns),
90 m_numberOfBlockRows(numberOfBlockRows),
91 m_numberOfBlockColumns(numberOfBlockColumns),
94 m_storageSize = GetRequiredStorageSize();
95 m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
96 Initialize(rowsPerBlock, columnsPerBlock);
99 template<
typename DataType,
typename InnerMatrixType>
100 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
unsigned int numberOfBlockRows,
unsigned int numberOfBlockColumns,
101 const Array<OneD, const unsigned int>& rowsPerBlock,
const Array<OneD, const unsigned int>& columnsPerBlock,
103 BaseType(std::accumulate(rowsPerBlock.data(), rowsPerBlock.data() + numberOfBlockRows, 0),
104 std::accumulate(columnsPerBlock.data(), columnsPerBlock.data() + numberOfBlockColumns, 0)),
106 m_rowSizes(numberOfBlockRows),
107 m_columnSizes(numberOfBlockColumns),
109 m_numberOfBlockRows(numberOfBlockRows),
110 m_numberOfBlockColumns(numberOfBlockColumns),
113 m_storageSize = GetRequiredStorageSize();
114 m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
115 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
118 template<
typename DataType,
typename InnerMatrixType>
119 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
const Array<OneD, const unsigned int>& rowsPerBlock,
120 const Array<OneD, const unsigned int>& columnsPerBlock,
122 BaseType(std::accumulate(rowsPerBlock.begin(), rowsPerBlock.end(), 0),
123 std::accumulate(columnsPerBlock.begin(), columnsPerBlock.end(), 0)),
125 m_rowSizes(rowsPerBlock.num_elements()),
126 m_columnSizes(columnsPerBlock.num_elements()),
128 m_numberOfBlockRows(rowsPerBlock.num_elements()),
129 m_numberOfBlockColumns(columnsPerBlock.num_elements()),
132 m_storageSize = GetRequiredStorageSize();
133 m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
134 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
137 template<
typename DataType,
typename InnerMatrixType>
138 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
const NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>& rhs) :
140 m_data(rhs.m_data.num_elements()),
141 m_rowSizes(rhs.m_rowSizes),
142 m_columnSizes(rhs.m_columnSizes),
143 m_storageSize(rhs.m_storageSize),
144 m_numberOfBlockRows(rhs.m_numberOfBlockRows),
145 m_numberOfBlockColumns(rhs.m_numberOfBlockColumns),
146 m_storageType(rhs.m_storageType)
148 for(
unsigned int i = 0; i < rhs.m_data.num_elements(); ++i)
150 m_data[i] = InnerType::CreateWrapper(rhs.m_data[i]);
154 template<
typename DataType,
typename InnerMatrixType>
155 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetRequiredStorageSize()
const
157 return BaseType::GetRequiredStorageSize(this->GetStorageType(), this->GetNumberOfBlockRows(),
158 this->GetNumberOfBlockColumns());
161 template<
typename DataType,
typename InnerMatrixType>
162 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CalculateBlockIndex(
unsigned int row,
unsigned int column)
const
164 unsigned int blockRows = GetNumberOfBlockRows();
165 unsigned int blockCols = GetNumberOfBlockColumns();
167 if( this->GetTransposeFlag() ==
'T' )
169 std::swap(blockRows, blockCols);
172 return BaseType::CalculateIndex(this->GetStorageType(),
173 row, column, blockRows, blockCols, this->GetTransposeFlag());
176 template<
typename DataType,
typename InnerMatrixType>
177 const typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::InnerType*
178 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetBlockPtr(
unsigned int row,
unsigned int column)
const
180 ASSERTL2(row < m_numberOfBlockRows, std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
181 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockRows) +
182 std::string(
" rows"));
183 ASSERTL2(column < m_numberOfBlockColumns, std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
184 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockColumns) +
185 std::string(
" columns"));
186 int x = CalculateBlockIndex(row,column);
193 return m_data[x].get();
197 template<
typename DataType,
typename InnerMatrixType>
198 boost::shared_ptr<const typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::InnerType>
199 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetBlock(
unsigned int row,
unsigned int column)
const
201 ASSERTL2(row < m_numberOfBlockRows, std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
202 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockRows) +
203 std::string(
" rows"));
204 ASSERTL2(column < m_numberOfBlockColumns, std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
205 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockColumns) +
206 std::string(
" columns"));
207 int x = CalculateBlockIndex(row,column);
210 return boost::shared_ptr<const InnerType>();
218 template<
typename DataType,
typename InnerMatrixType>
219 boost::shared_ptr<typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::InnerType>&
220 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetBlock(
unsigned int row,
unsigned int column)
222 ASSERTL2(row < m_numberOfBlockRows, std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
223 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockRows) +
224 std::string(
" rows"));
225 ASSERTL2(column < m_numberOfBlockColumns, std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
226 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockColumns) +
227 std::string(
" columns"));
228 int x = CalculateBlockIndex(row,column);
231 return m_nullBlockPtr;
239 template<
typename DataType,
typename InnerMatrixType>
240 void NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::SetBlock(
unsigned int row,
unsigned int column, boost::shared_ptr<InnerType>& m)
242 ASSERTL2(row < m_numberOfBlockRows, std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
243 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockRows) +
244 std::string(
" rows"));
245 ASSERTL2(column < m_numberOfBlockColumns, std::string(
"Column ") + boost::lexical_cast<std::string>(column) +
246 std::string(
" requested in a block matrix with a maximum of ") + boost::lexical_cast<std::string>(m_numberOfBlockColumns) +
247 std::string(
" columns"));
248 m_data[CalculateBlockIndex(row, column)] = InnerType::CreateWrapper(m);
253 template<
typename DataType,
typename InnerMatrixType>
254 typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::ConstGetValueType
255 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::operator()(
unsigned int row,
unsigned int col)
const
257 ASSERTL2(row < this->GetRows(), std::string(
"Row ") + boost::lexical_cast<std::string>(row) +
258 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
259 std::string(
" rows"));
260 ASSERTL2(col < this->GetColumns(), std::string(
"Column ") + boost::lexical_cast<std::string>(col) +
261 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
262 std::string(
" columns"));
265 const Array<OneD, unsigned int>* rowSizes = &m_rowSizes;
266 const Array<OneD, unsigned int>* columnSizes = &m_columnSizes;
268 if( this->GetTransposeFlag() ==
'T' )
270 std::swap(rowSizes, columnSizes);
273 unsigned int blockRow = std::lower_bound(rowSizes->begin(), rowSizes->end(), row) - rowSizes->begin();
274 unsigned int blockColumn = std::lower_bound(columnSizes->begin(), columnSizes->end(), col) - columnSizes->begin();
275 const boost::shared_ptr<const InnerType> block = GetBlock(blockRow, blockColumn);
277 unsigned int actualRow = row;
280 actualRow = row-((*rowSizes)[blockRow-1])-1;
283 unsigned int actualCol = col;
284 if( blockColumn > 0 )
286 actualCol = col-((*columnSizes)[blockColumn-1])-1;
292 return (*block)(actualRow, actualCol);
296 return m_zeroElement;
300 template<
typename DataType,
typename InnerMatrixType>
301 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetStorageSize()
const
303 return m_storageSize;
306 template<
typename DataType,
typename InnerMatrixType>
307 MatrixStorage NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetType()
const
309 return m_storageType;
312 template<
typename DataType,
typename InnerMatrixType>
313 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfBlockRows()
const
315 if( this->GetTransposeFlag() ==
'N' )
317 return m_numberOfBlockRows;
321 return m_numberOfBlockColumns;
325 template<
typename DataType,
typename InnerMatrixType>
326 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfBlockColumns()
const
328 if( this->GetTransposeFlag() ==
'N' )
330 return m_numberOfBlockColumns;
334 return m_numberOfBlockRows;
338 template<
typename DataType,
typename InnerMatrixType>
339 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfRowsInBlockRow(
unsigned int blockRow)
const
341 if( this->GetTransposeFlag() ==
'N' )
343 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockRows, m_rowSizes);
347 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockColumns, m_columnSizes);
351 template<
typename DataType,
typename InnerMatrixType>
352 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfColumnsInBlockColumn(
unsigned int blockCol)
const
354 if( this->GetTransposeFlag() ==
'T' )
356 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockRows, m_rowSizes);
360 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockColumns, m_columnSizes);
364 template<
typename DataType,
typename InnerMatrixType>
365 typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>
::iterator NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::begin() {
return iterator(*
this, 0, 0); }
367 template<
typename DataType,
typename InnerMatrixType>
368 typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>
::iterator NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::end() {
return iterator(*
this); }
370 template<
typename DataType,
typename InnerMatrixType>
371 typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>
::const_iterator NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::begin()
const {
return const_iterator(*
this, 0, 0); }
373 template<
typename DataType,
typename InnerMatrixType>
374 typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>
::const_iterator NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::end()
const {
return const_iterator(*
this); }
377 template<
typename DataType,
typename InnerMatrixType>
378 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag> NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CreateWrapper(
const NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>& rhs)
380 return ThisType(rhs);
383 template<
typename DataType,
typename InnerMatrixType>
384 boost::shared_ptr<NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag> > NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CreateWrapper(
const boost::shared_ptr<NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag> >& rhs)
386 return boost::shared_ptr<ThisType>(
new ThisType(*rhs));
390 template<
typename DataType,
typename InnerMatrixType>
391 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfElementsInBlock(
unsigned int block,
unsigned int totalBlocks,
const Array<OneD, unsigned int>& sizes)
393 ASSERTL2(block < totalBlocks, std::string(
"Block Element ") + boost::lexical_cast<std::string>(block) +
394 std::string(
" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(totalBlocks) +
395 std::string(
" blocks."));
398 return sizes[block]+1;
402 return sizes[block] - sizes[block-1];
406 template<
typename DataType,
typename InnerMatrixType>
407 void NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::Initialize(
const unsigned int* rowsPerBlock,
const unsigned int* columnsPerBlock)
409 m_storageSize = this->GetRows()*this->GetColumns();
410 if (this->GetRows() > 0)
412 m_rowSizes[0] = rowsPerBlock[0] - 1;
413 for(
unsigned int i = 1; i < m_numberOfBlockRows; ++i)
415 m_rowSizes[i] = rowsPerBlock[i] + m_rowSizes[i-1];
418 if (this->GetColumns() > 0)
420 m_columnSizes[0] = columnsPerBlock[0] - 1;
421 for(
unsigned int i = 1; i < m_numberOfBlockColumns; ++i)
423 m_columnSizes[i] = columnsPerBlock[i] + m_columnSizes[i-1];
428 template<
typename DataType,
typename InnerMatrixType>
429 typename boost::call_traits<typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NumberType>::value_type
430 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_GetValue(
unsigned int row,
unsigned int column)
const
432 return (*
this)(row, column);
435 template<
typename DataType,
typename InnerMatrixType>
436 unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_GetStorageSize()
const
438 return this->GetStorageSize();
441 template<
typename DataType,
typename InnerMatrixType>
442 MatrixStorage NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_GetStorageType()
const
444 return this->GetType();
447 template<
typename DataType,
typename InnerMatrixType>
448 void NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_Transpose()
450 BOOST_FOREACH(boost::shared_ptr<InnerType> ptr, m_data)
521 template LIB_UTILITIES_EXPORT class NekMatrix< NekMatrix<NekDouble, StandardMatrixTag>, BlockMatrixTag>;
522 template LIB_UTILITIES_EXPORT class NekMatrix<NekMatrix< NekMatrix<NekDouble, StandardMatrixTag>, BlockMatrixTag>, BlockMatrixTag>;
523 template LIB_UTILITIES_EXPORT class NekMatrix<NekMatrix< NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>, BlockMatrixTag>;
524 template LIB_UTILITIES_EXPORT class NekMatrix<NekMatrix<NekMatrix< NekMatrix<NekDouble, StandardMatrixTag>, ScaledMatrixTag>, BlockMatrixTag>, BlockMatrixTag>;