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