Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_BandedContainer_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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 
58 
59 namespace Ifpack2 {
60 
61 template<class MatrixType, class LocalScalarType>
63 BandedContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
64  const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
65  Container<MatrixType> (matrix, localRows),
66  numRows_ (localRows.size ()),
67  diagBlock_ (Teuchos::null),
68  ipiv_ (numRows_, 0),
69  kl_(-1),
70  ku_(-1)
71 {
72  typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
73  typedef typename Teuchos::ArrayView<const local_ordinal_type>::size_type size_type;
74  TEUCHOS_TEST_FOR_EXCEPTION(
75  ! matrix->hasColMap (), std::invalid_argument, "Ifpack2::BandedContainer: "
76  "The constructor's input matrix must have a column Map.");
77 
78  // Check whether the input set of local row indices is correct.
79  const map_type& rowMap = * (matrix->getRowMap ());
80  const size_type numRows = localRows.size ();
81  bool rowIndicesValid = true;
82  Teuchos::Array<local_ordinal_type> invalidLocalRowIndices;
83  for (size_type i = 0; i < numRows; ++i) {
84  if (! rowMap.isNodeLocalElement (localRows[i])) {
85  rowIndicesValid = false;
86  invalidLocalRowIndices.push_back (localRows[i]);
87  break;
88  }
89  }
90  TEUCHOS_TEST_FOR_EXCEPTION(
91  ! rowIndicesValid, std::invalid_argument, "Ifpack2::BandedContainer: "
92  "On process " << rowMap.getComm ()->getRank () << " of "
93  << rowMap.getComm ()->getSize () << ", in the given set of local row "
94  "indices localRows = " << Teuchos::toString (localRows) << ", the following "
95  "entries are not valid local row indices on the calling process: "
96  << Teuchos::toString (invalidLocalRowIndices) << ".");
97 
98 #ifdef HAVE_MPI
99  Teuchos::RCP<const Teuchos::Comm<int> > localComm =
100  Teuchos::rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF));
101 #else
102  Teuchos::RCP<const Teuchos::Comm<int> > localComm =
103  Teuchos::rcp (new Teuchos::SerialComm<int> ());
104 #endif // HAVE_MPI
105 
106  // FIXME (mfh 25 Aug 2013) What if the matrix's row Map has a
107  // different index base than zero?
108  const global_ordinal_type indexBase = 0;
109  localMap_ = rcp (new map_type (numRows_, indexBase, localComm));
110 }
111 
112 template<class MatrixType, class LocalScalarType>
114 BandedContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
115  const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
116  Container<MatrixType> (matrix, localRows)
117 {
118  TEUCHOS_TEST_FOR_EXCEPTION
119  (true, std::logic_error, "Ifpack2::BandedContainer: Not implemented for "
120  "LocalScalarType = " << Teuchos::TypeNameTraits<LocalScalarType>::name ()
121  << ".");
122 }
123 
124 
125 template<class MatrixType, class LocalScalarType>
128 {}
129 
130 template<class MatrixType, class LocalScalarType>
133 {}
134 
135 template<class MatrixType, class LocalScalarType>
137 setParameters (const Teuchos::ParameterList& List)
138 {
139  typedef typename Teuchos::ArrayView<const local_ordinal_type>::size_type size_type;
140 
141  if(List.isParameter("relaxation: banded container superdiagonals")) {
142  kl_ = List.get<int>("relaxation: banded container superdiagonals");
143  }
144  if(List.isParameter("relaxation: banded container subdiagonals")) {
145  ku_ = List.get<int>("relaxation: banded container subdiagonals");
146  }
147 
148  // The user provided insufficient information. If this is the case we check for the optimal values.
149  // User information may be overwritten only if necessary.
150  if (ku_ == -1 || kl_ == -1) {
151 
152  const Teuchos::ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
153 
154  const size_type numRows = localRows.size ();
155 
156  Teuchos::RCP<const row_matrix_type> matrix = this->getMatrix();
157 
158  // loop over local rows in current block
159  for (size_type i = 0; i < numRows; ++i) {
160  Teuchos::ArrayView< const local_ordinal_type > indices;
161  Teuchos::ArrayView< const scalar_type > values;
162  matrix->getLocalRowView(localRows[i], indices, values);
163 
164  size_type min_col_it = numRows > 0 ? numRows-1 : 0; // just a guess
165  size_type max_col_it = 0;
166 
167  size_type cntCols = 0;
168 
169  // loop over all column entries
170  for(size_type c = 0; c < indices.size(); c++) {
171  const local_ordinal_type lColIdx = indices[c]; // current column idx
172  // check whether lColIdx is contained in localRows[]
173  for (size_type j=0; j<numRows; ++j) {
174  if (localRows[j] == lColIdx) {
175  if (localRows[min_col_it] > lColIdx) {
176  min_col_it = j;
177  }
178  if (localRows[max_col_it] < lColIdx) {
179  max_col_it = j;
180  }
181  cntCols++;
182  }
183  }
184  if (cntCols == numRows) break; // skip remaining entries in column
185  }
186 
187  local_ordinal_type ku = Teuchos::as<local_ordinal_type>(max_col_it - i);
188  local_ordinal_type kl = Teuchos::as<local_ordinal_type>(i - min_col_it);
189 
190  if (ku > ku_) {
191  ku_ = ku; // update internal kl and ku with maximal values
192  }
193  if (kl > kl_) {
194  kl_ = kl;
195  }
196  }
197  }
198 
199  TEUCHOS_TEST_FOR_EXCEPTION
200  (kl_==-1 || ku_==-1, std::invalid_argument,
201  "Ifpack2::BandedContainer::setParameters: the user must provide the number"
202  " of sub- and superdiagonals in the 'kl' and 'ku' parameters.");
203 }
204 
205 template<class MatrixType, class LocalScalarType>
206 void
208 setParameters (const Teuchos::ParameterList& /* List */)
209 {
210  // Constructor of stub implementation throws, so no need to throw here.
211 }
212 
213 template<class MatrixType, class LocalScalarType>
214 void
217 {
218  using Teuchos::null;
219  using Teuchos::rcp;
220 
221  TEUCHOS_TEST_FOR_EXCEPTION
222  (kl_==-1 || ku_==-1, std::invalid_argument,
223  "Ifpack2::BandedContainer::initialize: the user must provide the number of"
224  " sub- and superdiagonals in the 'kl' and 'ku' parameters. Make sure that "
225  "you call BandedContainer<T>::setParameters!");
226 
227  //std::cout << "numRows_ = " << numRows_ << " kl_ = " << kl_ << " ku_ = " << ku_ << std::endl;
228  diagBlock_ = Teuchos::rcp(new Teuchos::SerialBandDenseMatrix<int, local_scalar_type>(numRows_, numRows_, kl_ /* lower bandwith */, kl_+ku_ /* upper bandwith. Don't forget to store kl extra superdiagonals for LU decomposition! */));
229 
230  // We assume that if you called this method, you intend to recompute
231  // everything.
232  IsInitialized_ = false;
233  IsComputed_ = false;
234 
235  // Fill the diagonal block and LU permutation array with zeros.
236  diagBlock_->putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
237  std::fill (ipiv_.begin (), ipiv_.end (), 0);
238 
239  IsInitialized_ = true;
240 }
241 
242 template<class MatrixType, class LocalScalarType>
243 void
246 {
247  // Constructor of stub implementation throws, so no need to throw here.
248 }
249 
250 template<class MatrixType, class LocalScalarType>
251 void
254 {
255  TEUCHOS_TEST_FOR_EXCEPTION(
256  static_cast<size_t> (ipiv_.size ()) != numRows_, std::logic_error,
257  "Ifpack2::BandedContainer::compute: ipiv_ array has the wrong size. "
258  "Please report this bug to the Ifpack2 developers.");
259 
260  IsComputed_ = false;
261  if (! this->isInitialized ()) {
262  this->initialize ();
263  }
264 
265  // Extract the submatrix.
266  extract (this->getMatrix ()); // extract the submatrix
267  factor (); // factor the submatrix
268 
269  IsComputed_ = true;
270 }
271 
272 template<class MatrixType, class LocalScalarType>
273 void
276 {
277  // Constructor of stub implementation throws, so no need to throw here.
278 }
279 
280 template<class MatrixType, class LocalScalarType>
281 void
283 factor ()
284 {
285  Teuchos::LAPACK<int, local_scalar_type> lapack;
286  int INFO = 0;
287 
288  // Plausibility checks for matrix
289  TEUCHOS_TEST_FOR_EXCEPTION(diagBlock_->values()==0, std::invalid_argument,
290  "BandedContainer<T>::factor: Diagonal block is an empty SerialBandDenseMatrix<T>!");
291  TEUCHOS_TEST_FOR_EXCEPTION(diagBlock_->upperBandwidth()<diagBlock_->lowerBandwidth(), std::invalid_argument,
292  "BandedContainer<T>::factor: Diagonal block needs kl additional superdiagonals for factorization! However, the number of superdiagonals is smaller than the number of subdiagonals!");
293 
294  lapack.GBTRF (diagBlock_->numRows(),
295  diagBlock_->numCols(),
296  diagBlock_->lowerBandwidth(),
297  diagBlock_->upperBandwidth()-diagBlock_->lowerBandwidth(), /* enter the real number of superdiagonals (see Teuchos_SerialBandDenseSolver)*/
298  diagBlock_->values(),
299  diagBlock_->stride(),
300  ipiv_.getRawPtr (),
301  &INFO);
302 
303  // INFO < 0 is a bug.
304  TEUCHOS_TEST_FOR_EXCEPTION(
305  INFO < 0, std::logic_error, "Ifpack2::BandedContainer::factor: "
306  "LAPACK's _GBTRF (LU factorization with partial pivoting) was called "
307  "incorrectly. INFO = " << INFO << " < 0. "
308  "Please report this bug to the Ifpack2 developers.");
309  // INFO > 0 means the matrix is singular. This is probably an issue
310  // either with the choice of rows the rows we extracted, or with the
311  // input matrix itself.
312  TEUCHOS_TEST_FOR_EXCEPTION(
313  INFO > 0, std::runtime_error, "Ifpack2::BandedContainer::factor: "
314  "LAPACK's _GBTRF (LU factorization with partial pivoting) reports that the "
315  "computed U factor is exactly singular. U(" << INFO << "," << INFO << ") "
316  "(one-based index i) is exactly zero. This probably means that the input "
317  "matrix has a singular diagonal block.");
318 }
319 
320 template<class MatrixType, class LocalScalarType>
321 void
323 factor ()
324 {
325  // Constructor of stub implementation throws, so no need to throw here.
326 }
327 
328 template<class MatrixType, class LocalScalarType>
329 void
331 applyImpl (const local_mv_type& X,
332  local_mv_type& Y,
333  Teuchos::ETransp mode,
334  LocalScalarType alpha,
335  LocalScalarType beta) const
336 {
337  using Teuchos::ArrayRCP;
338  using Teuchos::RCP;
339  using Teuchos::rcp;
340  using Teuchos::rcpFromRef;
341 
342  TEUCHOS_TEST_FOR_EXCEPTION(
343  X.getLocalLength () != Y.getLocalLength (),
344  std::logic_error, "Ifpack2::BandedContainer::applyImpl: X and Y have "
345  "incompatible dimensions (" << X.getLocalLength () << " resp. "
346  << Y.getLocalLength () << "). Please report this bug to "
347  "the Ifpack2 developers.");
348  TEUCHOS_TEST_FOR_EXCEPTION(
349  localMap_->getNodeNumElements () != X.getLocalLength (),
350  std::logic_error, "Ifpack2::BandedContainer::applyImpl: The inverse "
351  "operator and X have incompatible dimensions (" <<
352  localMap_->getNodeNumElements () << " resp. "
353  << X.getLocalLength () << "). Please report this bug to "
354  "the Ifpack2 developers.");
355  TEUCHOS_TEST_FOR_EXCEPTION(
356  localMap_->getNodeNumElements () != Y.getLocalLength (),
357  std::logic_error, "Ifpack2::BandedContainer::applyImpl: The inverse "
358  "operator and Y have incompatible dimensions (" <<
359  localMap_->getNodeNumElements () << " resp. "
360  << Y.getLocalLength () << "). Please report this bug to "
361  "the Ifpack2 developers.");
362  TEUCHOS_TEST_FOR_EXCEPTION(
363  X.getLocalLength() != static_cast<size_t> (mode == Teuchos::NO_TRANS ? diagBlock_->numCols() : diagBlock_->numRows()),
364  std::logic_error, "Ifpack2::BandedContainer::applyImpl: The input "
365  "multivector X has incompatible dimensions from those of the "
366  "inverse operator (" << X.getLocalLength () << " vs. "
367  << (mode == Teuchos::NO_TRANS ? diagBlock_->numCols() : diagBlock_->numRows())
368  << "). Please report this bug to the Ifpack2 developers.");
369  TEUCHOS_TEST_FOR_EXCEPTION(
370  Y.getLocalLength() != static_cast<size_t> (mode == Teuchos::NO_TRANS ? diagBlock_->numRows() : diagBlock_->numCols()),
371  std::logic_error, "Ifpack2::BandedContainer::applyImpl: The output "
372  "multivector Y has incompatible dimensions from those of the "
373  "inverse operator (" << Y.getLocalLength () << " vs. "
374  << (mode == Teuchos::NO_TRANS ? diagBlock_->numRows() : diagBlock_->numCols())
375  << "). Please report this bug to the Ifpack2 developers.");
376 
377  typedef Teuchos::ScalarTraits<local_scalar_type> STS;
378  const int numVecs = static_cast<int> (X.getNumVectors ());
379  if (alpha == STS::zero ()) { // don't need to solve the linear system
380  if (beta == STS::zero ()) {
381  // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
382  // any Inf or NaN values in Y (rather than multiplying them by
383  // zero, resulting in NaN values).
384  Y.putScalar (STS::zero ());
385  }
386  else { // beta != 0
387  Y.scale (beta);
388  }
389  }
390  else { // alpha != 0; must solve the linear system
391  Teuchos::LAPACK<int, local_scalar_type> lapack;
392  // If beta is nonzero or Y is not constant stride, we have to use
393  // a temporary output multivector. It gets a copy of X, since
394  // GBTRS overwrites its (multi)vector input with its output.
395  RCP<local_mv_type> Y_tmp;
396  if (beta == STS::zero () ){
397  Y = X;
398  Y_tmp = rcpFromRef (Y);
399  }
400  else {
401  Y_tmp = rcp (new local_mv_type (X)); // constructor copies X
402  }
403  const int Y_stride = static_cast<int> (Y_tmp->getStride ());
404  ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
405  local_scalar_type* const Y_ptr = Y_view.getRawPtr ();
406 
407  int INFO = 0;
408  const char trans =(mode == Teuchos::CONJ_TRANS ? 'C' : (mode == Teuchos::TRANS ? 'T' : 'N'));
409 
410  lapack.GBTRS(trans,
411  diagBlock_->numCols(),
412  diagBlock_->lowerBandwidth(),
413  diagBlock_->upperBandwidth()-diagBlock_->lowerBandwidth(), // /* enter the real number of superdiagonals (see Teuchos_SerialBandDenseSolver)*/
414  numVecs,
415  diagBlock_->values(),
416  diagBlock_->stride(),
417  ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
418 
419  TEUCHOS_TEST_FOR_EXCEPTION(
420  INFO != 0, std::runtime_error, "Ifpack2::BandedContainer::applyImpl: "
421  "LAPACK's _GBTRS (solve using LU factorization with partial pivoting) "
422  "failed with INFO = " << INFO << " != 0.");
423 
424  if (beta != STS::zero ()) {
425  Y.update (alpha, *Y_tmp, beta);
426  }
427  else if (! Y.isConstantStride ()) {
428  Tpetra::deep_copy (Y, *Y_tmp);
429  }
430  }
431 }
432 
433 template<class MatrixType, class LocalScalarType>
434 void
436 applyImpl (const local_mv_type& X,
437  local_mv_type& Y,
438  Teuchos::ETransp mode,
439  LocalScalarType alpha,
440  LocalScalarType beta) const
441 {
442  // Constructor of stub implementation throws, so no need to throw here.
443 }
444 
445 template<class MatrixType, class LocalScalarType>
446 void
448 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
449  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
450  Teuchos::ETransp mode,
451  scalar_type alpha,
452  scalar_type beta) const
453 {
454  using Teuchos::ArrayView;
455  using Teuchos::as;
456  using Teuchos::RCP;
457  using Teuchos::rcp;
458 
459  // The local operator might have a different Scalar type than
460  // MatrixType. This means that we might have to convert X and Y to
461  // the Tpetra::MultiVector specialization that the local operator
462  // wants. This class' X_ and Y_ internal fields are of the right
463  // type for the local operator, so we can use those as targets.
464 
465  // Tpetra::MultiVector specialization corresponding to MatrixType.
466  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
467  global_ordinal_type, node_type> global_mv_type;
469  const size_t numVecs = X.getNumVectors ();
470 
471  TEUCHOS_TEST_FOR_EXCEPTION(
472  ! IsComputed_, std::runtime_error, "Ifpack2::BandedContainer::apply: "
473  "You must have called the compute() method before you may call apply(). "
474  "You may call the apply() method as many times as you want after calling "
475  "compute() once, but you must have called compute() at least once.");
476  TEUCHOS_TEST_FOR_EXCEPTION(
477  numVecs != Y.getNumVectors (), std::runtime_error,
478  "Ifpack2::BandedContainer::apply: X and Y have different numbers of "
479  "vectors. X has " << X.getNumVectors ()
480  << ", but Y has " << X.getNumVectors () << ".");
481 
482  if (numVecs == 0) {
483  return; // done! nothing to do
484  }
485 
486  // The local operator works on a permuted subset of the local parts
487  // of X and Y. The subset and permutation are defined by the index
488  // array returned by getLocalRows(). If the permutation is trivial
489  // and the subset is exactly equal to the local indices, then we
490  // could use the local parts of X and Y exactly, without needing to
491  // permute. Otherwise, we have to use temporary storage to permute
492  // X and Y. For now, we always use temporary storage.
493  //
494  // Create temporary permuted versions of the input and output.
495  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
496  // store the permuted versions of X resp. Y. Note that X_local has
497  // the domain Map of the operator, which may be a permuted subset of
498  // the local Map corresponding to X.getMap(). Similarly, Y_local
499  // has the range Map of the operator, which may be a permuted subset
500  // of the local Map corresponding to Y.getMap(). numRows_ here
501  // gives the number of rows in the row Map of the local Inverse_
502  // operator.
503  //
504  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
505  // here that the row Map and the range Map of that operator are
506  // the same.
507  //
508  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
509  // really belongs in Tpetra.
510 
511  if (X_.is_null ()) {
512  X_ = rcp (new local_mv_type (localMap_, numVecs));
513  }
514  RCP<local_mv_type> X_local = X_;
515  TEUCHOS_TEST_FOR_EXCEPTION(
516  X_local->getLocalLength () != numRows_, std::logic_error,
517  "Ifpack2::BandedContainer::apply: "
518  "X_local has length " << X_local->getLocalLength () << ", which does "
519  "not match numRows_ = " << numRows_ << ". Please report this bug to "
520  "the Ifpack2 developers.");
521  ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
522  mvgs.gather (*X_local, X, localRows);
523 
524  // We must gather the contents of the output multivector Y even on
525  // input to applyImpl(), since the inverse operator might use it as
526  // an initial guess for a linear solve. We have no way of knowing
527  // whether it does or does not.
528 
529  if (Y_.is_null ()) {
530  Y_ = rcp (new local_mv_type (localMap_, numVecs));
531  }
532  RCP<local_mv_type> Y_local = Y_;
533  TEUCHOS_TEST_FOR_EXCEPTION(
534  Y_local->getLocalLength () != numRows_, std::logic_error,
535  "Ifpack2::BandedContainer::apply: "
536  "Y_local has length " << X_local->getLocalLength () << ", which does "
537  "not match numRows_ = " << numRows_ << ". Please report this bug to "
538  "the Ifpack2 developers.");
539  mvgs.gather (*Y_local, Y, localRows);
540 
541  // Apply the local operator:
542  // Y_local := beta*Y_local + alpha*M^{-1}*X_local
543  this->applyImpl (*X_local, *Y_local, mode, as<local_scalar_type> (alpha),
544  as<local_scalar_type> (beta));
545 
546  // Scatter the permuted subset output vector Y_local back into the
547  // original output multivector Y.
548  mvgs.scatter (Y, *Y_local, localRows);
549 }
550 
551 template<class MatrixType, class LocalScalarType>
552 void
554 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& /* X */,
555  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& /* Y */,
556  Teuchos::ETransp /* mode */,
557  scalar_type /* alpha */,
558  scalar_type /* beta */) const
559 {
560  // Constructor of stub implementation throws, so no need to throw here.
561 }
562 
563 template<class MatrixType, class LocalScalarType>
564 void
566 weightedApply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
567  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
568  const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
569  Teuchos::ETransp mode,
570  scalar_type alpha,
571  scalar_type beta) const
572 {
573  using Teuchos::ArrayRCP;
574  using Teuchos::ArrayView;
575  using Teuchos::Range1D;
576  using Teuchos::RCP;
577  using Teuchos::rcp;
578  using Teuchos::rcp_const_cast;
579  using std::cerr;
580  using std::endl;
581  typedef Teuchos::ScalarTraits<scalar_type> STS;
582 
583  TEUCHOS_TEST_FOR_EXCEPTION(
584  true, std::runtime_error, "Ifpack2::BandedContainer::"
585  "weightedApply: This code is not tested and not used. Expect bugs.");
586 
587  // The local operator template parameter might have a different
588  // Scalar type than MatrixType. This means that we might have to
589  // convert X and Y to the Tpetra::MultiVector specialization that
590  // the local operator wants. This class' X_ and Y_ internal fields
591  // are of the right type for the local operator, so we can use those
592  // as targets.
593 
594  // Tpetra::MultiVector specialization corresponding to the global operator.
595  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
596  global_ordinal_type, node_type> global_mv_type;
597  // Tpetra::Vector specialization corresponding to the local
598  // operator. We will use this for the subset permutation of the
599  // diagonal scaling D.
600  typedef Tpetra::Vector<local_scalar_type, local_ordinal_type,
601  global_ordinal_type, node_type> local_vec_type;
602 
604  const size_t numVecs = X.getNumVectors ();
605 
606  TEUCHOS_TEST_FOR_EXCEPTION(
607  ! IsComputed_, std::runtime_error, "Ifpack2::BandedContainer::"
608  "weightedApply: You must have called the compute() method before you may "
609  "call apply(). You may call the apply() method as many times as you want "
610  "after calling compute() once, but you must have called compute() at least "
611  "once.");
612  TEUCHOS_TEST_FOR_EXCEPTION(
613  numVecs != Y.getNumVectors (), std::runtime_error,
614  "Ifpack2::BandedContainer::weightedApply: X and Y have different numbers "
615  "of vectors. X has " << X.getNumVectors () << ", but Y has "
616  << X.getNumVectors () << ".");
617 
618  if (numVecs == 0) {
619  return; // done! nothing to do
620  }
621 
622  // The local operator works on a permuted subset of the local parts
623  // of X and Y. The subset and permutation are defined by the index
624  // array returned by getLocalRows(). If the permutation is trivial
625  // and the subset is exactly equal to the local indices, then we
626  // could use the local parts of X and Y exactly, without needing to
627  // permute. Otherwise, we have to use temporary storage to permute
628  // X and Y. For now, we always use temporary storage.
629  //
630  // Create temporary permuted versions of the input and output.
631  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
632  // store the permuted versions of X resp. Y. Note that X_local has
633  // the domain Map of the operator, which may be a permuted subset of
634  // the local Map corresponding to X.getMap(). Similarly, Y_local
635  // has the range Map of the operator, which may be a permuted subset
636  // of the local Map corresponding to Y.getMap(). numRows_ here
637  // gives the number of rows in the row Map of the local operator.
638  //
639  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
640  // here that the row Map and the range Map of that operator are
641  // the same.
642  //
643  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
644  // really belongs in Tpetra.
645 
646  if (X_.is_null ()) {
647  X_ = rcp (new local_mv_type (localMap_, numVecs));
648  }
649  RCP<local_mv_type> X_local = X_;
650  TEUCHOS_TEST_FOR_EXCEPTION(
651  X_local->getLocalLength () != numRows_, std::logic_error,
652  "Ifpack2::BandedContainer::weightedApply: "
653  "X_local has length " << X_local->getLocalLength () << ", which does "
654  "not match numRows_ = " << numRows_ << ". Please report this bug to "
655  "the Ifpack2 developers.");
656  ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
657  mvgs.gather (*X_local, X, localRows);
658 
659  // We must gather the output multivector Y even on input to
660  // applyImpl(), since the local operator might use it as an initial
661  // guess for a linear solve. We have no way of knowing whether it
662  // does or does not.
663 
664  if (Y_.is_null ()) {
665  Y_ = rcp (new local_mv_type (localMap_, numVecs));
666  }
667  RCP<local_mv_type> Y_local = Y_;
668  TEUCHOS_TEST_FOR_EXCEPTION(
669  Y_local->getLocalLength () != numRows_, std::logic_error,
670  "Ifpack2::BandedContainer::weightedApply: "
671  "Y_local has length " << X_local->getLocalLength () << ", which does "
672  "not match numRows_ = " << numRows_ << ". Please report this bug to "
673  "the Ifpack2 developers.");
674  mvgs.gather (*Y_local, Y, localRows);
675 
676  // Apply the diagonal scaling D to the input X. It's our choice
677  // whether the result has the original input Map of X, or the
678  // permuted subset Map of X_local. If the latter, we also need to
679  // gather D into the permuted subset Map. We choose the latter, to
680  // save memory and computation. Thus, we do the following:
681  //
682  // 1. Gather D into a temporary vector D_local.
683  // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
684  // 3. Compute X_scaled := diag(D_loca) * X_local.
685 
686  RCP<local_vec_type> D_local = rcp (new local_vec_type (localMap_));
687  TEUCHOS_TEST_FOR_EXCEPTION(
688  D_local->getLocalLength () != numRows_, std::logic_error,
689  "Ifpack2::BandedContainer::weightedApply: "
690  "D_local has length " << X_local->getLocalLength () << ", which does "
691  "not match numRows_ = " << numRows_ << ". Please report this bug to "
692  "the Ifpack2 developers.");
693  mvgs.gather (*D_local, D, localRows);
694  RCP<local_mv_type> X_scaled = rcp (new local_mv_type (localMap_, numVecs));
695  X_scaled->elementWiseMultiply (STS::one (), *D_local, *X_local, STS::zero ());
696 
697  // Y_temp will hold the result of M^{-1}*X_scaled. If beta == 0, we
698  // can write the result of Inverse_->apply() directly to Y_local, so
699  // Y_temp may alias Y_local. Otherwise, if beta != 0, we need
700  // temporary storage for M^{-1}*X_scaled, so Y_temp must be
701  // different than Y_local.
702  RCP<local_mv_type> Y_temp;
703  if (beta == STS::zero ()) {
704  Y_temp = Y_local;
705  } else {
706  Y_temp = rcp (new local_mv_type (localMap_, numVecs));
707  }
708 
709  // Apply the local operator: Y_temp := M^{-1} * X_scaled
710  applyImpl (*X_scaled, *Y_temp, mode, STS::one (), STS::zero ());
711  // Y_local := beta * Y_local + alpha * diag(D_local) * Y_temp.
712 
713  //
714  // Note that we still use the permuted subset scaling D_local here,
715  // because Y_temp has the same permuted subset Map. That's good, in
716  // fact, because it's a subset: less data to read and multiply.
717  Y_local->elementWiseMultiply (alpha, *D_local, *Y_temp, beta);
718 
719  // Copy the permuted subset output vector Y_local into the original
720  // output multivector Y.
721  mvgs.scatter (Y, *Y_local, localRows);
722 }
723 
724 template<class MatrixType, class LocalScalarType>
725 void
727 weightedApply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& /* X */,
728  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& /* Y */,
729  const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& /* D */,
730  Teuchos::ETransp /* mode */,
731  scalar_type /* alpha */,
732  scalar_type /* beta */) const
733 {
734  // Constructor of stub implementation throws, so no need to throw here.
735 }
736 
737 template<class MatrixType, class LocalScalarType>
738 std::ostream&
740 print (std::ostream& os) const
741 {
742  Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
743  fos.setOutputToRootOnly (0);
744  describe (fos);
745  return os;
746 }
747 
748 template<class MatrixType, class LocalScalarType>
749 std::ostream&
751 print (std::ostream& os) const
752 {
753  return os;
754  // Constructor of stub implementation throws, so no need to throw here.
755 }
756 
757 template<class MatrixType, class LocalScalarType>
758 std::string
760 description () const
761 {
762  std::ostringstream oss;
763  oss << Teuchos::Describable::description();
764  if (isInitialized()) {
765  if (isComputed()) {
766  oss << "{status = initialized, computed";
767  }
768  else {
769  oss << "{status = initialized, not computed";
770  }
771  }
772  else {
773  oss << "{status = not initialized, not computed";
774  }
775 
776  oss << "}";
777  return oss.str();
778 }
779 
780 template<class MatrixType, class LocalScalarType>
781 std::string
783 description () const
784 {
785  return "";
786  // Constructor of stub implementation throws, so no need to throw here.
787 }
788 
789 template<class MatrixType, class LocalScalarType>
790 void
792 describe (Teuchos::FancyOStream& os,
793  const Teuchos::EVerbosityLevel verbLevel) const
794 {
795  if(verbLevel==Teuchos::VERB_NONE) return;
796  os << "================================================================================" << std::endl;
797  os << "Ifpack2::BandedContainer" << std::endl;
798  os << "Number of rows = " << numRows_ << std::endl;
799  os << "Number of subdiagonals = " << diagBlock_->lowerBandwidth() << std::endl;
800  os << "Number of superdiagonals = " << diagBlock_->upperBandwidth() << std::endl;
801  os << "isInitialized() = " << IsInitialized_ << std::endl;
802  os << "isComputed() = " << IsComputed_ << std::endl;
803  os << "================================================================================" << std::endl;
804  os << std::endl;
805 }
806 
807 template<class MatrixType, class LocalScalarType>
808 void
810 describe (Teuchos::FancyOStream& os,
811  const Teuchos::EVerbosityLevel verbLevel) const
812 {
813  // Constructor of stub implementation throws, so no need to throw here.
814 }
815 
816 template<class MatrixType, class LocalScalarType>
817 void
819 extract (const Teuchos::RCP<const row_matrix_type>& globalMatrix)
820 {
821  using Teuchos::Array;
822  using Teuchos::ArrayView;
823  using Teuchos::toString;
824  typedef local_ordinal_type LO;
825  typedef global_ordinal_type GO;
826  typedef Tpetra::Map<LO, GO, node_type> map_type;
827  const size_t inputMatrixNumRows = globalMatrix->getNodeNumRows ();
828  // We only use the rank of the calling process and the number of MPI
829  // processes for generating error messages. Extraction itself is
830  // entirely local to each participating MPI process.
831  const int myRank = globalMatrix->getRowMap ()->getComm ()->getRank ();
832  const int numProcs = globalMatrix->getRowMap ()->getComm ()->getSize ();
833 
834  // Sanity check that the local row indices to extract fall within
835  // the valid range of local row indices for the input matrix.
836  ArrayView<const LO> localRows = this->getLocalRows ();
837  for (size_t j = 0; j < numRows_; ++j) {
838  TEUCHOS_TEST_FOR_EXCEPTION(
839  localRows[j] < 0 ||
840  static_cast<size_t> (localRows[j]) >= inputMatrixNumRows,
841  std::runtime_error, "Ifpack2::BandedContainer::extract: On process " <<
842  myRank << " of " << numProcs << ", localRows[j=" << j << "] = " <<
843  localRows[j] << ", which is out of the valid range of local row indices "
844  "indices [0, " << (inputMatrixNumRows - 1) << "] for the input matrix.");
845  }
846 
847  // Convert the local row indices we want into local column indices.
848  // For every local row ii_local = localRows[i] we take, we also want
849  // to take the corresponding column. To find the corresponding
850  // column, we use the row Map to convert the local row index
851  // ii_local into a global index ii_global, and then use the column
852  // Map to convert ii_global into a local column index jj_local. If
853  // the input matrix doesn't have a column Map, we need to be using
854  // global indices anyway...
855 
856  // We use the domain Map to exclude off-process global entries.
857  const map_type& globalRowMap = * (globalMatrix->getRowMap ());
858  const map_type& globalColMap = * (globalMatrix->getColMap ());
859  const map_type& globalDomMap = * (globalMatrix->getDomainMap ());
860 
861  bool rowIndsValid = true;
862  bool colIndsValid = true;
863  Array<LO> localCols (numRows_);
864  // For error messages, collect the sets of invalid row indices and
865  // invalid column indices. They are otherwise not useful.
866  Array<LO> invalidLocalRowInds;
867  Array<GO> invalidGlobalColInds;
868  for (size_t i = 0; i < numRows_; ++i) {
869  // ii_local is the (local) row index we want to look up.
870  const LO ii_local = localRows[i];
871  // Find the global index jj_global corresponding to ii_local.
872  // Global indices are the same (rather, are required to be the
873  // same) in all three Maps, which is why we use jj (suggesting a
874  // column index, which is how we will use it below).
875  const GO jj_global = globalRowMap.getGlobalElement (ii_local);
876  if (jj_global == Teuchos::OrdinalTraits<GO>::invalid ()) {
877  // If ii_local is not a local index in the row Map on the
878  // calling process, that means localRows is incorrect. We've
879  // already checked for this in the constructor, but we might as
880  // well check again here, since it's cheap to do so (just an
881  // integer comparison, since we need jj_global anyway).
882  rowIndsValid = false;
883  invalidLocalRowInds.push_back (ii_local);
884  break;
885  }
886  // Exclude "off-process" entries: that is, those in the column Map
887  // on this process that are not in the domain Map on this process.
888  if (globalDomMap.isNodeGlobalElement (jj_global)) {
889  // jj_global is not an off-process entry. Look up its local
890  // index in the column Map; we want to extract this column index
891  // from the input matrix. If jj_global is _not_ in the column
892  // Map on the calling process, that could mean that the column
893  // in question is empty on this process. That would be bad for
894  // solving linear systems with the extract submatrix. We could
895  // solve the resulting singular linear systems in a minimum-norm
896  // least-squares sense, but for now we simply raise an exception.
897  const LO jj_local = globalColMap.getLocalElement (jj_global);
898  if (jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
899  colIndsValid = false;
900  invalidGlobalColInds.push_back (jj_global);
901  break;
902  }
903  localCols[i] = jj_local;
904  }
905  }
906  TEUCHOS_TEST_FOR_EXCEPTION(
907  ! rowIndsValid, std::logic_error, "Ifpack2::BandedContainer::extract: "
908  "On process " << myRank << ", at least one row index in the set of local "
909  "row indices given to the constructor is not a valid local row index in "
910  "the input matrix's row Map on this process. This should be impossible "
911  "because the constructor checks for this case. Here is the complete set "
912  "of invalid local row indices: " << toString (invalidLocalRowInds) << ". "
913  "Please report this bug to the Ifpack2 developers.");
914  TEUCHOS_TEST_FOR_EXCEPTION(
915  ! colIndsValid, std::runtime_error, "Ifpack2::BandedContainer::extract: "
916  "On process " << myRank << ", "
917  "At least one row index in the set of row indices given to the constructor "
918  "does not have a corresponding column index in the input matrix's column "
919  "Map. This probably means that the column(s) in question is/are empty on "
920  "this process, which would make the submatrix to extract structurally "
921  "singular. Here is the compete set of invalid global column indices: "
922  << toString (invalidGlobalColInds) << ".");
923 
924  diagBlock_->putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
925 
926  const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries ();
927  Array<scalar_type> val (maxNumEntriesInRow);
928  Array<local_ordinal_type> ind (maxNumEntriesInRow);
929 
930  const local_ordinal_type INVALID =
931  Teuchos::OrdinalTraits<local_ordinal_type>::invalid ();
932  for (size_t i = 0; i < numRows_; ++i) {
933  const local_ordinal_type localRow = localRows[i];
934  size_t numEntries;
935  globalMatrix->getLocalRowCopy (localRow, ind (), val (), numEntries);
936 
937  for (size_t k = 0; k < numEntries; ++k) {
938  const local_ordinal_type localCol = ind[k];
939 
940  // Skip off-process elements
941  //
942  // FIXME (mfh 24 Aug 2013) This assumes the following:
943  //
944  // 1. The column and row Maps begin with the same set of
945  // on-process entries, in the same order. That is,
946  // on-process row and column indices are the same.
947  // 2. All off-process indices in the column Map of the input
948  // matrix occur after that initial set.
949  if (localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows) {
950  // for local column IDs, look for each ID in the list
951  // of columns hosted by this object
952  local_ordinal_type jj = INVALID;
953  for (size_t kk = 0; kk < numRows_; ++kk) {
954  if (localRows[kk] == localCol) {
955  jj = kk;
956  }
957  }
958  if (jj != INVALID) {
959  (*diagBlock_) (i, jj) += val[k]; // ???
960  }
961  }
962  }
963  }
964 }
965 
966 template<class MatrixType, class LocalScalarType>
967 void
969 extract (const Teuchos::RCP<const row_matrix_type>& /* globalMatrix */)
970 {
971  // Constructor of stub implementation throws, so no need to throw here.
972 }
973 
974 } // namespace Ifpack2
975 
976 #define IFPACK2_BANDEDCONTAINER_INSTANT(S,LO,GO,N) \
977  template class Ifpack2::BandedContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >;
978 
979 #endif // IFPACK2_BANDEDCONTAINER_HPP
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:131
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:340
MatrixType::node_type node_type
The Node type of the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:137
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_BandedContainer_decl.hpp:333
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:133
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:135
BandedContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::ArrayView< const local_ordinal_type > &localRows)
Constructor.
Definition: Ifpack2_BandedContainer_def.hpp:63
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:338
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:103
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_BandedContainer_decl.hpp:128
BELOS_DEPRECATED const char * toString(const StatusType status)
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Store and solve a local Banded linear problem.
Definition: Ifpack2_BandedContainer_decl.hpp:108
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_BandedContainer_decl.hpp:336