Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_SparseContainer_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_SPARSECONTAINER_DEF_HPP
44#define IFPACK2_SPARSECONTAINER_DEF_HPP
45
47#ifdef HAVE_MPI
48#include <mpi.h>
49#include "Teuchos_DefaultMpiComm.hpp"
50#else
51#include "Teuchos_DefaultSerialComm.hpp"
52#endif
53#include "Teuchos_TestForException.hpp"
54
55namespace Ifpack2 {
56
57//==============================================================================
58template<class MatrixType, class InverseType>
60SparseContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
61 const Teuchos::Array<Teuchos::Array<LO> >& partitions,
62 const Teuchos::RCP<const import_type>&,
63 bool pointIndexed) :
64 ContainerImpl<MatrixType, InverseScalar> (matrix, partitions, pointIndexed),
65#ifdef HAVE_MPI
66 localComm_ (Teuchos::rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF)))
67#else
68 localComm_ (Teuchos::rcp (new Teuchos::SerialComm<int> ()))
69#endif // HAVE_MPI
70{}
71
72//==============================================================================
73template<class MatrixType, class InverseType>
76
77//==============================================================================
78template<class MatrixType, class InverseType>
79void SparseContainer<MatrixType,InverseType>::setParameters(const Teuchos::ParameterList& List)
80{
81 List_ = List;
82}
83
84//==============================================================================
85template<class MatrixType, class InverseType>
87{
88 // We assume that if you called this method, you intend to recompute
89 // everything. Thus, we release references to all the internal
90 // objects. We do this first to save memory. (In an RCP
91 // assignment, the right-hand side and left-hand side coexist before
92 // the left-hand side's reference count gets updated.)
93
94 //Will create the diagonal blocks and their inverses
95 //in extract()
96 diagBlocks_.assign(this->numBlocks_, Teuchos::null);
97 Inverses_.assign(this->numBlocks_, Teuchos::null);
98
99 this->IsInitialized_ = true;
100 this->IsComputed_ = false;
101}
102
103//==============================================================================
104template<class MatrixType, class InverseType>
106{
107 this->IsComputed_ = false;
108 if (!this->isInitialized ()) {
109 this->initialize ();
110 }
111
112 // Extract the submatrices.
113 this->extract();
114
115 // The inverse operator already has a pointer to the submatrix. Now
116 // that the submatrix is filled in, we can initialize and compute
117 // the inverse operator.
118 for(int i = 0; i < this->numBlocks_; i++)
119 {
120 Inverses_[i]->setParameters(List_);
121 Inverses_[i]->initialize ();
122 }
123 for(int i = 0; i < this->numBlocks_; i++)
124 Inverses_[i]->compute ();
125 this->IsComputed_ = true;
126}
127
128//==============================================================================
129template<class MatrixType, class InverseType>
131{
132 for(auto inv : Inverses_)
133 delete inv.get();
134 Inverses_.clear();
135 diagBlocks_.clear();
137}
138
139//==============================================================================
140template<class MatrixType, class InverseType>
142solveBlockMV(const inverse_mv_type& X,
143 inverse_mv_type& Y,
144 int blockIndex,
145 Teuchos::ETransp mode,
146 InverseScalar alpha,
147 InverseScalar beta) const
148{
149 TEUCHOS_TEST_FOR_EXCEPTION(
150 Inverses_[blockIndex]->getDomainMap()->getLocalNumElements() != X.getLocalLength(),
151 std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
152 "operator and X have incompatible dimensions (" <<
153 Inverses_[blockIndex]->getDomainMap()->getLocalNumElements() << " resp. "
154 << X.getLocalLength() << "). Please report this bug to "
155 "the Ifpack2 developers.");
156 TEUCHOS_TEST_FOR_EXCEPTION(
157 Inverses_[blockIndex]->getRangeMap()->getLocalNumElements() != Y.getLocalLength(),
158 std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
159 "operator and Y have incompatible dimensions (" <<
160 Inverses_[blockIndex]->getRangeMap()->getLocalNumElements() << " resp. "
161 << Y.getLocalLength() << "). Please report this bug to "
162 "the Ifpack2 developers.");
163 Inverses_[blockIndex]->apply(X, Y, mode, alpha, beta);
164}
165
166template<class MatrixType, class InverseType>
168apply (ConstHostView X,
169 HostView Y,
170 int blockIndex,
171 Teuchos::ETransp mode,
172 SC alpha,
173 SC beta) const
174{
175 using Teuchos::ArrayView;
176 using Teuchos::as;
177
178 // The InverseType template parameter might have different template
179 // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
180 // example, MatrixType (a global object) might use a bigger GO
181 // (global ordinal) type than InverseType (which deals with the
182 // diagonal block, a local object). This means that we might have
183 // to convert X and Y to the Tpetra::MultiVector specialization that
184 // InverseType wants. This class' X_ and Y_ internal fields are of
185 // the right type for InverseType, so we can use those as targets.
186
187 // Tpetra::MultiVector specialization corresponding to InverseType.
189 size_t numVecs = X.extent(1);
190
191 TEUCHOS_TEST_FOR_EXCEPTION(
192 ! this->IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply: "
193 "You must have called the compute() method before you may call apply(). "
194 "You may call the apply() method as many times as you want after calling "
195 "compute() once, but you must have called compute() at least once.");
196 TEUCHOS_TEST_FOR_EXCEPTION(
197 X.extent(1) != Y.extent(1), std::runtime_error,
198 "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
199 "vectors. X has " << X.extent(1)
200 << ", but Y has " << Y.extent(1) << ".");
201
202 if (numVecs == 0) {
203 return; // done! nothing to do
204 }
205
206 const LO numRows = this->blockSizes_[blockIndex];
207
208 // The operator Inverse_ works on a permuted subset of the local
209 // parts of X and Y. The subset and permutation are defined by the
210 // index array returned by getBlockRows(). If the permutation is
211 // trivial and the subset is exactly equal to the local indices,
212 // then we could use the local parts of X and Y exactly, without
213 // needing to permute. Otherwise, we have to use temporary storage
214 // to permute X and Y. For now, we always use temporary storage.
215 //
216 // Create temporary permuted versions of the input and output.
217 // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
218 // store the permuted versions of X resp. Y. Note that X_local has
219 // the domain Map of the operator, which may be a permuted subset of
220 // the local Map corresponding to X.getMap(). Similarly, Y_local
221 // has the range Map of the operator, which may be a permuted subset
222 // of the local Map corresponding to Y.getMap(). numRows here
223 // gives the number of rows in the row Map of the local Inverse_
224 // operator.
225 //
226 // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
227 // here that the row Map and the range Map of that operator are
228 // the same.
229 //
230 // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
231 // really belongs in Tpetra.
232
233 if(invX.size() == 0)
234 {
235 for(LO i = 0; i < this->numBlocks_; i++)
236 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
237 for(LO i = 0; i < this->numBlocks_; i++)
238 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
239 }
240 inverse_mv_type& X_local = invX[blockIndex];
241 TEUCHOS_TEST_FOR_EXCEPTION(
242 X_local.getLocalLength() != size_t(numRows * this->scalarsPerRow_), std::logic_error,
243 "Ifpack2::SparseContainer::apply: "
244 "X_local has length " << X_local.getLocalLength() << ", which does "
245 "not match numRows = " << numRows * this->scalarsPerRow_ << ". Please report this bug to "
246 "the Ifpack2 developers.");
247 const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
248 if(this->scalarsPerRow_ == 1)
249 mvgs.gatherMVtoView(X_local, X, blockRows);
250 else
251 mvgs.gatherMVtoViewBlock(X_local, X, blockRows, this->scalarsPerRow_);
252
253 // We must gather the output multivector Y even on input to
254 // Inverse_->apply(), since the Inverse_ operator might use it as an
255 // initial guess for a linear solve. We have no way of knowing
256 // whether it does or does not.
257
258 inverse_mv_type& Y_local = invY[blockIndex];
259 TEUCHOS_TEST_FOR_EXCEPTION(
260 Y_local.getLocalLength () != size_t(numRows * this->scalarsPerRow_), std::logic_error,
261 "Ifpack2::SparseContainer::apply: "
262 "Y_local has length " << Y_local.getLocalLength () << ", which does "
263 "not match numRows = " << numRows * this->scalarsPerRow_ << ". Please report this bug to "
264 "the Ifpack2 developers.");
265
266 if(this->scalarsPerRow_ == 1)
267 mvgs.gatherMVtoView(Y_local, Y, blockRows);
268 else
269 mvgs.gatherMVtoViewBlock(Y_local, Y, blockRows, this->scalarsPerRow_);
270
271 // Apply the local operator:
272 // Y_local := beta*Y_local + alpha*M^{-1}*X_local
273 this->solveBlockMV(X_local, Y_local, blockIndex, mode,
274 InverseScalar(alpha), InverseScalar(beta));
275
276
277 // Scatter the permuted subset output vector Y_local back into the
278 // original output multivector Y.
279 if(this->scalarsPerRow_ == 1)
280 mvgs.scatterMVtoView(Y, Y_local, blockRows);
281 else
282 mvgs.scatterMVtoViewBlock(Y, Y_local, blockRows, this->scalarsPerRow_);
283}
284
285//==============================================================================
286template<class MatrixType, class InverseType>
288weightedApply (ConstHostView X,
289 HostView Y,
290 ConstHostView D,
291 int blockIndex,
292 Teuchos::ETransp mode,
293 SC alpha,
294 SC beta) const
295{
296 using Teuchos::ArrayView;
297 using Teuchos::Range1D;
298 using std::cerr;
299 using std::endl;
300 typedef Teuchos::ScalarTraits<InverseScalar> STS;
301
302 // The InverseType template parameter might have different template
303 // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
304 // example, MatrixType (a global object) might use a bigger GO
305 // (global ordinal) type than InverseType (which deals with the
306 // diagonal block, a local object). This means that we might have
307 // to convert X and Y to the Tpetra::MultiVector specialization that
308 // InverseType wants. This class' X_ and Y_ internal fields are of
309 // the right type for InverseType, so we can use those as targets.
310
311 // Tpetra::Vector specialization corresponding to InverseType.
312 typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal, InverseGlobalOrdinal, InverseNode> inverse_vector_type;
313
315 const size_t numVecs = X.extent(1);
316
317 TEUCHOS_TEST_FOR_EXCEPTION(
318 ! this->IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::"
319 "weightedApply: You must have called the compute() method before you may "
320 "call apply(). You may call the apply() method as many times as you want "
321 "after calling compute() once, but you must have called compute() at least "
322 "once.");
323 TEUCHOS_TEST_FOR_EXCEPTION(
324 X.extent(1) != Y.extent(1), std::runtime_error,
325 "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
326 "of vectors. X has " << X.extent(1) << ", but Y has "
327 << Y.extent(1) << ".");
328
329 //bmk 7-2019: BlockRelaxation already checked this, but if that changes...
330 TEUCHOS_TEST_FOR_EXCEPTION(
331 this->scalarsPerRow_ > 1, std::logic_error, "Ifpack2::SparseContainer::weightedApply: "
332 "Use of block rows isn't allowed in overlapping Jacobi (the only method that uses weightedApply");
333
334 if (numVecs == 0) {
335 return; // done! nothing to do
336 }
337
338 // The operator Inverse_ works on a permuted subset of the local
339 // parts of X and Y. The subset and permutation are defined by the
340 // index array returned by getBlockRows(). If the permutation is
341 // trivial and the subset is exactly equal to the local indices,
342 // then we could use the local parts of X and Y exactly, without
343 // needing to permute. Otherwise, we have to use temporary storage
344 // to permute X and Y. For now, we always use temporary storage.
345 //
346 // Create temporary permuted versions of the input and output.
347 // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
348 // store the permuted versions of X resp. Y. Note that X_local has
349 // the domain Map of the operator, which may be a permuted subset of
350 // the local Map corresponding to X.getMap(). Similarly, Y_local
351 // has the range Map of the operator, which may be a permuted subset
352 // of the local Map corresponding to Y.getMap(). numRows here
353 // gives the number of rows in the row Map of the local Inverse_
354 // operator.
355 //
356 // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
357 // here that the row Map and the range Map of that operator are
358 // the same.
359 //
360 // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
361 // really belongs in Tpetra.
362 const LO numRows = this->blockSizes_[blockIndex];
363
364 if(invX.size() == 0)
365 {
366 for(LO i = 0; i < this->numBlocks_; i++)
367 invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
368 for(LO i = 0; i < this->numBlocks_; i++)
369 invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
370 }
371 inverse_mv_type& X_local = invX[blockIndex];
372 const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
373 mvgs.gatherMVtoView(X_local, X, blockRows);
374
375 // We must gather the output multivector Y even on input to
376 // Inverse_->apply(), since the Inverse_ operator might use it as an
377 // initial guess for a linear solve. We have no way of knowing
378 // whether it does or does not.
379
380 inverse_mv_type Y_local = invY[blockIndex];
381 TEUCHOS_TEST_FOR_EXCEPTION(
382 Y_local.getLocalLength() != size_t(numRows), std::logic_error,
383 "Ifpack2::SparseContainer::weightedApply: "
384 "Y_local has length " << X_local.getLocalLength() << ", which does "
385 "not match numRows = " << numRows << ". Please report this bug to "
386 "the Ifpack2 developers.");
387 mvgs.gatherMVtoView(Y_local, Y, blockRows);
388
389 // Apply the diagonal scaling D to the input X. It's our choice
390 // whether the result has the original input Map of X, or the
391 // permuted subset Map of X_local. If the latter, we also need to
392 // gather D into the permuted subset Map. We choose the latter, to
393 // save memory and computation. Thus, we do the following:
394 //
395 // 1. Gather D into a temporary vector D_local.
396 // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
397 // 3. Compute X_scaled := diag(D_loca) * X_local.
398
399 inverse_vector_type D_local(Inverses_[blockIndex]->getDomainMap());
400 TEUCHOS_TEST_FOR_EXCEPTION(
401 D_local.getLocalLength() != size_t(this->blockSizes_[blockIndex]), std::logic_error,
402 "Ifpack2::SparseContainer::weightedApply: "
403 "D_local has length " << X_local.getLocalLength () << ", which does "
404 "not match numRows = " << this->blockSizes_[blockIndex] << ". Please report this bug to "
405 "the Ifpack2 developers.");
406 mvgs.gatherMVtoView(D_local, D, blockRows);
407 inverse_mv_type X_scaled(Inverses_[blockIndex]->getDomainMap(), numVecs);
408 X_scaled.elementWiseMultiply(STS::one(), D_local, X_local, STS::zero());
409
410 // Y_temp will hold the result of M^{-1}*X_scaled. If beta == 0, we
411 // can write the result of Inverse_->apply() directly to Y_local, so
412 // Y_temp may alias Y_local. Otherwise, if beta != 0, we need
413 // temporary storage for M^{-1}*X_scaled, so Y_temp must be
414 // different than Y_local.
415 inverse_mv_type* Y_temp;
416 if (InverseScalar(beta) == STS::zero ()) {
417 Y_temp = &Y_local;
418 } else {
419 Y_temp = new inverse_mv_type(Inverses_[blockIndex]->getRangeMap(), numVecs);
420 }
421 // Apply the local operator: Y_temp := M^{-1} * X_scaled
422 Inverses_[blockIndex]->apply(X_scaled, *Y_temp, mode);
423 // Y_local := beta * Y_local + alpha * diag(D_local) * Y_tmp.
424 //
425 // Note that we still use the permuted subset scaling D_local here,
426 // because Y_temp has the same permuted subset Map. That's good, in
427 // fact, because it's a subset: less data to read and multiply.
428 Y_local.elementWiseMultiply(alpha, D_local, *Y_temp, beta);
429 if(Y_temp != &Y_local)
430 delete Y_temp;
431 // Copy the permuted subset output vector Y_local into the original
432 // output multivector Y.
433 mvgs.scatterMVtoView(Y, Y_local, blockRows);
434}
435
436//==============================================================================
437template<class MatrixType, class InverseType>
438std::ostream& SparseContainer<MatrixType,InverseType>::print(std::ostream& os) const
439{
440 Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
441 fos.setOutputToRootOnly(0);
442 describe(fos);
443 return(os);
444}
445
446//==============================================================================
447template<class MatrixType, class InverseType>
449{
450 std::ostringstream oss;
451 oss << "\"Ifpack2::SparseContainer\": {";
452 if (this->isInitialized()) {
453 if (this->isComputed()) {
454 oss << "status = initialized, computed";
455 }
456 else {
457 oss << "status = initialized, not computed";
458 }
459 }
460 else {
461 oss << "status = not initialized, not computed";
462 }
463 for(int i = 0; i < this->numBlocks_; i++)
464 {
465 oss << ", Block Inverse " << i << ": {";
466 oss << Inverses_[i]->description();
467 oss << "}";
468 }
469 oss << "}";
470 return oss.str();
471}
472
473//==============================================================================
474template<class MatrixType, class InverseType>
475void SparseContainer<MatrixType,InverseType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
476{
477 using std::endl;
478 if(verbLevel==Teuchos::VERB_NONE) return;
479 os << "================================================================================" << endl;
480 os << "Ifpack2::SparseContainer" << endl;
481 for(int i = 0; i < this->numBlocks_; i++)
482 {
483 os << "Block " << i << " rows: = " << this->blockSizes_[i] << endl;
484 }
485 os << "isInitialized() = " << this->IsInitialized_ << endl;
486 os << "isComputed() = " << this->IsComputed_ << endl;
487 os << "================================================================================" << endl;
488 os << endl;
489}
490
491//==============================================================================
492template<class MatrixType, class InverseType>
494extract ()
495{
496 using Teuchos::Array;
497 using Teuchos::ArrayView;
498 using Teuchos::RCP;
499 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
500 //To extract diagonal blocks, need to translate local rows to local columns.
501 //Strategy: make a lookup table that translates local cols in the matrix to offsets in blockRows_:
502 //blockOffsets_[b] <= offset < blockOffsets_[b+1]: tests whether the column is in block b.
503 //offset - blockOffsets_[b]: gives the column within block b.
504 //
505 //This provides the block and col within a block in O(1).
506 Teuchos::Array<InverseGlobalOrdinal> indicesToInsert;
507 Teuchos::Array<InverseScalar> valuesToInsert;
508 if(this->scalarsPerRow_ > 1)
509 {
510 Array<LO> colToBlockOffset(this->inputBlockMatrix_->getLocalNumCols(), INVALID);
511 for(int i = 0; i < this->numBlocks_; i++)
512 {
513 //Get the interval where block i is defined in blockRows_
514 LO blockStart = this->blockOffsets_[i];
515 LO blockSize = this->blockSizes_[i];
516 LO blockPointSize = this->blockSizes_[i] * this->scalarsPerRow_;
517 LO blockEnd = blockStart + blockSize;
518 ArrayView<const LO> blockRows = this->getBlockRows(i);
519 //Set the lookup table entries for the columns appearing in block i.
520 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
521 //this is OK. The values updated here are only needed to process block i's entries.
522 for(LO j = 0; j < blockSize; j++)
523 {
524 LO localCol = this->translateRowToCol(blockRows[j]);
525 colToBlockOffset[localCol] = blockStart + j;
526 }
527 //First, count the number of entries in each row of block i
528 //(to allocate it with StaticProfile)
529 Array<size_t> rowEntryCounts(blockPointSize, 0);
530 //blockRow counts the BlockCrs LIDs that are going into this block
531 //Rows are inserted into the CrsMatrix in sequential order
532 using inds_type = typename block_crs_matrix_type::local_inds_host_view_type;
533 using vals_type = typename block_crs_matrix_type::values_host_view_type;
534 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
535 {
536 //get a raw view of the whole block row
537 inds_type indices;
538 vals_type values;
539 LO inputRow = this->blockRows_[blockStart + blockRow];
540 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
541 LO numEntries = (LO) indices.size();
542 for(LO br = 0; br < this->bcrsBlockSize_; br++)
543 {
544 for(LO k = 0; k < numEntries; k++)
545 {
546 LO colOffset = colToBlockOffset[indices[k]];
547 if(blockStart <= colOffset && colOffset < blockEnd)
548 {
549 rowEntryCounts[blockRow * this->bcrsBlockSize_ + br] += this->bcrsBlockSize_;
550 }
551 }
552 }
553 }
554 //now that row sizes are known, can allocate the diagonal matrix
555 RCP<InverseMap> tempMap(new InverseMap(blockPointSize, 0, this->localComm_));
556 diagBlocks_[i] = rcp(new InverseCrs(tempMap, rowEntryCounts));
557 Inverses_[i] = rcp(new InverseType(diagBlocks_[i]));
558 //insert the actual entries, one row at a time
559 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
560 {
561 //get a raw view of the whole block row
562 inds_type indices;
563 vals_type values;
564 LO inputRow = this->blockRows_[blockStart + blockRow];
565 this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values);
566 LO numEntries = (LO) indices.size();
567 for(LO br = 0; br < this->bcrsBlockSize_; br++)
568 {
569 indicesToInsert.clear();
570 valuesToInsert.clear();
571 for(LO k = 0; k < numEntries; k++)
572 {
573 LO colOffset = colToBlockOffset[indices[k]];
574 if(blockStart <= colOffset && colOffset < blockEnd)
575 {
576 LO blockCol = colOffset - blockStart;
577 //bc iterates over the columns in (block) entry k
578 for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
579 {
580 indicesToInsert.push_back(blockCol * this->bcrsBlockSize_ + bc);
581 valuesToInsert.push_back(values[k * this->bcrsBlockSize_ * this->bcrsBlockSize_ + bc * this->bcrsBlockSize_ + br]);
582 }
583 }
584 }
585 InverseGlobalOrdinal rowToInsert = blockRow * this->bcrsBlockSize_ + br;
586 if(indicesToInsert.size())
587 diagBlocks_[i]->insertGlobalValues(rowToInsert, indicesToInsert(), valuesToInsert());
588 }
589 }
590 diagBlocks_[i]->fillComplete();
591 }
592 }
593 else
594 {
595 //get the mapping from point-indexed matrix columns to offsets in blockRows_
596 //(this includes regular CrsMatrix columns, in which case scalarsPerRow_ == 1)
597 Array<LO> colToBlockOffset(this->inputMatrix_->getLocalNumCols() * this->bcrsBlockSize_, INVALID);
598 for(int i = 0; i < this->numBlocks_; i++)
599 {
600 //Get the interval where block i is defined in blockRows_
601 LO blockStart = this->blockOffsets_[i];
602 LO blockSize = this->blockSizes_[i];
603 LO blockEnd = blockStart + blockSize;
604 ArrayView<const LO> blockRows = this->getBlockRows(i);
605 //Set the lookup table entries for the columns appearing in block i.
606 //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
607 //this is OK. The values updated here are only needed to process block i's entries.
608 for(LO j = 0; j < blockSize; j++)
609 {
610 //translateRowToCol will return the corresponding split column
611 LO localCol = this->translateRowToCol(blockRows[j]);
612 colToBlockOffset[localCol] = blockStart + j;
613 }
614 Teuchos::Array<size_t> rowEntryCounts(blockSize, 0);
615 for(LO j = 0; j < blockSize; j++)
616 {
617 rowEntryCounts[j] = this->getInputRowView(this->blockRows_[blockStart + j]).size();
618 }
619 RCP<InverseMap> tempMap(new InverseMap(blockSize, 0, this->localComm_));
620 diagBlocks_[i] = rcp(new InverseCrs(tempMap, rowEntryCounts));
621 Inverses_[i] = rcp(new InverseType(diagBlocks_[i]));
622 for(LO blockRow = 0; blockRow < blockSize; blockRow++)
623 {
624 valuesToInsert.clear();
625 indicesToInsert.clear();
626 //get a view of the split row
627 LO inputSplitRow = this->blockRows_[blockStart + blockRow];
628 auto rowView = this->getInputRowView(inputSplitRow);
629 for(size_t k = 0; k < rowView.size(); k++)
630 {
631 LO colOffset = colToBlockOffset[rowView.ind(k)];
632 if(blockStart <= colOffset && colOffset < blockEnd)
633 {
634 LO blockCol = colOffset - blockStart;
635 indicesToInsert.push_back(blockCol);
636 valuesToInsert.push_back(rowView.val(k));
637 }
638 }
639 if(indicesToInsert.size())
640 diagBlocks_[i]->insertGlobalValues(blockRow, indicesToInsert(), valuesToInsert());
641 }
642 diagBlocks_[i]->fillComplete ();
643 }
644 }
645}
646
647template<typename MatrixType, typename InverseType>
649{
650 typedef ILUT<Tpetra::RowMatrix<SC, LO, GO, NO> > ILUTInverse;
651#ifdef HAVE_IFPACK2_AMESOS2
653 if(std::is_same<InverseType, ILUTInverse>::value)
654 {
655 return "SparseILUT";
656 }
657 else if(std::is_same<InverseType, AmesosInverse>::value)
658 {
659 return "SparseAmesos";
660 }
661 else
662 {
663 throw std::logic_error("InverseType for SparseContainer must be Ifpack2::ILUT or Details::Amesos2Wrapper");
664 }
665#else
666 // Macros can't have commas in their arguments, so we have to
667 // compute the bool first argument separately.
668 constexpr bool inverseTypeIsILUT = std::is_same<InverseType, ILUTInverse>::value;
669 TEUCHOS_TEST_FOR_EXCEPTION(! inverseTypeIsILUT, std::logic_error,
670 "InverseType for SparseContainer must be Ifpack2::ILUT<ROW>");
671 return "SparseILUT"; //the only supported sparse container specialization if no Amesos2
672#endif
673}
674
675} // namespace Ifpack2
676
677// For ETI
678#include "Ifpack2_ILUT.hpp"
679#ifdef HAVE_IFPACK2_AMESOS2
680#include "Ifpack2_Details_Amesos2Wrapper.hpp"
681#endif
682
683// There's no need to instantiate for CrsMatrix too. All Ifpack2
684// preconditioners can and should do dynamic casts if they need a type
685// more specific than RowMatrix.
686
687#ifdef HAVE_IFPACK2_AMESOS2
688# define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
689 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
690 Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; \
691 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
692 Ifpack2::Details::Amesos2Wrapper<Tpetra::RowMatrix<S,LO,GO,N> > >;
693#else
694# define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
695 template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S,LO,GO,N>, \
696 Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> > >;
697#endif
698#endif // IFPACK2_SPARSECONTAINER_DEF_HPP
Ifpack2::SparseContainer class declaration.
The implementation of the numerical features of Container (Jacobi, Gauss-Seidel, SGS)....
Definition Ifpack2_Container_decl.hpp:344
typename mv_type::dual_view_type::t_host HostView
Definition Ifpack2_Container_decl.hpp:139
Wrapper class for direct solvers in Amesos2.
Definition Ifpack2_Details_Amesos2Wrapper_decl.hpp:111
Implementation detail of Ifpack2::Container subclasses.
Definition Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition Ifpack2_ILUT_decl.hpp:100
Store and solve a local sparse linear problem.
Definition Ifpack2_SparseContainer_decl.hpp:136
virtual void apply(ConstHostView X, HostView Y, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * M^{-1} X + beta*Y.
Definition Ifpack2_SparseContainer_def.hpp:168
SparseContainer(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_SparseContainer_def.hpp:60
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_SparseContainer_def.hpp:475
virtual void weightedApply(ConstHostView X, HostView Y, ConstHostView W, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Definition Ifpack2_SparseContainer_def.hpp:288
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition Ifpack2_SparseContainer_def.hpp:648
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition Ifpack2_SparseContainer_def.hpp:86
void clearBlocks()
Definition Ifpack2_SparseContainer_def.hpp:130
virtual void compute()
Initialize and compute all blocks.
Definition Ifpack2_SparseContainer_def.hpp:105
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition Ifpack2_SparseContainer_def.hpp:75
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition Ifpack2_SparseContainer_def.hpp:438
virtual std::string description() const
A one-line description of this object.
Definition Ifpack2_SparseContainer_def.hpp:448
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition Ifpack2_SparseContainer_def.hpp:79
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74