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
35#include <boost/core/ignore_unused.hpp>
36
40
41namespace Nektar
42{
43template <typename DataType, typename InnerMatrixType>
45 MatrixStorage type)
46 : BaseType(0, 0, type), m_data(), m_rowSizes(), m_columnSizes(),
47 m_storageSize(), m_numberOfBlockRows(0), m_numberOfBlockColumns(0)
48{
49}
50
51template <typename DataType, typename InnerMatrixType>
53 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
54 unsigned int rowsPerBlock, unsigned int columnsPerBlock, MatrixStorage type)
55 : BaseType(numberOfBlockRows * rowsPerBlock,
56 numberOfBlockColumns * columnsPerBlock, type),
57 m_data(), m_rowSizes(numberOfBlockRows),
58 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
59 m_numberOfBlockRows(numberOfBlockRows),
60 m_numberOfBlockColumns(numberOfBlockColumns)
61{
62 m_storageSize = GetRequiredStorageSize();
64 m_storageSize, std::shared_ptr<InnerType>());
65 for (unsigned int i = 1; i <= numberOfBlockRows; ++i)
66 {
67 m_rowSizes[i - 1] = i * rowsPerBlock - 1;
68 }
69
70 for (unsigned int i = 1; i <= numberOfBlockColumns; ++i)
71 {
72 m_columnSizes[i - 1] = i * columnsPerBlock - 1;
73 }
74}
75
76template <typename DataType, typename InnerMatrixType>
78 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
79 const unsigned int *rowsPerBlock, const unsigned int *columnsPerBlock,
80 MatrixStorage type)
81 : BaseType(
82 std::accumulate(rowsPerBlock, rowsPerBlock + numberOfBlockRows, 0),
83 std::accumulate(columnsPerBlock,
84 columnsPerBlock + numberOfBlockColumns, 0),
85 type),
86 m_data(), m_rowSizes(numberOfBlockRows),
87 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
88 m_numberOfBlockRows(numberOfBlockRows),
89 m_numberOfBlockColumns(numberOfBlockColumns)
90{
91 m_storageSize = GetRequiredStorageSize();
93 m_storageSize, std::shared_ptr<InnerType>());
94 Initialize(rowsPerBlock, columnsPerBlock);
95}
96
97template <typename DataType, typename InnerMatrixType>
99 unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
100 const Array<OneD, const unsigned int> &rowsPerBlock,
101 const Array<OneD, const unsigned int> &columnsPerBlock, MatrixStorage type)
102 : BaseType(std::accumulate(rowsPerBlock.data(),
103 rowsPerBlock.data() + numberOfBlockRows, 0),
104 std::accumulate(columnsPerBlock.data(),
105 columnsPerBlock.data() + numberOfBlockColumns,
106 0),
107 type),
108 m_data(), m_rowSizes(numberOfBlockRows),
109 m_columnSizes(numberOfBlockColumns), m_storageSize(0),
110 m_numberOfBlockRows(numberOfBlockRows),
111 m_numberOfBlockColumns(numberOfBlockColumns)
112{
113 m_storageSize = GetRequiredStorageSize();
115 m_storageSize, std::shared_ptr<InnerType>());
116 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
117}
118
119template <typename DataType, typename InnerMatrixType>
121 const Array<OneD, const unsigned int> &rowsPerBlock,
122 const Array<OneD, const unsigned int> &columnsPerBlock, MatrixStorage type)
123 : BaseType(
124 std::accumulate(rowsPerBlock.begin(), rowsPerBlock.end(), 0),
125 std::accumulate(columnsPerBlock.begin(), columnsPerBlock.end(), 0),
126 type),
127 m_data(), m_rowSizes(rowsPerBlock.size()),
128 m_columnSizes(columnsPerBlock.size()), m_storageSize(0),
129 m_numberOfBlockRows(rowsPerBlock.size()),
130 m_numberOfBlockColumns(columnsPerBlock.size())
131{
132 m_storageSize = GetRequiredStorageSize();
134 m_storageSize, std::shared_ptr<InnerType>());
135 Initialize(rowsPerBlock.data(), columnsPerBlock.data());
136}
137
138template <typename DataType, typename InnerMatrixType>
139NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(
140 const NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag> &rhs)
141 : BaseType(rhs), m_data(rhs.m_data.size()), m_rowSizes(rhs.m_rowSizes),
142 m_columnSizes(rhs.m_columnSizes), m_storageSize(rhs.m_storageSize),
143 m_numberOfBlockRows(rhs.m_numberOfBlockRows),
144 m_numberOfBlockColumns(rhs.m_numberOfBlockColumns)
145{
146 for (unsigned int i = 0; i < rhs.m_data.size(); ++i)
147 {
148 m_data[i] = InnerType::CreateWrapper(rhs.m_data[i]);
149 }
150}
151
152template <typename DataType, typename InnerMatrixType>
153unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>,
154 BlockMatrixTag>::GetRequiredStorageSize() const
155{
156 return BaseType::GetRequiredStorageSize(this->GetStorageType(),
157 this->GetNumberOfBlockRows(),
158 this->GetNumberOfBlockColumns());
159}
160
161template <typename DataType, typename InnerMatrixType>
163 BlockMatrixTag>::CalculateBlockIndex(unsigned int row,
164 unsigned int column)
165 const
166{
167 return BaseType::CalculateIndex(this->GetStorageType(), row, column,
168 m_numberOfBlockRows, m_numberOfBlockColumns,
169 this->GetTransposeFlag());
170}
171
172template <typename DataType, typename InnerMatrixType>
174 BlockMatrixTag>::InnerType *
176 unsigned int row, unsigned int column) const
177{
178 ASSERTL2(
179 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
180 : row < m_numberOfBlockColumns,
181 std::string("Row ") + std::to_string(row) +
182 std::string(" requested in a block matrix with a maximum of ") +
183 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
184 ASSERTL2(
185 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
186 : column < m_numberOfBlockColumns,
187 std::string("Column ") + std::to_string(column) +
188 std::string(" requested in a block matrix with a maximum of ") +
189 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
190 int x = CalculateBlockIndex(row, column);
191 if (x == -1)
192 {
193 return 0;
194 }
195 else
196 {
197 return m_data[x].get();
198 }
199}
200
201template <typename DataType, typename InnerMatrixType>
202std::shared_ptr<const typename NekMatrix<NekMatrix<DataType, InnerMatrixType>,
203 BlockMatrixTag>::InnerType>
205 unsigned int row, unsigned int column) const
206{
207 ASSERTL2(
208 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
209 : row < m_numberOfBlockColumns,
210 std::string("Row ") + std::to_string(row) +
211 std::string(" requested in a block matrix with a maximum of ") +
212 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
213 ASSERTL2(
214 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
215 : column < m_numberOfBlockRows,
216 std::string("Column ") + std::to_string(column) +
217 std::string(" requested in a block matrix with a maximum of ") +
218 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
219 int x = CalculateBlockIndex(row, column);
220 if (x < 0)
221 {
222 return std::shared_ptr<const InnerType>();
223 }
224 else
225 {
226 return m_data[x];
227 }
228}
229
230template <typename DataType, typename InnerMatrixType>
231std::shared_ptr<typename NekMatrix<NekMatrix<DataType, InnerMatrixType>,
232 BlockMatrixTag>::InnerType>
234 unsigned int row, unsigned int column)
235{
236 ASSERTL2(
237 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
238 : row < m_numberOfBlockColumns,
239 std::string("Row ") + std::to_string(row) +
240 std::string(" requested in a block matrix with a maximum of ") +
241 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
242 ASSERTL2(
243 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
244 : column < m_numberOfBlockRows,
245 std::string("Column ") + std::to_string(column) +
246 std::string(" requested in a block matrix with a maximum of ") +
247 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
248 int x = CalculateBlockIndex(row, column);
249 if (x == -1)
250 {
251 return m_nullBlockPtr;
252 }
253 else
254 {
255 return m_data[x];
256 }
257}
258
259template <typename DataType, typename InnerMatrixType>
261 unsigned int row, unsigned int column, std::shared_ptr<InnerType> &m)
262{
263 ASSERTL2(
264 this->GetTransposeFlag() == 'N' ? row < m_numberOfBlockRows
265 : row < m_numberOfBlockColumns,
266 std::string("Row ") + std::to_string(row) +
267 std::string(" requested in a block matrix with a maximum of ") +
268 std::to_string(m_numberOfBlockRows) + std::string(" rows"));
269 ASSERTL2(
270 this->GetTransposeFlag() == 'N' ? column < m_numberOfBlockColumns
271 : column < m_numberOfBlockRows,
272 std::string("Column ") + std::to_string(column) +
273 std::string(" requested in a block matrix with a maximum of ") +
274 std::to_string(m_numberOfBlockColumns) + std::string(" columns"));
275 m_data[CalculateBlockIndex(row, column)] = InnerType::CreateWrapper(m);
276}
277
278template <typename DataType, typename InnerMatrixType>
280 BlockMatrixTag>::ConstGetValueType
282 unsigned int row, unsigned int col) const
283{
284 ASSERTL2(row < this->GetRows(),
285 std::string("Row ") + std::to_string(row) +
286 std::string(" requested in a matrix with a maximum of ") +
287 std::to_string(this->GetRows()) + std::string(" rows"));
288 ASSERTL2(col < this->GetColumns(),
289 std::string("Column ") + std::to_string(col) +
290 std::string(" requested in a matrix with a maximum of ") +
291 std::to_string(this->GetColumns()) + std::string(" columns"));
292
293 const Array<OneD, unsigned int> *rowSizes = &m_rowSizes;
294 const Array<OneD, unsigned int> *columnSizes = &m_columnSizes;
295
296 if (this->GetTransposeFlag() == 'T')
297 {
298 std::swap(rowSizes, columnSizes);
299 }
300
301 unsigned int blockRow =
302 std::lower_bound(rowSizes->begin(), rowSizes->end(), row) -
303 rowSizes->begin();
304 unsigned int blockColumn =
305 std::lower_bound(columnSizes->begin(), columnSizes->end(), col) -
306 columnSizes->begin();
307 const std::shared_ptr<const InnerType> block =
308 GetBlock(blockRow, blockColumn);
309
310 unsigned int actualRow = row;
311 if (blockRow > 0)
312 {
313 actualRow = row - ((*rowSizes)[blockRow - 1]) - 1;
314 }
315
316 unsigned int actualCol = col;
317 if (blockColumn > 0)
318 {
319 actualCol = col - ((*columnSizes)[blockColumn - 1]) - 1;
320 }
321
322 if (block)
323 {
324 return (*block)(actualRow, actualCol);
325 }
326 else
327 {
328 return m_zeroElement;
329 }
330}
331
332template <typename DataType, typename InnerMatrixType>
334 BlockMatrixTag>::GetStorageSize() const
335{
336 return m_storageSize;
337}
338
339template <typename DataType, typename InnerMatrixType>
341 BlockMatrixTag>::GetNumberOfBlockRows() const
342{
343 if (this->GetTransposeFlag() == 'N')
344 {
345 return m_numberOfBlockRows;
346 }
347 else
348 {
349 return m_numberOfBlockColumns;
350 }
351}
352
353template <typename DataType, typename InnerMatrixType>
355 BlockMatrixTag>::GetNumberOfBlockColumns() const
356{
357 if (this->GetTransposeFlag() == 'N')
358 {
359 return m_numberOfBlockColumns;
360 }
361 else
362 {
363 return m_numberOfBlockRows;
364 }
365}
366
367template <typename DataType, typename InnerMatrixType>
369 GetNumberOfRowsInBlockRow(unsigned int blockRow) const
370{
371 if (this->GetTransposeFlag() == 'N')
372 {
373 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockRows,
374 m_rowSizes);
375 }
376 else
377 {
378 return GetNumberOfElementsInBlock(blockRow, m_numberOfBlockColumns,
379 m_columnSizes);
380 }
381}
382
383template <typename DataType, typename InnerMatrixType>
385 GetNumberOfColumnsInBlockColumn(unsigned int blockCol) const
386{
387 if (this->GetTransposeFlag() == 'T')
388 {
389 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockRows,
390 m_rowSizes);
391 }
392 else
393 {
394 return GetNumberOfElementsInBlock(blockCol, m_numberOfBlockColumns,
395 m_columnSizes);
396 }
397}
398
399template <typename DataType, typename InnerMatrixType>
401 GetBlockSizes(Array<OneD, unsigned int> &rowSizes,
402 Array<OneD, unsigned int> &colSizes) const
403{
404 if (this->GetTransposeFlag() == 'T')
405 {
406 rowSizes = m_columnSizes;
407 colSizes = m_rowSizes;
408 }
409 else
410 {
411 rowSizes = m_rowSizes;
412 colSizes = m_columnSizes;
413 }
414}
415
416template <typename DataType, typename InnerMatrixType>
418 BlockMatrixTag>::iterator
420{
421 return iterator(*this, 0, 0);
422}
423
424template <typename DataType, typename InnerMatrixType>
426 BlockMatrixTag>::iterator
428{
429 return iterator(*this);
430}
431
432template <typename DataType, typename InnerMatrixType>
434 BlockMatrixTag>::const_iterator
436{
437 return const_iterator(*this, 0, 0);
438}
439
440template <typename DataType, typename InnerMatrixType>
442 BlockMatrixTag>::const_iterator
444{
445 return const_iterator(*this);
446}
447
448template <typename DataType, typename InnerMatrixType>
450 NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::
452 BlockMatrixTag> &rhs)
453{
454 return ThisType(rhs);
455}
456
457template <typename DataType, typename InnerMatrixType>
458std::shared_ptr<NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>>
459NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CreateWrapper(
460 const std::shared_ptr<
461 NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>> &rhs)
462{
463 return std::shared_ptr<ThisType>(new ThisType(*rhs));
464}
465
466template <typename DataType, typename InnerMatrixType>
468 GetNumberOfElementsInBlock(unsigned int block, unsigned int totalBlocks,
469 const Array<OneD, unsigned int> &sizes)
470{
471 boost::ignore_unused(totalBlocks);
472 ASSERTL2(block < totalBlocks,
473 std::string("Block Element ") + std::to_string(block) +
474 std::string(" requested in a matrix with a maximum of ") +
475 std::to_string(totalBlocks) + std::string(" blocks."));
476 if (block == 0)
477 {
478 return sizes[block] + 1;
479 }
480 else
481 {
482 return sizes[block] - sizes[block - 1];
483 }
484}
485
486template <typename DataType, typename InnerMatrixType>
488 BlockMatrixTag>::Initialize(const unsigned int *rowsPerBlock,
489 const unsigned int *columnsPerBlock)
490{
491 m_storageSize = this->GetRows() * this->GetColumns();
492 if (this->GetRows() > 0)
493 {
494 m_rowSizes[0] = rowsPerBlock[0] - 1;
495 for (unsigned int i = 1; i < m_numberOfBlockRows; ++i)
496 {
497 m_rowSizes[i] = rowsPerBlock[i] + m_rowSizes[i - 1];
498 }
499 }
500 if (this->GetColumns() > 0)
501 {
502 m_columnSizes[0] = columnsPerBlock[0] - 1;
503 for (unsigned int i = 1; i < m_numberOfBlockColumns; ++i)
504 {
505 m_columnSizes[i] = columnsPerBlock[i] + m_columnSizes[i - 1];
506 }
507 }
508}
509
510template <typename DataType, typename InnerMatrixType>
511typename boost::call_traits<
513 BlockMatrixTag>::NumberType>::value_type
515 unsigned int row, unsigned int column) const
516{
517 return (*this)(row, column);
518}
519
520template <typename DataType, typename InnerMatrixType>
522 BlockMatrixTag>::v_GetStorageSize() const
523{
524 return this->GetStorageSize();
525}
526
527template <typename DataType, typename InnerMatrixType>
529 BlockMatrixTag>::v_Transpose()
530{
531 for (auto &ptr : m_data)
532 {
533 if (ptr.get() != 0)
534 {
535 ptr->Transpose();
536 }
537 }
538}
539
540template LIB_UTILITIES_EXPORT class NekMatrix<
542template LIB_UTILITIES_EXPORT class NekMatrix<
544 BlockMatrixTag>;
545template LIB_UTILITIES_EXPORT class NekMatrix<
547 BlockMatrixTag>;
548template LIB_UTILITIES_EXPORT class NekMatrix<
549 NekMatrix<
551 BlockMatrixTag>,
552 BlockMatrixTag>;
553
554template LIB_UTILITIES_EXPORT class NekMatrix<
556template LIB_UTILITIES_EXPORT class NekMatrix<
558 BlockMatrixTag>;
559template LIB_UTILITIES_EXPORT class NekMatrix<
561 BlockMatrixTag>;
562template LIB_UTILITIES_EXPORT class NekMatrix<
563 NekMatrix<
565 BlockMatrixTag>,
566 BlockMatrixTag>;
567} // 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:272
#define LIB_UTILITIES_EXPORT
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2