Nektar++
BlockMatrix.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: BlockMatrix.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33///////////////////////////////////////////////////////////////////////////////
34
38
39namespace Nektar
40{
41template <typename DataType, typename InnerMatrixType>
43 MatrixStorage type)
44 : BaseType(0, 0, type), m_data(), m_rowSizes(), m_columnSizes(),
45 m_storageSize(), m_numberOfBlockRows(0), m_numberOfBlockColumns(0)
46{
47}
48
49template <typename DataType, typename InnerMatrixType>
51 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
52 unsigned int rowsPerBlock, unsigned int columnsPerBlock, MatrixStorage type)
53 : BaseType(numberOfBlockRows * rowsPerBlock,
54 numberOfBlockColumns * columnsPerBlock, type),
55 m_data(), m_rowSizes(numberOfBlockRows),
56 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
57 m_numberOfBlockRows(numberOfBlockRows),
58 m_numberOfBlockColumns(numberOfBlockColumns)
59{
60 m_storageSize = GetRequiredStorageSize();
62 m_storageSize, std::shared_ptr<InnerType>());
63 for (unsigned int i = 1; i <= numberOfBlockRows; ++i)
64 {
65 m_rowSizes[i - 1] = i * rowsPerBlock - 1;
66 }
67
68 for (unsigned int i = 1; i <= numberOfBlockColumns; ++i)
69 {
70 m_columnSizes[i - 1] = i * columnsPerBlock - 1;
71 }
72}
73
74template <typename DataType, typename InnerMatrixType>
76 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
77 const unsigned int *rowsPerBlock, const unsigned int *columnsPerBlock,
78 MatrixStorage type)
79 : BaseType(
80 std::accumulate(rowsPerBlock, rowsPerBlock + numberOfBlockRows, 0),
81 std::accumulate(columnsPerBlock,
82 columnsPerBlock + numberOfBlockColumns, 0),
83 type),
84 m_data(), m_rowSizes(numberOfBlockRows),
85 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
86 m_numberOfBlockRows(numberOfBlockRows),
87 m_numberOfBlockColumns(numberOfBlockColumns)
88{
89 m_storageSize = GetRequiredStorageSize();
91 m_storageSize, std::shared_ptr<InnerType>());
92 Initialize(rowsPerBlock, columnsPerBlock);
93}
94
95template <typename DataType, typename InnerMatrixType>
97 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
98 const Array<OneD, const unsigned int> &rowsPerBlock,
99 const Array<OneD, const unsigned int> &columnsPerBlock, MatrixStorage type)
100 : BaseType(std::accumulate(rowsPerBlock.data(),
101 rowsPerBlock.data() + numberOfBlockRows, 0),
102 std::accumulate(columnsPerBlock.data(),
103 columnsPerBlock.data() + numberOfBlockColumns,
104 0),
105 type),
106 m_data(), m_rowSizes(numberOfBlockRows),
107 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
108 m_numberOfBlockRows(numberOfBlockRows),
109 m_numberOfBlockColumns(numberOfBlockColumns)
110{
111 m_storageSize = GetRequiredStorageSize();
113 m_storageSize, std::shared_ptr<InnerType>());
114 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
115}
116
117template <typename DataType, typename InnerMatrixType>
119 const Array<OneD, const unsigned int> &rowsPerBlock,
120 const Array<OneD, const unsigned int> &columnsPerBlock, MatrixStorage type)
121 : BaseType(
122 std::accumulate(rowsPerBlock.begin(), rowsPerBlock.end(), 0),
123 std::accumulate(columnsPerBlock.begin(), columnsPerBlock.end(), 0),
124 type),
125 m_data(), m_rowSizes(rowsPerBlock.size()),
126 m_columnSizes(columnsPerBlock.size()), m_storageSize(0),
127 m_numberOfBlockRows(rowsPerBlock.size()),
128 m_numberOfBlockColumns(columnsPerBlock.size())
129{
130 m_storageSize = GetRequiredStorageSize();
132 m_storageSize, std::shared_ptr<InnerType>());
133 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
134}
135
136template <typename DataType, typename InnerMatrixType>
137NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
138 const NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag> &rhs)
139 : BaseType(rhs), m_data(rhs.m_data.size()), m_rowSizes(rhs.m_rowSizes),
140 m_columnSizes(rhs.m_columnSizes), m_storageSize(rhs.m_storageSize),
141 m_numberOfBlockRows(rhs.m_numberOfBlockRows),
142 m_numberOfBlockColumns(rhs.m_numberOfBlockColumns)
143{
144 for (unsigned int i = 0; i < rhs.m_data.size(); ++i)
145 {
146 m_data[i] = InnerType::CreateWrapper(rhs.m_data[i]);
147 }
148}
149
150template <typename DataType, typename InnerMatrixType>
151unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>,
152 BlockMatrixTag>::GetRequiredStorageSize() const
153{
154 return BaseType::GetRequiredStorageSize(this->GetStorageType(),
155 this->GetNumberOfBlockRows(),
156 this->GetNumberOfBlockColumns());
157}
158
159template <typename DataType, typename InnerMatrixType>
161 BlockMatrixTag>::CalculateBlockIndex(unsigned int row,
162 unsigned int column)
163 const
164{
165 return BaseType::CalculateIndex(this->GetStorageType(), row, column,
166 m_numberOfBlockRows, m_numberOfBlockColumns,
167 this->GetTransposeFlag());
168}
169
170template <typename DataType, typename InnerMatrixType>
172 BlockMatrixTag>::InnerType *
174 unsigned int row, unsigned int column) const
175{
176 ASSERTL2(
177 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
178 : row < m_numberOfBlockColumns,
179 std::string("Row ") + std::to_string(row) +
180 std::string(" requested in a block matrix with a maximum of ") +
181 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
182 ASSERTL2(
183 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
184 : column < m_numberOfBlockColumns,
185 std::string("Column ") + std::to_string(column) +
186 std::string(" requested in a block matrix with a maximum of ") +
187 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
188 int x = CalculateBlockIndex(row, column);
189 if (x == -1)
190 {
191 return nullptr;
192 }
193 else
194 {
195 return m_data[x].get();
196 }
197}
198
199template <typename DataType, typename InnerMatrixType>
200std::shared_ptr<const typename NekMatrix<NekMatrix<DataType, InnerMatrixType>,
201 BlockMatrixTag>::InnerType>
203 unsigned int row, unsigned int column) const
204{
205 ASSERTL2(
206 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
207 : row < m_numberOfBlockColumns,
208 std::string("Row ") + std::to_string(row) +
209 std::string(" requested in a block matrix with a maximum of ") +
210 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
211 ASSERTL2(
212 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
213 : column < m_numberOfBlockRows,
214 std::string("Column ") + std::to_string(column) +
215 std::string(" requested in a block matrix with a maximum of ") +
216 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
217 int x = CalculateBlockIndex(row, column);
218 if (x < 0)
219 {
220 return std::shared_ptr<const InnerType>();
221 }
222 else
223 {
224 return m_data[x];
225 }
226}
227
228template <typename DataType, typename InnerMatrixType>
229std::shared_ptr<typename NekMatrix<NekMatrix<DataType, InnerMatrixType>,
230 BlockMatrixTag>::InnerType> &
232 unsigned int row, unsigned int column)
233{
234 ASSERTL2(
235 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
236 : row < m_numberOfBlockColumns,
237 std::string("Row ") + std::to_string(row) +
238 std::string(" requested in a block matrix with a maximum of ") +
239 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
240 ASSERTL2(
241 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
242 : column < m_numberOfBlockRows,
243 std::string("Column ") + std::to_string(column) +
244 std::string(" requested in a block matrix with a maximum of ") +
245 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
246 int x = CalculateBlockIndex(row, column);
247 if (x == -1)
248 {
249 return m_nullBlockPtr;
250 }
251 else
252 {
253 return m_data[x];
254 }
255}
256
257template <typename DataType, typename InnerMatrixType>
259 unsigned int row, unsigned int column, std::shared_ptr<InnerType> &m)
260{
261 ASSERTL2(
262 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
263 : row < m_numberOfBlockColumns,
264 std::string("Row ") + std::to_string(row) +
265 std::string(" requested in a block matrix with a maximum of ") +
266 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
267 ASSERTL2(
268 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
269 : column < m_numberOfBlockRows,
270 std::string("Column ") + std::to_string(column) +
271 std::string(" requested in a block matrix with a maximum of ") +
272 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
273 m_data[CalculateBlockIndex(row, column)] = InnerType::CreateWrapper(m);
274}
275
276template <typename DataType, typename InnerMatrixType>
278 BlockMatrixTag>::ConstGetValueType
280 unsigned int row, unsigned int col) const
281{
282 ASSERTL2(row < this->GetRows(),
283 std::string("Row ") + std::to_string(row) +
284 std::string(" requested in a matrix with a maximum of ") +
285 std::to_string(this->GetRows()) + std::string(" rows"));
286 ASSERTL2(col < this->GetColumns(),
287 std::string("Column ") + std::to_string(col) +
288 std::string(" requested in a matrix with a maximum of ") +
289 std::to_string(this->GetColumns()) + std::string(" columns"));
290
291 const Array<OneD, unsigned int> *rowSizes = &m_rowSizes;
292 const Array<OneD, unsigned int> *columnSizes = &m_columnSizes;
293
294 if (this->GetTransposeFlag() == 'T')
295 {
296 std::swap(rowSizes, columnSizes);
297 }
298
299 unsigned int blockRow =
300 std::lower_bound(rowSizes->begin(), rowSizes->end(), row) -
301 rowSizes->begin();
302 unsigned int blockColumn =
303 std::lower_bound(columnSizes->begin(), columnSizes->end(), col) -
304 columnSizes->begin();
305 const std::shared_ptr<const InnerType> block =
306 GetBlock(blockRow, blockColumn);
307
308 unsigned int actualRow = row;
309 if (blockRow > 0)
310 {
311 actualRow = row - ((*rowSizes)[blockRow - 1]) - 1;
312 }
313
314 unsigned int actualCol = col;
315 if (blockColumn > 0)
316 {
317 actualCol = col - ((*columnSizes)[blockColumn - 1]) - 1;
318 }
319
320 if (block)
321 {
322 return (*block)(actualRow, actualCol);
323 }
324 else
325 {
326 return m_zeroElement;
327 }
328}
329
330template <typename DataType, typename InnerMatrixType>
332 BlockMatrixTag>::GetStorageSize() const
333{
334 return m_storageSize;
335}
336
337template <typename DataType, typename InnerMatrixType>
339 BlockMatrixTag>::GetNumberOfBlockRows() const
340{
341 if (this->GetTransposeFlag() == 'N')
342 {
343 return m_numberOfBlockRows;
344 }
345 else
346 {
347 return m_numberOfBlockColumns;
348 }
349}
350
351template <typename DataType, typename InnerMatrixType>
353 BlockMatrixTag>::GetNumberOfBlockColumns() const
354{
355 if (this->GetTransposeFlag() == 'N')
356 {
357 return m_numberOfBlockColumns;
358 }
359 else
360 {
361 return m_numberOfBlockRows;
362 }
363}
364
365template <typename DataType, typename InnerMatrixType>
367 GetNumberOfRowsInBlockRow(unsigned int blockRow) const
368{
369 if (this->GetTransposeFlag() == 'N')
370 {
371 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockRows,
372 m_rowSizes);
373 }
374 else
375 {
376 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockColumns,
377 m_columnSizes);
378 }
379}
380
381template <typename DataType, typename InnerMatrixType>
383 GetNumberOfColumnsInBlockColumn(unsigned int blockCol) const
384{
385 if (this->GetTransposeFlag() == 'T')
386 {
387 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockRows,
388 m_rowSizes);
389 }
390 else
391 {
392 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockColumns,
393 m_columnSizes);
394 }
395}
396
397template <typename DataType, typename InnerMatrixType>
399 GetBlockSizes(Array<OneD, unsigned int> &rowSizes,
400 Array<OneD, unsigned int> &colSizes) const
401{
402 if (this->GetTransposeFlag() == 'T')
403 {
404 rowSizes = m_columnSizes;
405 colSizes = m_rowSizes;
406 }
407 else
408 {
409 rowSizes = m_rowSizes;
410 colSizes = m_columnSizes;
411 }
412}
413
414template <typename DataType, typename InnerMatrixType>
416 BlockMatrixTag>::iterator
418{
419 return iterator(*this, 0, 0);
420}
421
422template <typename DataType, typename InnerMatrixType>
424 BlockMatrixTag>::iterator
426{
427 return iterator(*this);
428}
429
430template <typename DataType, typename InnerMatrixType>
432 BlockMatrixTag>::const_iterator
434{
435 return const_iterator(*this, 0, 0);
436}
437
438template <typename DataType, typename InnerMatrixType>
440 BlockMatrixTag>::const_iterator
442{
443 return const_iterator(*this);
444}
445
446template <typename DataType, typename InnerMatrixType>
448 NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::
450 BlockMatrixTag> &rhs)
451{
452 return ThisType(rhs);
453}
454
455template <typename DataType, typename InnerMatrixType>
456std::shared_ptr<NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>>
457NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CreateWrapper(
458 const std::shared_ptr<
459 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>> &rhs)
460{
461 return std::shared_ptr<ThisType>(new ThisType(*rhs));
462}
463
464template <typename DataType, typename InnerMatrixType>
466 GetNumberOfElementsInBlock(unsigned int block,
467 [[maybe_unused]] unsigned int totalBlocks,
468 const Array<OneD, unsigned int> &sizes)
469{
470 ASSERTL2(block < totalBlocks,
471 std::string("Block Element ") + std::to_string(block) +
472 std::string(" requested in a matrix with a maximum of ") +
473 std::to_string(totalBlocks) + std::string(" blocks."));
474 if (block == 0)
475 {
476 return sizes[block] + 1;
477 }
478 else
479 {
480 return sizes[block] - sizes[block - 1];
481 }
482}
483
484template <typename DataType, typename InnerMatrixType>
486 BlockMatrixTag>::Initialize(const unsigned int *rowsPerBlock,
487 const unsigned int *columnsPerBlock)
488{
489 m_storageSize = this->GetRows() * this->GetColumns();
490 if (this->GetRows() > 0)
491 {
492 m_rowSizes[0] = rowsPerBlock[0] - 1;
493 for (unsigned int i = 1; i < m_numberOfBlockRows; ++i)
494 {
495 m_rowSizes[i] = rowsPerBlock[i] + m_rowSizes[i - 1];
496 }
497 }
498 if (this->GetColumns() > 0)
499 {
500 m_columnSizes[0] = columnsPerBlock[0] - 1;
501 for (unsigned int i = 1; i < m_numberOfBlockColumns; ++i)
502 {
503 m_columnSizes[i] = columnsPerBlock[i] + m_columnSizes[i - 1];
504 }
505 }
506}
507
508template <typename DataType, typename InnerMatrixType>
509typename boost::call_traits<
511 BlockMatrixTag>::NumberType>::value_type
513 unsigned int row, unsigned int column) const
514{
515 return (*this)(row, column);
516}
517
518template <typename DataType, typename InnerMatrixType>
520 BlockMatrixTag>::v_GetStorageSize() const
521{
522 return this->GetStorageSize();
523}
524
525template <typename DataType, typename InnerMatrixType>
527 BlockMatrixTag>::v_Transpose()
528{
529 for (auto &ptr : m_data)
530 {
531 if (ptr.get() != nullptr)
532 {
533 ptr->Transpose();
534 }
535 }
536}
537
538template LIB_UTILITIES_EXPORT class NekMatrix<
540template LIB_UTILITIES_EXPORT class NekMatrix<
542 BlockMatrixTag>;
543template LIB_UTILITIES_EXPORT class NekMatrix<
545 BlockMatrixTag>;
546template LIB_UTILITIES_EXPORT class NekMatrix<
547 NekMatrix<
549 BlockMatrixTag>,
550 BlockMatrixTag>;
551
552template LIB_UTILITIES_EXPORT class NekMatrix<
554template LIB_UTILITIES_EXPORT class NekMatrix<
556 BlockMatrixTag>;
557template LIB_UTILITIES_EXPORT class NekMatrix<
559 BlockMatrixTag>;
560template LIB_UTILITIES_EXPORT class NekMatrix<
561 NekMatrix<
563 BlockMatrixTag>,
564 BlockMatrixTag>;
565} // namespace Nektar
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
#define LIB_UTILITIES_EXPORT
STL namespace.