43 #ifndef IFPACK2_DENSECONTAINER_DEF_HPP
44 #define IFPACK2_DENSECONTAINER_DEF_HPP
46 #include "Tpetra_CrsMatrix.hpp"
48 #include "Tpetra_BlockMultiVector.hpp"
54 # include "Teuchos_DefaultSerialComm.hpp"
59 template<
class MatrixType,
class LocalScalarType>
65 ContainerImpl<MatrixType, LocalScalarType> (matrix, partitions, pointIndexed),
66 scalarOffsets_ (this->numBlocks_)
69 !matrix->hasColMap(), std::invalid_argument,
"Ifpack2::DenseContainer: "
70 "The constructor's input matrix must have a column Map.");
76 scalarOffsets_[i] = totalScalars;
80 scalars_.
resize(totalScalars);
85 diagBlocks_.emplace_back(
Teuchos::View, scalars_.
data() + scalarOffsets_[i], denseRows, denseRows, denseRows);
91 template<
class MatrixType,
class LocalScalarType>
95 template<
class MatrixType,
class LocalScalarType>
101 for(
int i = 0; i < this->numBlocks_; i++)
103 std::fill (ipiv_.begin (), ipiv_.end (), 0);
105 this->IsInitialized_ =
true;
108 this->IsComputed_ =
false;
111 template<
class MatrixType,
class LocalScalarType>
122 this->IsComputed_ =
false;
123 if (!this->isInitialized ()) {
129 this->IsComputed_ =
true;
132 template<
class MatrixType,
class LocalScalarType>
138 SC zero = STS::zero();
146 if(this->scalarsPerRow_ > 1)
148 Array<LO> colToBlockOffset(this->inputBlockMatrix_->getNodeNumCols(), INVALID);
149 for(
int i = 0; i < this->numBlocks_; i++)
152 LO blockStart = this->blockOffsets_[i];
153 LO blockEnd = blockStart + this->blockSizes_[i];
154 ArrayView<const LO> blockRows = this->getBlockRows(i);
158 for(
size_t j = 0; j < size_t(blockRows.size()); j++)
160 LO localCol = this->translateRowToCol(blockRows[j]);
161 colToBlockOffset[localCol] = blockStart + j;
163 for(LO blockRow = 0; blockRow < LO(blockRows.size()); blockRow++)
169 LO inputRow = this->blockRows_[blockStart + blockRow];
170 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values, numEntries);
171 for(LO k = 0; k < numEntries; k++)
173 LO colOffset = colToBlockOffset[indices[k]];
174 if(blockStart <= colOffset && colOffset < blockEnd)
179 LO blockCol = colOffset - blockStart;
180 for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
182 for(LO br = 0; br < this->bcrsBlockSize_; br++)
184 LO r = this->bcrsBlockSize_ * blockRow + br;
185 LO c = this->bcrsBlockSize_ * blockCol + bc;
186 auto val = values[k * (this->bcrsBlockSize_ * this->bcrsBlockSize_) + (br + this->bcrsBlockSize_ * bc)];
188 diagBlocks_[i](r, c) = val;
200 Array<LO> colToBlockOffset(this->inputMatrix_->getNodeNumCols() * this->bcrsBlockSize_, INVALID);
201 for(
int i = 0; i < this->numBlocks_; i++)
204 LO blockStart = this->blockOffsets_[i];
205 LO blockEnd = blockStart + this->blockSizes_[i];
206 ArrayView<const LO> blockRows = this->getBlockRows(i);
210 for(
size_t j = 0; j < size_t(blockRows.size()); j++)
213 LO localCol = this->translateRowToCol(blockRows[j]);
214 colToBlockOffset[localCol] = blockStart + j;
216 for(
size_t blockRow = 0; blockRow < size_t(blockRows.size()); blockRow++)
219 LO inputPointRow = this->blockRows_[blockStart + blockRow];
220 auto rowView = this->getInputRowView(inputPointRow);
221 for(
size_t k = 0; k < rowView.size(); k++)
223 LO colOffset = colToBlockOffset[rowView.ind(k)];
224 if(blockStart <= colOffset && colOffset < blockEnd)
226 LO blockCol = colOffset - blockStart;
227 auto val = rowView.val(k);
229 diagBlocks_[i](blockRow, blockCol) = rowView.val(k);
238 template<
class MatrixType,
class LocalScalarType>
240 DenseContainer<MatrixType, LocalScalarType>::
244 for(
int i = 0; i < this->numBlocks_; i++)
247 int* blockIpiv = &ipiv_[this->blockOffsets_[i] * this->scalarsPerRow_];
248 lapack.
GETRF(diagBlocks_[i].numRows(),
249 diagBlocks_[i].numCols(),
250 diagBlocks_[i].values(),
251 diagBlocks_[i].stride(),
255 INFO < 0, std::logic_error,
"Ifpack2::DenseContainer::factor: "
256 "LAPACK's _GETRF (LU factorization with partial pivoting) was called "
257 "incorrectly. INFO = " << INFO <<
" < 0. "
258 "Please report this bug to the Ifpack2 developers.");
263 INFO > 0, std::runtime_error,
"Ifpack2::DenseContainer::factor: "
264 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the "
265 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") "
266 "(one-based index i) is exactly zero. This probably means that the input "
267 "matrix has a singular diagonal block.");
271 template<
class MatrixType,
class LocalScalarType>
273 DenseContainer<MatrixType, LocalScalarType>::
274 solveBlock(HostSubviewLocal X,
281 #ifdef HAVE_IFPACK2_DEBUG
283 X.extent (0) != Y.extent (0),
284 std::logic_error,
"Ifpack2::DenseContainer::solveBlock: X and Y have "
285 "incompatible dimensions (" << X.extent (0) <<
" resp. "
286 << Y.extent (0) <<
"). Please report this bug to "
287 "the Ifpack2 developers.");
290 X.extent (1) != Y.extent(1),
291 std::logic_error,
"Ifpack2::DenseContainer::solveBlock: X and Y have "
292 "incompatible numbers of vectors (" << X.extent (1) <<
" resp. "
293 << Y.extent (1) <<
"). Please report this bug to "
294 "the Ifpack2 developers.");
298 size_t numRows = X.extent(0);
299 size_t numVecs = X.extent(1);
300 if(alpha == STS::zero()) {
301 if(beta == STS::zero()) {
305 for(
size_t j = 0; j < numVecs; j++)
307 for(
size_t i = 0; i < numRows; i++)
308 Y(i, j) = STS::zero();
312 for(
size_t j = 0; j < numVecs; j++)
314 for(
size_t i = 0; i < numRows; i++)
315 Y(i, j) = beta * Y(i, j);
323 std::vector<LSC> yTemp(numVecs * numRows);
324 for(
size_t j = 0; j < numVecs; j++)
326 for(
size_t i = 0; i < numRows; i++)
327 yTemp[j * numRows + i] = X(i, j);
331 int* blockIpiv = &ipiv_[this->blockOffsets_[blockIndex] * this->scalarsPerRow_];
335 diagBlocks_[blockIndex].numRows (),
337 diagBlocks_[blockIndex].values (),
338 diagBlocks_[blockIndex].stride (),
344 if (beta != STS::zero ()) {
345 for(
size_t j = 0; j < numVecs; j++)
347 for(
size_t i = 0; i < numRows; i++)
349 Y(i, j) *= ISC(beta);
350 Y(i, j) += ISC(alpha * yTemp[j * numRows + i]);
355 for(
size_t j = 0; j < numVecs; j++)
357 for(
size_t i = 0; i < numRows; i++)
358 Y(i, j) = ISC(yTemp[j * numRows + i]);
363 INFO != 0, std::runtime_error,
"Ifpack2::DenseContainer::solveBlock: "
364 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) "
365 "failed with INFO = " << INFO <<
" != 0.");
369 template<
class MatrixType,
class LocalScalarType>
376 this->describe (fos);
380 template<
class MatrixType,
class LocalScalarType>
385 std::ostringstream oss;
386 oss <<
"Ifpack::DenseContainer: ";
387 if (this->isInitialized()) {
388 if (this->isComputed()) {
389 oss <<
"{status = initialized, computed";
392 oss <<
"{status = initialized, not computed";
396 oss <<
"{status = not initialized, not computed";
403 template<
class MatrixType,
class LocalScalarType>
411 os <<
"================================================================================" << endl;
412 os <<
"Ifpack2::DenseContainer" << endl;
413 for(
int i = 0; i < this->numBlocks_; i++)
415 os <<
"Block " << i <<
" number of rows = " << this->blockSizes_[i] << endl;
417 os <<
"isInitialized() = " << this->IsInitialized_ << endl;
418 os <<
"isComputed() = " << this->IsComputed_ << endl;
419 os <<
"================================================================================" << endl;
423 template<
class MatrixType,
class LocalScalarType>
431 template<
class MatrixType,
class LocalScalarType>
443 #define IFPACK2_DENSECONTAINER_INSTANT(S,LO,GO,N) \
444 template class Ifpack2::DenseContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >;
446 #endif // IFPACK2_DENSECONTAINER_DEF_HPP
Teuchos::Array< LO > blockRows_
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container_decl.hpp:293
Teuchos::Array< LO > blockSizes_
Number of rows in each block.
Definition: Ifpack2_Container_decl.hpp:295
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container_decl.hpp:291
The implementation of the numerical features of Container (Jacobi, Gauss-Seidel, SGS). This class allows a custom scalar type (LocalScalarType) to be used for storing blocks and solving the block systems. Hiding this template parameter from the Container interface simplifies the BlockRelaxation and ContainerFactory classes.
Definition: Ifpack2_Container_decl.hpp:342
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void compute()
Extract the local diagonal block and prepare the solver.
Definition: Ifpack2_DenseContainer_def.hpp:114
LO scalarsPerRow_
Definition: Ifpack2_Container_decl.hpp:316
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_DenseContainer_def.hpp:432
DenseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< LO > > &partitions, const Teuchos::RCP< const import_type > &importer, bool pointIndexed)
Constructor.
Definition: Ifpack2_DenseContainer_def.hpp:61
void GETRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, OrdinalType *info) const
Store and solve a local dense linear problem.
Definition: Ifpack2_DenseContainer_decl.hpp:103
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_DenseContainer_def.hpp:372
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to the given FancyOStream.
Definition: Ifpack2_DenseContainer_def.hpp:406
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void resize(size_type new_size, const value_type &x=value_type())
virtual ~DenseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_DenseContainer_def.hpp:93
void GETRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_DenseContainer_def.hpp:383
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_DenseContainer_def.hpp:98