Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_TriDiContainer_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_TRIDICONTAINER_DEF_HPP
44 #define IFPACK2_TRIDICONTAINER_DEF_HPP
45 
46 #include "Teuchos_LAPACK.hpp"
47 
48 #ifdef HAVE_MPI
49 # include <mpi.h>
51 #else
52 # include "Teuchos_DefaultSerialComm.hpp"
53 #endif // HAVE_MPI
54 
55 
56 namespace Ifpack2 {
57 
58 template<class MatrixType, class LocalScalarType>
61  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
63  bool pointIndexed) :
64  ContainerImpl<MatrixType, LocalScalarType>(matrix, partitions, pointIndexed),
65  ipiv_(this->blockRows_.size() * this->scalarsPerRow_),
66  scalarOffsets_ (this->numBlocks_)
67 {
69  ! matrix->hasColMap (), std::invalid_argument, "Ifpack2::TriDiContainer: "
70  "The constructor's input matrix must have a column Map.");
71  // FIXME (mfh 25 Aug 2013) What if the matrix's row Map has a
72  // different index base than zero?
73  //compute scalar array offsets (probably different from blockOffsets_)
74  LO scalarTotal = 0;
75  for(LO i = 0; i < this->numBlocks_; i++)
76  {
77  scalarOffsets_[i] = scalarTotal;
78  //actualBlockRows is how many scalars it takes to store the diagonal for block i
79  LO actualBlockRows = this->scalarsPerRow_ * this->blockSizes_[i];
80  if(actualBlockRows == 1)
81  {
82  scalarTotal += 1;
83  }
84  else
85  {
86  //this formula is exact for any block size of 1 or more
87  //includes 1 subdiagonal and 2 superdiagonals.
88  scalarTotal += 4 * (actualBlockRows - 1);
89  }
90  }
91  //Allocate scalar arrays
92  scalars_.resize(scalarTotal);
93  diagBlocks_.reserve(this->numBlocks_);
94  for(int i = 0; i < this->numBlocks_; i++)
95  {
96  diagBlocks_.emplace_back(Teuchos::View, scalars_.data() + scalarOffsets_[i], this->blockSizes_[i] * this->scalarsPerRow_);
97  }
98 }
99 
100 template<class MatrixType, class LocalScalarType>
103 
104 template<class MatrixType, class LocalScalarType>
106 {
107  for(int i = 0; i < this->numBlocks_; i++)
108  diagBlocks_[i].putScalar(Teuchos::ScalarTraits<LSC>::zero());
109  this->IsInitialized_ = true;
110  // We assume that if you called this method, you intend to recompute
111  // everything.
112  this->IsComputed_ = false;
113 }
114 
115 template<class MatrixType, class LocalScalarType>
117 {
118  #ifdef HAVE_IFPACK2_DEBUG
120  ipiv_.size () != this->blockRows_.size( * scalarsPerRow_), std::logic_error,
121  "Ifpack2::TriDiContainer::compute: ipiv_ array has the wrong size. "
122  "Please report this bug to the Ifpack2 developers.");
123  #endif
124 
125  this->IsComputed_ = false;
126  if (! this->isInitialized ()) {
127  this->initialize ();
128  }
129  extract (); // extract the submatrices
130  factor (); // factor them
131  this->IsComputed_ = true;
132 }
133 
134 template<class MatrixType, class LocalScalarType>
136 {
137  using Teuchos::Array;
138  using Teuchos::ArrayView;
140  const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
141  //To extract diagonal blocks, need to translate local rows to local columns.
142  //Strategy: make a lookup table that translates local cols in the matrix to offsets in blockRows_:
143  //blockOffsets_[b] <= offset < blockOffsets_[b+1]: tests whether the column is in block b.
144  //offset - blockOffsets_[b]: gives the column within block b.
145  //
146  //This provides the block and col within a block in O(1).
147  if(this->scalarsPerRow_ > 1)
148  {
149  Array<LO> colToBlockOffset(this->inputBlockMatrix_->getNodeNumCols(), INVALID);
150  for(int i = 0; i < this->numBlocks_; i++)
151  {
152  //Get the interval where block i is defined in blockRows_
153  LO blockStart = this->blockOffsets_[i];
154  LO blockEnd = blockStart + this->blockSizes_[i];
155  ArrayView<const LO> blockRows = this->getBlockRows(i);
156  //Set the lookup table entries for the columns appearing in block i.
157  //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
158  //this is OK. The values updated here are only needed to process block i's entries.
159  for(size_t j = 0; j < size_t(blockRows.size()); j++)
160  {
161  LO localCol = this->translateRowToCol(blockRows[j]);
162  colToBlockOffset[localCol] = blockStart + j;
163  }
164  for(LO blockRow = 0; blockRow < LO(blockRows.size()); blockRow++)
165  {
166  //get a raw view of the whole block row
167  const LO* indices;
168  SC* values;
169  LO numEntries;
170  LO inputRow = this->blockRows_[blockStart + blockRow];
171  this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values, numEntries);
172  for(LO k = 0; k < numEntries; k++)
173  {
174  LO colOffset = colToBlockOffset[indices[k]];
175  if(blockStart <= colOffset && colOffset < blockEnd)
176  {
177  //This entry does appear in the diagonal block.
178  //(br, bc) identifies the scalar's position in the BlockCrs block.
179  //Convert this to (r, c) which is its position in the container block.
180  LO blockCol = colOffset - blockStart;
181  for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
182  {
183  for(LO br = 0; br < this->bcrsBlockSize_; br++)
184  {
185  LO r = this->bcrsBlockSize_ * blockRow + br;
186  LO c = this->bcrsBlockSize_ * blockCol + bc;
187  auto val = values[k * (this->bcrsBlockSize_ * this->bcrsBlockSize_) + (br + this->bcrsBlockSize_ * bc)];
188  if(val != zero)
189  diagBlocks_[i](r, c) = val;
190  }
191  }
192  }
193  }
194  }
195  }
196  }
197  else
198  {
199  //get the mapping from point-indexed matrix columns to offsets in blockRows_
200  //(this includes regular CrsMatrix columns, in which case bcrsBlockSize_ == 1)
201  Array<LO> colToBlockOffset(this->inputMatrix_->getNodeNumCols() * this->bcrsBlockSize_, INVALID);
202  for(int i = 0; i < this->numBlocks_; i++)
203  {
204  //Get the interval where block i is defined in blockRows_
205  LO blockStart = this->blockOffsets_[i];
206  LO blockEnd = blockStart + this->blockSizes_[i];
207  ArrayView<const LO> blockRows = this->getBlockRows(i);
208  //Set the lookup table entries for the columns appearing in block i.
209  //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
210  //this is OK. The values updated here are only needed to process block i's entries.
211  for(size_t j = 0; j < size_t(blockRows.size()); j++)
212  {
213  //translateRowToCol will return the corresponding split column
214  LO localCol = this->translateRowToCol(blockRows[j]);
215  colToBlockOffset[localCol] = blockStart + j;
216  }
217  for(size_t blockRow = 0; blockRow < size_t(blockRows.size()); blockRow++)
218  {
219  //get a view of the split row
220  LO inputPointRow = this->blockRows_[blockStart + blockRow];
221  auto rowView = this->getInputRowView(inputPointRow);
222  for(size_t k = 0; k < rowView.size(); k++)
223  {
224  LO colOffset = colToBlockOffset[rowView.ind(k)];
225  if(blockStart <= colOffset && colOffset < blockEnd)
226  {
227  LO blockCol = colOffset - blockStart;
228  auto val = rowView.val(k);
229  if(val != zero)
230  diagBlocks_[i](blockRow, blockCol) = rowView.val(k);
231  }
232  }
233  }
234  }
235  }
236 }
237 
238 template<class MatrixType, class LocalScalarType>
239 void TriDiContainer<MatrixType, LocalScalarType>::clearBlocks ()
240 {
241  diagBlocks_.clear();
242  scalars_.clear();
243  ContainerImpl<MatrixType, LocalScalarType>::clearBlocks();
244 }
245 
246 template<class MatrixType, class LocalScalarType>
247 void TriDiContainer<MatrixType, LocalScalarType>::factor ()
248 {
249  for(int i = 0; i < this->numBlocks_; i++)
250  {
252  int INFO = 0;
253  int* blockIpiv = &ipiv_[this->blockOffsets_[i] * this->scalarsPerRow_];
254  lapack.GTTRF (diagBlocks_[i].numRowsCols (),
255  diagBlocks_[i].DL(),
256  diagBlocks_[i].D(),
257  diagBlocks_[i].DU(),
258  diagBlocks_[i].DU2(),
259  blockIpiv, &INFO);
260  // INFO < 0 is a bug.
262  INFO < 0, std::logic_error, "Ifpack2::TriDiContainer::factor: "
263  "LAPACK's _GTTRF (LU factorization with partial pivoting) was called "
264  "incorrectly. INFO = " << INFO << " < 0. "
265  "Please report this bug to the Ifpack2 developers.");
266  // INFO > 0 means the matrix is singular. This is probably an issue
267  // either with the choice of rows the rows we extracted, or with the
268  // input matrix itself.
270  INFO > 0, std::runtime_error, "Ifpack2::TriDiContainer::factor: "
271  "LAPACK's _GTTRF (LU factorization with partial pivoting) reports that the "
272  "computed U factor is exactly singular. U(" << INFO << "," << INFO << ") "
273  "(one-based index i) is exactly zero. This probably means that the input "
274  "matrix has a singular diagonal block.");
275  }
276 }
277 
278 template<class MatrixType, class LocalScalarType>
280 solveBlock(HostSubviewLocal X,
281  HostSubviewLocal Y,
282  int blockIndex,
283  Teuchos::ETransp mode,
284  LSC alpha,
285  LSC beta) const
286 {
288  size_t numVecs = X.extent(1);
289  size_t numRows = X.extent(0);
290 
291  #ifdef HAVE_IFPACK2_DEBUG
293  X.extent (0) != Y.extent (0),
294  std::logic_error, "Ifpack2::TriDiContainer::solveBlock: X and Y have "
295  "incompatible dimensions (" << X.extent (0) << " resp. "
296  << Y.extent (0) << "). Please report this bug to "
297  "the Ifpack2 developers.");
299  X.extent (0) != static_cast<size_t> (diagBlocks_[blockIndex].numRowsCols()),
300  std::logic_error, "Ifpack2::TriDiContainer::solveBlock: The input "
301  "multivector X has incompatible dimensions from those of the "
302  "inverse operator (" << X.extent (0) << " vs. "
303  << (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numRowsCols () : diagBlocks_[blockIndex].numRowsCols())
304  << "). Please report this bug to the Ifpack2 developers.");
305  TEUCHOS_TEST_FOR_EXCEPTION(
306  Y.extent (0) != static_cast<size_t> (diagBlocks_[blockIndex].numRowsCols()),
307  std::logic_error, "Ifpack2::TriDiContainer::solveBlock: The output "
308  "multivector Y has incompatible dimensions from those of the "
309  "inverse operator (" << Y.extent (0) << " vs. "
310  << (mode == Teuchos::NO_TRANS ? diagBlocks_[blockIndex].numRowsCols() : diagBlocks_[blockIndex].numRowsCols ())
311  << "). Please report this bug to the Ifpack2 developers.");
312  #endif
313 
314  if(alpha == zero) { // don't need to solve the linear system
315  if(beta == zero) {
316  // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
317  // any Inf or NaN values in Y (rather than multiplying them by
318  // zero, resulting in NaN values).
319  for(size_t j = 0; j < numVecs; j++)
320  for(size_t i = 0; i < numRows; i++)
321  Y(i, j) = zero;
322  }
323  else { // beta != 0
324  for(size_t j = 0; j < numVecs; j++)
325  for(size_t i = 0; i < numRows; i++)
326  Y(i, j) *= ISC(beta);
327  }
328  }
329  else { // alpha != 0; must solve the linear system
331  // If beta is nonzero or Y is not constant stride, we have to use
332  // a temporary output multivector. It gets a copy of X, since
333  // GETRS overwrites its (multi)vector input with its output.
334 
335  std::vector<LSC> yTemp(numVecs * numRows);
336  for(size_t j = 0; j < numVecs; j++)
337  {
338  for(size_t i = 0; i < numRows; i++)
339  yTemp[j * numRows + i] = X(i, j);
340  }
341 
342  int INFO = 0;
343  const char trans =
344  (mode == Teuchos::CONJ_TRANS ? 'C' : (mode == Teuchos::TRANS ? 'T' : 'N'));
345  int* blockIpiv = &ipiv_[this->blockOffsets_[blockIndex] * this->scalarsPerRow_];
346  lapack.GTTRS (trans,
347  diagBlocks_[blockIndex].numRowsCols(),
348  numVecs,
349  diagBlocks_[blockIndex].DL(),
350  diagBlocks_[blockIndex].D(),
351  diagBlocks_[blockIndex].DU(),
352  diagBlocks_[blockIndex].DU2(),
353  blockIpiv,
354  yTemp.data(),
355  numRows,
356  &INFO);
357 
358  if (beta != zero) {
359  for(size_t j = 0; j < numVecs; j++)
360  {
361  for(size_t i = 0; i < numRows; i++)
362  {
363  Y(i, j) *= ISC(beta);
364  Y(i, j) += ISC(alpha * yTemp[j * numRows + i]);
365  }
366  }
367  }
368  else {
369  for(size_t j = 0; j < numVecs; j++)
370  {
371  for(size_t i = 0; i < numRows; i++)
372  Y(i, j) = ISC(yTemp[j * numRows + i]);
373  }
374  }
375 
376  TEUCHOS_TEST_FOR_EXCEPTION(
377  INFO != 0, std::runtime_error, "Ifpack2::TriDiContainer::solveBlock: "
378  "LAPACK's _GETRS (solve using LU factorization with partial pivoting) "
379  "failed with INFO = " << INFO << " != 0.");
380  }
381 }
382 
383 template<class MatrixType, class LocalScalarType>
384 std::ostream& TriDiContainer<MatrixType, LocalScalarType>::print(std::ostream& os) const
385 {
386  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
387  fos.setOutputToRootOnly(0);
388  describe(fos);
389  return(os);
390 }
391 
392 template<class MatrixType, class LocalScalarType>
394 {
395  std::ostringstream oss;
397  if (this->isInitialized()) {
398  if (this->isComputed()) {
399  oss << "{status = initialized, computed";
400  }
401  else {
402  oss << "{status = initialized, not computed";
403  }
404  }
405  else {
406  oss << "{status = not initialized, not computed";
407  }
408 
409  oss << "}";
410  return oss.str();
411 }
412 
413 template<class MatrixType, class LocalScalarType>
414 void
417  const Teuchos::EVerbosityLevel verbLevel) const
418 {
419  using std::endl;
420  if(verbLevel==Teuchos::VERB_NONE) return;
421  os << "================================================================================" << endl;
422  os << "Ifpack2::TriDiContainer" << endl;
423  os << "Number of blocks = " << this->numBlocks_ << endl;
424  os << "isInitialized() = " << this->IsInitialized_ << endl;
425  os << "isComputed() = " << this->IsComputed_ << endl;
426  os << "================================================================================" << endl;
427  os << endl;
428 }
429 
430 template<class MatrixType, class LocalScalarType>
432 {
433  return "TriDi";
434 }
435 
436 #define IFPACK2_TRIDICONTAINER_INSTANT(S,LO,GO,N) \
437  template class Ifpack2::TriDiContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >;
438 
439 } // namespace Ifpack2
440 
441 #endif
Store and solve a local TriDi linear problem.
Definition: Ifpack2_TriDiContainer_decl.hpp:105
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_TriDiContainer_def.hpp:393
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_TriDiContainer_def.hpp:431
Teuchos::Array< LO > blockSizes_
Number of rows in each block.
Definition: Ifpack2_Container_decl.hpp:295
virtual void compute()
Extract the local diagonal block and prepare the solver.
Definition: Ifpack2_TriDiContainer_def.hpp:116
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container_decl.hpp:291
void GTTRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *dl, const ScalarType *d, const ScalarType *du, const ScalarType *du2, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
void GTTRF(const OrdinalType &n, ScalarType *dl, ScalarType *d, ScalarType *du, ScalarType *du2, OrdinalType *IPIV, OrdinalType *info) const
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)
LO scalarsPerRow_
Definition: Ifpack2_Container_decl.hpp:316
virtual ~TriDiContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_TriDiContainer_def.hpp:102
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_TriDiContainer_def.hpp:416
typename Kokkos::Details::ArithTraits< SC >::val_type ISC
Internal representation of Scalar in Kokkos::View.
Definition: Ifpack2_Container_decl.hpp:135
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_TriDiContainer_def.hpp:384
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual std::string description() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void resize(size_type new_size, const value_type &x=value_type())
TriDiContainer(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_TriDiContainer_def.hpp:60
void solveBlock(HostSubviewLocal X, HostSubviewLocal Y, int blockIndex, Teuchos::ETransp mode, LSC alpha, LSC beta) const
Definition: Ifpack2_TriDiContainer_def.hpp:280
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_TriDiContainer_def.hpp:105