Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_BandedContainer_def.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#ifndef IFPACK2_BANDEDCONTAINER_DEF_HPP
44#define IFPACK2_BANDEDCONTAINER_DEF_HPP
45
46#include "Teuchos_LAPACK.hpp"
47#include "Tpetra_CrsMatrix.hpp"
48#include <iostream>
49#include <sstream>
50
51#ifdef HAVE_MPI
52# include <mpi.h>
53# include "Teuchos_DefaultMpiComm.hpp"
54#else
55# include "Teuchos_DefaultSerialComm.hpp"
56#endif // HAVE_MPI
57
58namespace Ifpack2 {
59
60template<class MatrixType, class LocalScalarType>
62BandedContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
63 const Teuchos::Array<Teuchos::Array<LO> >& partitions,
64 const Teuchos::RCP<const import_type>&,
65 bool pointIndexed) :
66 ContainerImpl<MatrixType, LocalScalarType>(matrix, partitions, pointIndexed),
67 ipiv_(this->blockRows_.size() * this->scalarsPerRow_),
68 kl_(this->numBlocks_, -1),
69 ku_(this->numBlocks_, -1),
70 scalarOffsets_(this->numBlocks_)
71{
72 TEUCHOS_TEST_FOR_EXCEPTION(
73 ! matrix->hasColMap (), std::invalid_argument, "Ifpack2::BandedContainer: "
74 "The constructor's input matrix must have a column Map.");
75}
76
77template<class MatrixType, class LocalScalarType>
80
81template<class MatrixType, class LocalScalarType>
83setParameters (const Teuchos::ParameterList& List)
84{
85 if(List.isParameter("relaxation: banded container superdiagonals"))
86 {
87 int ku = List.get<int>("relaxation: banded container superdiagonals");
88 for(LO b = 0; b < this->numBlocks_; b++)
89 ku_[b] = ku;
90 }
91 if(List.isParameter("relaxation: banded container subdiagonals"))
92 {
93 int kl = List.get<int>("relaxation: banded container subdiagonals");
94 for(LO b = 0; b < this->numBlocks_; b++)
95 kl_[b] = kl;
96 }
97}
98
99template<class MatrixType, class LocalScalarType>
102{
103 using Teuchos::Array;
104 using Teuchos::ArrayView;
105 //now, for any block where kl_ or ku_ has not already been set, compute the actual bandwidth
106 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
107 size_t colToOffsetSize = this->inputMatrix_->getLocalNumCols();
108 if(this->pointIndexed_)
109 colToOffsetSize *= this->bcrsBlockSize_;
110 Array<LO> colToBlockOffset(colToOffsetSize, INVALID);
111 //Same logic as extract() to find entries in blocks efficiently
112 //(but here, just to find bandwidth, no scalars used)
113 for(int i = 0; i < this->numBlocks_; i++)
114 {
115 //maxSub, maxSuper are the maximum lower and upper bandwidth
116 LO maxSub = 0;
117 LO maxSuper = 0;
118 if(this->scalarsPerRow_ > 1)
119 {
120 //Get the interval where block i is defined in blockRows_
121 LO blockStart = this->blockOffsets_[i];
122 LO blockEnd = (i == this->numBlocks_ - 1) ? this->blockRows_.size() : this->blockOffsets_[i + 1];
123 ArrayView<const LO> blockRows = this->getBlockRows(i);
124 //Set the lookup table entries for the columns appearing in block i.
125 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
126 //this is OK. The values updated here are only needed to process block i's entries.
127 for(size_t j = 0; j < size_t(blockRows.size()); j++)
128 {
129 LO localCol = this->translateRowToCol(blockRows[j]);
130 colToBlockOffset[localCol] = blockStart + j;
131 }
132
133 using h_inds_type = typename block_crs_matrix_type::local_inds_host_view_type;
134 using h_vals_type = typename block_crs_matrix_type::values_host_view_type;
135 for(LO blockRow = 0; blockRow < LO(blockRows.size()); blockRow++)
136 {
137 //get a raw view of the whole block row
138 h_inds_type indices;
139 h_vals_type values;
140 LO inputRow = this->blockRows_[blockStart + blockRow];
141 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
142 LO numEntries = (LO) indices.size();
143 for(LO k = 0; k < numEntries; k++)
144 {
145 LO colOffset = colToBlockOffset[indices[k]];
146 if(blockStart <= colOffset && colOffset < blockEnd)
147 {
148 //This entry does appear in the diagonal block.
149 //(br, bc) identifies the scalar's position in the BlockCrs block.
150 //Convert this to (r, c) which is its position in the container block.
151 LO blockCol = colOffset - blockStart;
152 for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
153 {
154 for(LO br = 0; br < this->bcrsBlockSize_; br++)
155 {
156 LO r = this->bcrsBlockSize_ * blockRow + br;
157 LO c = this->bcrsBlockSize_ * blockCol + bc;
158 if(r - c > maxSub)
159 maxSub = r - c;
160 if(c - r > maxSuper)
161 maxSuper = c - r;
162 }
163 }
164 }
165 }
166 }
167 }
168 else
169 {
170 //Get the interval where block i is defined in blockRows_
171 LO blockStart = this->blockOffsets_[i];
172 LO blockEnd = (i == this->numBlocks_ - 1) ? this->blockRows_.size() : this->blockOffsets_[i + 1];
173 ArrayView<const LO> blockRows = this->getBlockRows(i);
174 //Set the lookup table entries for the columns appearing in block i.
175 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
176 //this is OK. The values updated here are only needed to process block i's entries.
177 for(size_t j = 0; j < size_t(blockRows.size()); j++)
178 {
179 //translateRowToCol will return the corresponding split column
180 LO localCol = this->translateRowToCol(blockRows[j]);
181 colToBlockOffset[localCol] = blockStart + j;
182 }
183 for(LO blockRow = 0; blockRow < LO(blockRows.size()); blockRow++)
184 {
185 //get a view of the general row
186 LO inputSplitRow = this->blockRows_[blockStart + blockRow];
187 auto rowView = this->getInputRowView(inputSplitRow);
188 for(size_t k = 0; k < rowView.size(); k++)
189 {
190 LO colOffset = colToBlockOffset[rowView.ind(k)];
191 if(blockStart <= colOffset && colOffset < blockEnd)
192 {
193 LO blockCol = colOffset - blockStart;
194 maxSub = std::max(maxSub, blockRow - blockCol);
195 maxSuper = std::max(maxSuper, blockCol - blockRow);
196 }
197 }
198 }
199 }
200 kl_[i] = maxSub;
201 ku_[i] = maxSuper;
202 }
203}
204
205template<class MatrixType, class LocalScalarType>
206void
208initialize ()
209{
210 //note: in BlockRelaxation, Container_->setParameters() immediately
211 //precedes Container_->initialize(). So no matter what, either all
212 //the block bandwidths were set (in setParameters()) or none of them were.
213 //If none were they must be computed individually.
214 if(kl_[0] == -1)
215 computeBandwidth();
216 GO totalScalars = 0;
217 for(LO b = 0; b < this->numBlocks_; b++)
218 {
219 LO stride = 2 * kl_[b] + ku_[b] + 1;
220 scalarOffsets_[b] = totalScalars;
221 totalScalars += stride * this->blockSizes_[b] * this->scalarsPerRow_;
222 }
223 scalars_.resize(totalScalars);
224 for(int b = 0; b < this->numBlocks_; b++)
225 {
226 //NOTE: the stride and upper bandwidth used to construct the SerialBandDenseMatrix looks
227 //too large, but the extra kl_ in upper band space is needed by the LAPACK factorization routine.
228 LO nrows = this->blockSizes_[b] * this->scalarsPerRow_;
229 diagBlocks_.emplace_back(Teuchos::View, scalars_.data() + scalarOffsets_[b], 2 * kl_[b] + ku_[b] + 1, nrows, nrows, kl_[b], kl_[b] + ku_[b]);
230 diagBlocks_[b].putScalar(Teuchos::ScalarTraits<LSC>::zero());
231 }
232 std::fill (ipiv_.begin (), ipiv_.end (), 0);
233 // We assume that if you called this method, you intend to recompute
234 // everything.
235 this->IsComputed_ = false;
236 this->IsInitialized_ = true;
237}
238
239template<class MatrixType, class LocalScalarType>
240void
242compute ()
243{
244 TEUCHOS_TEST_FOR_EXCEPTION(
245 ipiv_.size () != this->blockRows_.size() * this->scalarsPerRow_, std::logic_error,
246 "Ifpack2::BandedContainer::compute: ipiv_ array has the wrong size. "
247 "Please report this bug to the Ifpack2 developers.");
248
249 this->IsComputed_ = false;
250 if (!this->isInitialized ()) {
251 this->initialize ();
252 }
253
254 extract (); // Extract the submatrices
255 factor (); // Factor them
256
257 this->IsComputed_ = true;
258}
259
260template<class MatrixType, class LocalScalarType>
262{
263 using Teuchos::Array;
264 using Teuchos::ArrayView;
265 using STS = Teuchos::ScalarTraits<SC>;
266 SC zero = STS::zero();
267 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
268 //To extract diagonal blocks, need to translate local rows to local columns.
269 //Strategy: make a lookup table that translates local cols in the matrix to offsets in blockRows_:
270 //blockOffsets_[b] <= offset < blockOffsets_[b+1]: tests whether the column is in block b.
271 //offset - blockOffsets_[b]: gives the column within block b.
272 //
273 //This provides the block and col within a block in O(1).
274 if(this->scalarsPerRow_ > 1)
275 {
276 Array<LO> colToBlockOffset(this->inputBlockMatrix_->getLocalNumCols(), INVALID);
277 for(int i = 0; i < this->numBlocks_; i++)
278 {
279 //Get the interval where block i is defined in blockRows_
280 LO blockStart = this->blockOffsets_[i];
281 LO blockEnd = blockStart + this->blockSizes_[i];
282 ArrayView<const LO> blockRows = this->getBlockRows(i);
283 //Set the lookup table entries for the columns appearing in block i.
284 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
285 //this is OK. The values updated here are only needed to process block i's entries.
286 for(size_t j = 0; j < size_t(blockRows.size()); j++)
287 {
288 LO localCol = this->translateRowToCol(blockRows[j]);
289 colToBlockOffset[localCol] = blockStart + j;
290 }
291 using h_inds_type = typename block_crs_matrix_type::local_inds_host_view_type;
292 using h_vals_type = typename block_crs_matrix_type::values_host_view_type;
293 for(LO blockRow = 0; blockRow < LO(blockRows.size()); blockRow++)
294 {
295 //get a raw view of the whole block row
296 h_inds_type indices;
297 h_vals_type values;
298 LO inputRow = this->blockRows_[blockStart + blockRow];
299 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
300 LO numEntries = (LO) indices.size();
301 for(LO k = 0; k < numEntries; k++)
302 {
303 LO colOffset = colToBlockOffset[indices[k]];
304 if(blockStart <= colOffset && colOffset < blockEnd)
305 {
306 //This entry does appear in the diagonal block.
307 //(br, bc) identifies the scalar's position in the BlockCrs block.
308 //Convert this to (r, c) which is its position in the container block.
309 LO blockCol = colOffset - blockStart;
310 for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
311 {
312 for(LO br = 0; br < this->bcrsBlockSize_; br++)
313 {
314 LO r = this->bcrsBlockSize_ * blockRow + br;
315 LO c = this->bcrsBlockSize_ * blockCol + bc;
316 auto val = values[k * (this->bcrsBlockSize_ * this->bcrsBlockSize_) + (br + this->bcrsBlockSize_ * bc)];
317 if(val != zero)
318 diagBlocks_[i](r, c) = val;
319 }
320 }
321 }
322 }
323 }
324 }
325 }
326 else
327 {
328 //get the mapping from point-indexed matrix columns to offsets in blockRows_
329 //(this includes regular CrsMatrix columns, in which case bcrsBlockSize_ == 1)
330 Array<LO> colToBlockOffset(this->inputMatrix_->getLocalNumCols() * this->bcrsBlockSize_, INVALID);
331 for(int i = 0; i < this->numBlocks_; i++)
332 {
333 //Get the interval where block i is defined in blockRows_
334 LO blockStart = this->blockOffsets_[i];
335 LO blockEnd = blockStart + this->blockSizes_[i];
336 ArrayView<const LO> blockRows = this->getBlockRows(i);
337 //Set the lookup table entries for the columns appearing in block i.
338 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
339 //this is OK. The values updated here are only needed to process block i's entries.
340 for(size_t j = 0; j < size_t(blockRows.size()); j++)
341 {
342 //translateRowToCol will return the corresponding split column
343 LO localCol = this->translateRowToCol(blockRows[j]);
344 colToBlockOffset[localCol] = blockStart + j;
345 }
346 for(size_t blockRow = 0; blockRow < size_t(blockRows.size()); blockRow++)
347 {
348 //get a view of the split row
349 LO inputPointRow = this->blockRows_[blockStart + blockRow];
350 auto rowView = this->getInputRowView(inputPointRow);
351 for(size_t k = 0; k < rowView.size(); k++)
352 {
353 LO colOffset = colToBlockOffset[rowView.ind(k)];
354 if(blockStart <= colOffset && colOffset < blockEnd)
355 {
356 LO blockCol = colOffset - blockStart;
357 auto val = rowView.val(k);
358 if(val != zero)
359 diagBlocks_[i](blockRow, blockCol) = rowView.val(k);
360 }
361 }
362 }
363 }
364 }
365}
366
367template<class MatrixType, class LocalScalarType>
368void
369BandedContainer<MatrixType, LocalScalarType>::
370clearBlocks ()
371{
372 diagBlocks_.clear();
373 scalars_.clear();
374 ContainerImpl<MatrixType, LocalScalarType>::clearBlocks();
375}
376
377template<class MatrixType, class LocalScalarType>
378void
379BandedContainer<MatrixType, LocalScalarType>::
380factor ()
381{
382 Teuchos::LAPACK<int, LSC> lapack;
383 int INFO = 0;
384
385 for(int i = 0; i < this->numBlocks_; i++)
386 {
387 // Plausibility checks for matrix
388 TEUCHOS_TEST_FOR_EXCEPTION(diagBlocks_[i].values()==0, std::invalid_argument,
389 "BandedContainer<T>::factor: Diagonal block is an empty SerialBandDenseMatrix<T>!");
390 TEUCHOS_TEST_FOR_EXCEPTION(diagBlocks_[i].upperBandwidth() < diagBlocks_[i].lowerBandwidth(), std::invalid_argument,
391 "BandedContainer<T>::factor: Diagonal block needs kl additional superdiagonals for factorization! However, the number of superdiagonals is smaller than the number of subdiagonals!");
392 int* blockIpiv = &ipiv_[this->blockOffsets_[i] * this->scalarsPerRow_];
393 lapack.GBTRF (diagBlocks_[i].numRows(),
394 diagBlocks_[i].numCols(),
395 diagBlocks_[i].lowerBandwidth(),
396 diagBlocks_[i].upperBandwidth() - diagBlocks_[i].lowerBandwidth(), /* enter the real number of superdiagonals (see Teuchos_SerialBandDenseSolver)*/
397 diagBlocks_[i].values(),
398 diagBlocks_[i].stride(),
399 blockIpiv,
400 &INFO);
401
402 // INFO < 0 is a bug.
403 TEUCHOS_TEST_FOR_EXCEPTION(
404 INFO < 0, std::logic_error, "Ifpack2::BandedContainer::factor: "
405 "LAPACK's _GBTRF (LU factorization with partial pivoting) was called "
406 "incorrectly. INFO = " << INFO << " < 0. "
407 "Please report this bug to the Ifpack2 developers.");
408 // INFO > 0 means the matrix is singular. This is probably an issue
409 // either with the choice of rows the rows we extracted, or with the
410 // input matrix itself.
411 TEUCHOS_TEST_FOR_EXCEPTION(
412 INFO > 0, std::runtime_error, "Ifpack2::BandedContainer::factor: "
413 "LAPACK's _GBTRF (LU factorization with partial pivoting) reports that the "
414 "computed U factor is exactly singular. U(" << INFO << "," << INFO << ") "
415 "(one-based index i) is exactly zero. This probably means that the input "
416 "matrix has a singular diagonal block.");
417 }
418}
419
420template<class MatrixType, class LocalScalarType>
421void
422BandedContainer<MatrixType, LocalScalarType>::
423solveBlock(ConstHostSubviewLocal X,
424 HostSubviewLocal Y,
425 int blockIndex,
426 Teuchos::ETransp mode,
427 const LSC alpha,
428 const LSC beta) const
429{
430 #ifdef HAVE_IFPACK2_DEBUG
431 TEUCHOS_TEST_FOR_EXCEPTION(
432 X.extent (0) != Y.extent (0),
433 std::logic_error, "Ifpack2::BandedContainer::solveBlock: X and Y have "
434 "incompatible dimensions (" << X.extent (0) << " resp. "
435 << Y.extent (0) << "). Please report this bug to "
436 "the Ifpack2 developers.");
437 TEUCHOS_TEST_FOR_EXCEPTION(
438 X.extent (0) != static_cast<size_t> (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numCols() : diagBlocks_[blockIndex].numRows()),
439 std::logic_error, "Ifpack2::BandedContainer::solveBlock: The input "
440 "multivector X has incompatible dimensions from those of the "
441 "inverse operator (" << X.extent (0) << " vs. "
442 << (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numCols() : diagBlocks_[blockIndex].numRows())
443 << "). Please report this bug to the Ifpack2 developers.");
444 TEUCHOS_TEST_FOR_EXCEPTION(
445 Y.extent (0) != static_cast<size_t> (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numRows() : diagBlocks_[blockIndex].numCols()),
446 std::logic_error, "Ifpack2::BandedContainer::solveBlock: The output "
447 "multivector Y has incompatible dimensions from those of the "
448 "inverse operator (" << Y.extent (0) << " vs. "
449 << (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numRows() : diagBlocks_[blockIndex].numCols())
450 << "). Please report this bug to the Ifpack2 developers.");
451 #endif
452
453 size_t numRows = X.extent (0);
454 size_t numVecs = X.extent (1);
455
456 LSC zero = Teuchos::ScalarTraits<LSC>::zero ();
457 if (alpha == zero) { // don't need to solve the linear system
458 if (beta == zero) {
459 // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
460 // any Inf or NaN values in Y (rather than multiplying them by
461 // zero, resulting in NaN values).
462 for(size_t j = 0; j < numVecs; j++)
463 for(size_t i = 0; i < numRows; i++)
464 Y(i, j) = zero;
465 }
466 else { // beta != 0
467 for(size_t j = 0; j < numVecs; j++)
468 for(size_t i = 0; i < numRows; i++)
469 Y(i, j) = beta * Y(i, j);
470 }
471 }
472 else { // alpha != 0; must solve the linear system
473 Teuchos::LAPACK<int, LSC> lapack;
474 // If beta is nonzero or Y is not constant stride, we have to use
475 // a temporary output multivector. It gets a copy of X, since
476 // GBTRS overwrites its (multi)vector input with its output.
477 std::vector<LSC> yTemp(numVecs * numRows);
478 for(size_t j = 0; j < numVecs; j++)
479 {
480 for(size_t i = 0; i < numRows; i++)
481 yTemp[j * numRows + i] = X(i, j);
482 }
483
484 int INFO = 0;
485 const char trans = (mode == Teuchos::CONJ_TRANS ? 'C' : (mode == Teuchos::TRANS ? 'T' : 'N'));
486
487 const int* blockIpiv = &ipiv_[this->blockOffsets_[blockIndex] * this->scalarsPerRow_];
488
489 lapack.GBTRS(trans,
490 diagBlocks_[blockIndex].numCols(),
491 diagBlocks_[blockIndex].lowerBandwidth(),
492 /* enter the real number of superdiagonals (see Teuchos_SerialBandDenseSolver)*/
493 diagBlocks_[blockIndex].upperBandwidth() - diagBlocks_[blockIndex].lowerBandwidth(),
494 numVecs,
495 diagBlocks_[blockIndex].values(),
496 diagBlocks_[blockIndex].stride(),
497 blockIpiv,
498 yTemp.data(),
499 numRows,
500 &INFO);
501
502 if (beta != zero) {
503 for(size_t j = 0; j < numVecs; j++)
504 {
505 for(size_t i = 0; i < numRows; i++)
506 {
507 Y(i, j) *= ISC(beta);
508 Y(i, j) += ISC(alpha * yTemp[j * numRows + i]);
509 }
510 }
511 }
512 else {
513 for(size_t j = 0; j < numVecs; j++)
514 {
515 for(size_t i = 0; i < numRows; i++)
516 Y(i, j) = ISC(yTemp[j * numRows + i]);
517 }
518 }
519
520 TEUCHOS_TEST_FOR_EXCEPTION(
521 INFO != 0, std::runtime_error, "Ifpack2::BandedContainer::solveBlock: "
522 "LAPACK's _GBTRS (solve using LU factorization with partial pivoting) "
523 "failed with INFO = " << INFO << " != 0.");
524 }
525}
526
527template<class MatrixType, class LocalScalarType>
528std::ostream&
530print (std::ostream& os) const
531{
532 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
533 fos.setOutputToRootOnly (0);
534 describe (fos);
535 return os;
536}
537
538template<class MatrixType, class LocalScalarType>
539std::string
541description () const
542{
543 std::ostringstream oss;
544 oss << Teuchos::Describable::description();
545 if (this->isInitialized()) {
546 if (this->isComputed()) {
547 oss << "{status = initialized, computed";
548 }
549 else {
550 oss << "{status = initialized, not computed";
551 }
552 }
553 else {
554 oss << "{status = not initialized, not computed";
555 }
556 oss << "}";
557 return oss.str();
558}
559
560template<class MatrixType, class LocalScalarType>
561void
563describe (Teuchos::FancyOStream& os,
564 const Teuchos::EVerbosityLevel verbLevel) const
565{
566 if(verbLevel==Teuchos::VERB_NONE) return;
567 os << "================================================================================" << std::endl;
568 os << "Ifpack2::BandedContainer" << std::endl;
569 for(int i = 0; i < this->numBlocks_; i++)
570 {
571 os << "Block " << i << ": Number of rows = " << this->blockSizes_[i] << std::endl;
572 os << "Block " << i << ": Number of subdiagonals = " << diagBlocks_[i].lowerBandwidth() << std::endl;
573 os << "Block " << i << ": Number of superdiagonals = " << diagBlocks_[i].upperBandwidth() << std::endl;
574 }
575 os << "isInitialized() = " << this->IsInitialized_ << std::endl;
576 os << "isComputed() = " << this->IsComputed_ << std::endl;
577 os << "================================================================================" << std::endl;
578 os << std::endl;
579}
580
581template<class MatrixType, class LocalScalarType>
583{
584 return "Banded";
585}
586
587} // namespace Ifpack2
588
589#define IFPACK2_BANDEDCONTAINER_INSTANT(S,LO,GO,N) \
590 template class Ifpack2::BandedContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >;
591
592#endif // IFPACK2_BANDEDCONTAINER_HPP
Store and solve a local Banded linear problem.
Definition Ifpack2_BandedContainer_decl.hpp:104
BandedContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< LO > > &partitions, const Teuchos::RCP< const import_type > &, bool pointIndexed)
Constructor.
Definition Ifpack2_BandedContainer_def.hpp:62
virtual void compute() override
Initialize and compute each block.
Definition Ifpack2_BandedContainer_def.hpp:242
virtual void initialize() override
Do all set-up operations that only require matrix structure.
Definition Ifpack2_BandedContainer_def.hpp:208
virtual std::ostream & print(std::ostream &os) const override
Print information about this object to the given output stream.
Definition Ifpack2_BandedContainer_def.hpp:530
virtual ~BandedContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition Ifpack2_BandedContainer_def.hpp:79
virtual void setParameters(const Teuchos::ParameterList &List) override
Set all necessary parameters (super- and sub-diagonal bandwidth)
Definition Ifpack2_BandedContainer_def.hpp:83
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition Ifpack2_BandedContainer_def.hpp:582
virtual std::string description() const override
A one-line description of this object.
Definition Ifpack2_BandedContainer_def.hpp:541
void computeBandwidth()
Definition Ifpack2_BandedContainer_def.hpp:101
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const override
Print the object with some verbosity level to the given FancyOStream.
Definition Ifpack2_BandedContainer_def.hpp:563
The implementation of the numerical features of Container (Jacobi, Gauss-Seidel, SGS)....
Definition Ifpack2_Container_decl.hpp:344
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74