Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_OverlappingRowMatrix_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #ifndef IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
00044 #define IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
00045 
00046 #include <Ifpack2_OverlappingRowMatrix_decl.hpp>
00047 #include <Ifpack2_Details_OverlappingRowGraph.hpp>
00048 #include <Tpetra_CrsMatrix.hpp>
00049 #include <Teuchos_CommHelpers.hpp>
00050 
00051 namespace Ifpack2 {
00052 
00053 template<class MatrixType>
00054 OverlappingRowMatrix<MatrixType>::
00055 OverlappingRowMatrix (const Teuchos::RCP<const row_matrix_type>& A,
00056                       const int overlapLevel) :
00057   A_ (A),
00058   OverlapLevel_ (overlapLevel),
00059   UseSubComm_ (false)
00060 {
00061   using Teuchos::RCP;
00062   using Teuchos::rcp;
00063   using Teuchos::Array;
00064   using Teuchos::outArg;
00065   using Teuchos::rcp_const_cast;
00066   using Teuchos::rcp_dynamic_cast;
00067   using Teuchos::rcp_implicit_cast;
00068   using Teuchos::REDUCE_SUM;
00069   using Teuchos::reduceAll;
00070   typedef Tpetra::global_size_t GST;
00071   typedef Tpetra::CrsGraph<local_ordinal_type,
00072                            global_ordinal_type, node_type> crs_graph_type;
00073   TEUCHOS_TEST_FOR_EXCEPTION(
00074     OverlapLevel_ <= 0, std::runtime_error,
00075     "Ifpack2::OverlappingRowMatrix: OverlapLevel must be > 0.");
00076   TEUCHOS_TEST_FOR_EXCEPTION(
00077     A_->getComm()->getSize() == 1, std::runtime_error,
00078     "Ifpack2::OverlappingRowMatrix: Matrix must be "
00079     "distributed over more than one MPI process.");
00080 
00081   RCP<const crs_matrix_type> ACRS =
00082     rcp_dynamic_cast<const crs_matrix_type, const row_matrix_type> (A_);
00083   TEUCHOS_TEST_FOR_EXCEPTION(
00084     ACRS.is_null (), std::runtime_error,
00085     "Ifpack2::OverlappingRowMatrix: The input matrix must be a Tpetra::"
00086     "CrsMatrix with matching template parameters.  This class currently "
00087     "requires that CrsMatrix's fifth template parameter be the default.");
00088   RCP<const crs_graph_type> A_crsGraph = ACRS->getCrsGraph ();
00089 
00090   const size_t numMyRowsA = A_->getNodeNumRows ();
00091   const global_ordinal_type global_invalid =
00092     Teuchos::OrdinalTraits<global_ordinal_type>::invalid ();
00093 
00094   // Temp arrays
00095   Array<global_ordinal_type> ExtElements;
00096   RCP<map_type>        TmpMap;
00097   RCP<crs_graph_type>  TmpGraph;
00098   RCP<import_type>     TmpImporter;
00099   RCP<const map_type>  RowMap, ColMap;
00100 
00101   // The big import loop
00102   for (int overlap = 0 ; overlap < OverlapLevel_ ; ++overlap) {
00103     // Get the current maps
00104     if (overlap == 0) {
00105       RowMap = A_->getRowMap ();
00106       ColMap = A_->getColMap ();
00107     }
00108     else {
00109       RowMap = TmpGraph->getRowMap ();
00110       ColMap = TmpGraph->getColMap ();
00111     }
00112 
00113     const size_t size = ColMap->getNodeNumElements () - RowMap->getNodeNumElements ();
00114     Array<global_ordinal_type> mylist (size);
00115     size_t count = 0;
00116 
00117     // define the set of rows that are in ColMap but not in RowMap
00118     for (local_ordinal_type i = 0 ; (size_t) i < ColMap->getNodeNumElements() ; ++i) {
00119       const global_ordinal_type GID = ColMap->getGlobalElement (i);
00120       if (A_->getRowMap ()->getLocalElement (GID) == global_invalid) {
00121         typedef typename Array<global_ordinal_type>::iterator iter_type;
00122         const iter_type end = ExtElements.end ();
00123         const iter_type pos = std::find (ExtElements.begin (), end, GID);
00124         if (pos == end) {
00125           ExtElements.push_back (GID);
00126           mylist[count] = GID;
00127           ++count;
00128         }
00129       }
00130     }
00131 
00132     // mfh 24 Nov 2013: We don't need TmpMap, TmpGraph, or
00133     // TmpImporter after this loop, so we don't have to construct them
00134     // on the last round.
00135     if (overlap + 1 < OverlapLevel_) {
00136       // Allocate & import new matrices, maps, etc.
00137       //
00138       // FIXME (mfh 24 Nov 2013) Do we always want to use index base
00139       // zero?  It doesn't really matter, since the actual index base
00140       // (in the current implementation of Map) will always be the
00141       // globally least GID.
00142       TmpMap = rcp (new map_type (global_invalid, mylist (0, count),
00143                                   Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
00144                                   A_->getComm (), A_->getNode ()));
00145       TmpGraph = rcp (new crs_graph_type (TmpMap, 0));
00146       TmpImporter = rcp (new import_type (A_->getRowMap (), TmpMap));
00147 
00148       TmpGraph->doImport (*A_crsGraph, *TmpImporter, Tpetra::INSERT);
00149       TmpGraph->fillComplete (A_->getDomainMap (), TmpMap);
00150     }
00151   }
00152 
00153   // build the map containing all the nodes (original
00154   // matrix + extended matrix)
00155   Array<global_ordinal_type> mylist (numMyRowsA + ExtElements.size ());
00156   for (local_ordinal_type i = 0; (size_t)i < numMyRowsA; ++i) {
00157     mylist[i] = A_->getRowMap ()->getGlobalElement (i);
00158   }
00159   for (local_ordinal_type i = 0; i < ExtElements.size (); ++i) {
00160     mylist[i + numMyRowsA] = ExtElements[i];
00161   }
00162 
00163   RowMap_ = rcp (new map_type (global_invalid, mylist (),
00164                                Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
00165                                A_->getComm (), A_->getNode ()));
00166   ColMap_ = RowMap_;
00167 
00168   // now build the map corresponding to all the external nodes
00169   // (with respect to A().RowMatrixRowMap().
00170   ExtMap_ = rcp (new map_type (global_invalid, ExtElements (),
00171                                Teuchos::OrdinalTraits<global_ordinal_type>::zero (),
00172                                A_->getComm (), A_->getNode ()));
00173   ExtMatrix_ = rcp (new crs_matrix_type (ExtMap_, ColMap_, 0));
00174   ExtImporter_ = rcp (new import_type (A_->getRowMap (), ExtMap_));
00175 
00176   RCP<crs_matrix_type> ExtMatrixCRS =
00177     rcp_dynamic_cast<crs_matrix_type, row_matrix_type> (ExtMatrix_);
00178   ExtMatrixCRS->doImport (*ACRS, *ExtImporter_, Tpetra::INSERT);
00179   ExtMatrixCRS->fillComplete (A_->getDomainMap (), RowMap_);
00180 
00181   Importer_ = rcp (new import_type (A_->getRowMap (), RowMap_));
00182 
00183   // fix indices for overlapping matrix
00184   const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
00185 
00186   GST NumMyNonzeros_tmp = A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
00187   GST NumMyRows_tmp = numMyRowsA + numMyRowsB;
00188   {
00189     GST inArray[2], outArray[2];
00190     inArray[0] = NumMyNonzeros_tmp;
00191     inArray[1] = NumMyRows_tmp;
00192     outArray[0] = 0;
00193     outArray[1] = 0;
00194     reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, 2, inArray, outArray);
00195     NumGlobalNonzeros_ = outArray[0];
00196     NumGlobalRows_ = outArray[1];
00197   }
00198   // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyNonzeros_tmp,
00199   //                      outArg (NumGlobalNonzeros_));
00200   // reduceAll<int, GST> (* (A_->getComm ()), REDUCE_SUM, NumMyRows_tmp,
00201   //                      outArg (NumGlobalRows_));
00202 
00203   MaxNumEntries_ = A_->getNodeMaxNumRowEntries ();
00204   if (MaxNumEntries_ < ExtMatrix_->getNodeMaxNumRowEntries ()) {
00205     MaxNumEntries_ = ExtMatrix_->getNodeMaxNumRowEntries ();
00206   }
00207 
00208   // Create the graph (returned by getGraph()).
00209   typedef Details::OverlappingRowGraph<row_graph_type> row_graph_impl_type;
00210   RCP<row_graph_impl_type> graph =
00211     rcp (new row_graph_impl_type (A_->getGraph (),
00212                                   ExtMatrix_->getGraph (),
00213                                   RowMap_,
00214                                   ColMap_,
00215                                   NumGlobalRows_,
00216                                   NumGlobalRows_, // # global cols == # global rows
00217                                   NumGlobalNonzeros_,
00218                                   MaxNumEntries_,
00219                                   rcp_const_cast<const import_type> (Importer_),
00220                                   rcp_const_cast<const import_type> (ExtImporter_)));
00221   graph_ = rcp_const_cast<const row_graph_type> (rcp_implicit_cast<row_graph_type> (graph));
00222   // Resize temp arrays
00223   Indices_.resize (MaxNumEntries_);
00224   Values_.resize (MaxNumEntries_);
00225 }
00226 
00227 
00228 template<class MatrixType>
00229 OverlappingRowMatrix<MatrixType>::
00230 OverlappingRowMatrix (const Teuchos::RCP<const row_matrix_type>& A,
00231                       const int overlapLevel,
00232                       const int subdomainID)
00233 {
00234   //FIXME
00235   TEUCHOS_TEST_FOR_EXCEPTION(
00236     true, std::logic_error,
00237     "Ifpack2::OverlappingRowMatrix: Subdomain code not implemented yet.");
00238 }
00239 
00240 
00241 template<class MatrixType>
00242 OverlappingRowMatrix<MatrixType>::~OverlappingRowMatrix() {}
00243 
00244 
00245 template<class MatrixType>
00246 Teuchos::RCP<const Teuchos::Comm<int> >
00247 OverlappingRowMatrix<MatrixType>::getComm () const
00248 {
00249   return A_->getComm ();
00250 }
00251 
00252 
00253 template<class MatrixType>
00254 Teuchos::RCP<typename MatrixType::node_type>
00255 OverlappingRowMatrix<MatrixType>::getNode () const
00256 {
00257   return A_->getNode();
00258 }
00259 
00260 
00261 template<class MatrixType>
00262 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00263 OverlappingRowMatrix<MatrixType>::getRowMap () const
00264 {
00265   // FIXME (mfh 12 July 2013) Is this really the right Map to return?
00266   return RowMap_;
00267 }
00268 
00269 
00270 template<class MatrixType>
00271 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00272 OverlappingRowMatrix<MatrixType>::getColMap () const
00273 {
00274   // FIXME (mfh 12 July 2013) Is this really the right Map to return?
00275   return ColMap_;
00276 }
00277 
00278 
00279 template<class MatrixType>
00280 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00281 OverlappingRowMatrix<MatrixType>::getDomainMap () const
00282 {
00283   return A_->getDomainMap();
00284 }
00285 
00286 
00287 template<class MatrixType>
00288 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00289 OverlappingRowMatrix<MatrixType>::getRangeMap () const
00290 {
00291   return A_->getRangeMap ();
00292 }
00293 
00294 
00295 template<class MatrixType>
00296 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00297 OverlappingRowMatrix<MatrixType>::getGraph() const
00298 {
00299   return graph_;
00300 }
00301 
00302 
00303 template<class MatrixType>
00304 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumRows() const
00305 {
00306   return NumGlobalRows_;
00307 }
00308 
00309 
00310 template<class MatrixType>
00311 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumCols() const
00312 {
00313   return NumGlobalRows_;
00314 }
00315 
00316 
00317 template<class MatrixType>
00318 size_t OverlappingRowMatrix<MatrixType>::getNodeNumRows() const
00319 {
00320   return A_->getNodeNumRows () + ExtMatrix_->getNodeNumRows ();
00321 }
00322 
00323 
00324 template<class MatrixType>
00325 size_t OverlappingRowMatrix<MatrixType>::getNodeNumCols() const
00326 {
00327   return this->getNodeNumRows ();
00328 }
00329 
00330 
00331 template<class MatrixType>
00332 typename MatrixType::global_ordinal_type
00333 OverlappingRowMatrix<MatrixType>::getIndexBase () const
00334 {
00335   return A_->getIndexBase();
00336 }
00337 
00338 
00339 template<class MatrixType>
00340 Tpetra::global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumEntries() const
00341 {
00342   return NumGlobalNonzeros_;
00343 }
00344 
00345 
00346 template<class MatrixType>
00347 size_t OverlappingRowMatrix<MatrixType>::getNodeNumEntries() const
00348 {
00349   return A_->getNodeNumEntries () + ExtMatrix_->getNodeNumEntries ();
00350 }
00351 
00352 
00353 template<class MatrixType>
00354 size_t
00355 OverlappingRowMatrix<MatrixType>::
00356 getNumEntriesInGlobalRow (global_ordinal_type globalRow) const
00357 {
00358   const local_ordinal_type localRow = RowMap_->getLocalElement (globalRow);
00359   if (localRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
00360     return Teuchos::OrdinalTraits<size_t>::invalid();
00361   } else {
00362     return getNumEntriesInLocalRow (localRow);
00363   }
00364 }
00365 
00366 
00367 template<class MatrixType>
00368 size_t
00369 OverlappingRowMatrix<MatrixType>::
00370 getNumEntriesInLocalRow (local_ordinal_type localRow) const
00371 {
00372   using Teuchos::as;
00373   const size_t numMyRowsA = A_->getNodeNumRows ();
00374   if (as<size_t> (localRow) < numMyRowsA) {
00375     return A_->getNumEntriesInLocalRow (localRow);
00376   } else {
00377     return ExtMatrix_->getNumEntriesInLocalRow (as<local_ordinal_type> (localRow - numMyRowsA));
00378   }
00379 }
00380 
00381 
00382 template<class MatrixType>
00383 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumDiags() const
00384 {
00385   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalNumDiags() not supported.");
00386 }
00387 
00388 
00389 template<class MatrixType>
00390 size_t OverlappingRowMatrix<MatrixType>::getNodeNumDiags() const
00391 {
00392   return A_->getNodeNumDiags();
00393 }
00394 
00395 
00396 template<class MatrixType>
00397 size_t OverlappingRowMatrix<MatrixType>::getGlobalMaxNumRowEntries() const
00398 {
00399   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalMaxNumRowEntries() not supported.");
00400 }
00401 
00402 
00403 template<class MatrixType>
00404 size_t OverlappingRowMatrix<MatrixType>::getNodeMaxNumRowEntries() const
00405 {
00406   return MaxNumEntries_;
00407 }
00408 
00409 
00410 template<class MatrixType>
00411 bool OverlappingRowMatrix<MatrixType>::hasColMap() const
00412 {
00413   return true;
00414 }
00415 
00416 
00417 template<class MatrixType>
00418 bool OverlappingRowMatrix<MatrixType>::isLowerTriangular() const
00419 {
00420   return A_->isLowerTriangular();
00421 }
00422 
00423 
00424 template<class MatrixType>
00425 bool OverlappingRowMatrix<MatrixType>::isUpperTriangular() const
00426 {
00427   return A_->isUpperTriangular();
00428 }
00429 
00430 
00431 template<class MatrixType>
00432 bool OverlappingRowMatrix<MatrixType>::isLocallyIndexed() const
00433 {
00434   return true;
00435 }
00436 
00437 
00438 template<class MatrixType>
00439 bool OverlappingRowMatrix<MatrixType>::isGloballyIndexed() const
00440 {
00441   return false;
00442 }
00443 
00444 
00445 template<class MatrixType>
00446 bool OverlappingRowMatrix<MatrixType>::isFillComplete() const
00447 {
00448   return true;
00449 }
00450 
00451 
00452 template<class MatrixType>
00453 void
00454 OverlappingRowMatrix<MatrixType>::
00455 getGlobalRowCopy (global_ordinal_type GlobalRow,
00456                   const Teuchos::ArrayView<global_ordinal_type> &Indices,
00457                   const Teuchos::ArrayView<scalar_type>& Values,
00458                   size_t& NumEntries) const
00459 {
00460   const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
00461   if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
00462     NumEntries = Teuchos::OrdinalTraits<size_t>::invalid ();
00463   } else {
00464     if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
00465       A_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
00466     } else {
00467       ExtMatrix_->getGlobalRowCopy (GlobalRow, Indices, Values, NumEntries);
00468     }
00469   }
00470 }
00471 
00472 
00473 template<class MatrixType>
00474 void
00475 OverlappingRowMatrix<MatrixType>::
00476 getLocalRowCopy (local_ordinal_type LocalRow,
00477                  const Teuchos::ArrayView<local_ordinal_type> &Indices,
00478                  const Teuchos::ArrayView<scalar_type> &Values,
00479                  size_t &NumEntries) const
00480 {
00481   using Teuchos::as;
00482   const size_t numMyRowsA = A_->getNodeNumRows ();
00483   if (as<size_t> (LocalRow) < numMyRowsA) {
00484     A_->getLocalRowCopy (LocalRow, Indices, Values, NumEntries);
00485   } else {
00486     ExtMatrix_->getLocalRowCopy (LocalRow - as<local_ordinal_type> (numMyRowsA),
00487                                  Indices, Values, NumEntries);
00488   }
00489 }
00490 
00491 
00492 template<class MatrixType>
00493 void
00494 OverlappingRowMatrix<MatrixType>::
00495 getGlobalRowView (global_ordinal_type GlobalRow,
00496                   Teuchos::ArrayView<const global_ordinal_type>& indices,
00497                   Teuchos::ArrayView<const scalar_type>& values) const
00498 {
00499   const local_ordinal_type LocalRow = RowMap_->getLocalElement (GlobalRow);
00500   if (LocalRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid())  {
00501     indices = Teuchos::null;
00502     values = Teuchos::null;
00503   } else {
00504     if (Teuchos::as<size_t> (LocalRow) < A_->getNodeNumRows ()) {
00505       A_->getGlobalRowView (GlobalRow, indices, values);
00506     } else {
00507       ExtMatrix_->getGlobalRowView (GlobalRow, indices, values);
00508     }
00509   }
00510 }
00511 
00512 
00513 template<class MatrixType>
00514 void
00515 OverlappingRowMatrix<MatrixType>::
00516 getLocalRowView (local_ordinal_type LocalRow,
00517                  Teuchos::ArrayView<const local_ordinal_type>& indices,
00518                  Teuchos::ArrayView<const scalar_type>& values) const
00519 {
00520   using Teuchos::as;
00521   const size_t numMyRowsA = A_->getNodeNumRows ();
00522   if (as<size_t> (LocalRow) < numMyRowsA) {
00523     A_->getLocalRowView (LocalRow, indices, values);
00524   } else {
00525     ExtMatrix_->getLocalRowView (LocalRow - as<local_ordinal_type> (numMyRowsA),
00526                                  indices, values);
00527   }
00528 }
00529 
00530 
00531 template<class MatrixType>
00532 void
00533 OverlappingRowMatrix<MatrixType>::
00534 getLocalDiagCopy (Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& diag) const
00535 {
00536   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getLocalDiagCopy not supported.");
00537 }
00538 
00539 
00540 template<class MatrixType>
00541 void
00542 OverlappingRowMatrix<MatrixType>::
00543 leftScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
00544 {
00545   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
00546 }
00547 
00548 
00549 template<class MatrixType>
00550 void
00551 OverlappingRowMatrix<MatrixType>::
00552 rightScale (const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& x)
00553 {
00554   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
00555 }
00556 
00557 
00558 template<class MatrixType>
00559 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00560 OverlappingRowMatrix<MatrixType>::getFrobeniusNorm () const
00561 {
00562   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
00563 }
00564 
00565 
00566 template<class MatrixType>
00567 void
00568 OverlappingRowMatrix<MatrixType>::
00569 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
00570        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
00571        Teuchos::ETransp mode,
00572        scalar_type alpha,
00573        scalar_type beta) const
00574 {
00575   using Teuchos::ArrayRCP;
00576   using Teuchos::as;
00577   typedef scalar_type RangeScalar;
00578   typedef scalar_type DomainScalar;
00579 
00580   typedef Teuchos::ScalarTraits<RangeScalar> STRS;
00581   TEUCHOS_TEST_FOR_EXCEPTION(
00582     X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00583     "Ifpack2::OverlappingRowMatrix::apply: The input X and the output Y must "
00584     "have the same number of columns.  X.getNumVectors() = "
00585     << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
00586     << ".");
00587 
00588   // FIXME (mfh 13 July 2013) This would be a good candidate for a
00589   // Kokkos local parallel operator implementation.  That would
00590   // obviate the need for getting views of the data and make the code
00591   // below a lot simpler.
00592 
00593   const RangeScalar zero = STRS::zero ();
00594   ArrayRCP<ArrayRCP<const DomainScalar> > x_ptr = X.get2dView();
00595   ArrayRCP<ArrayRCP<RangeScalar> >        y_ptr = Y.get2dViewNonConst();
00596   Y.putScalar(zero);
00597   size_t NumVectors = Y.getNumVectors();
00598 
00599   const size_t numMyRowsA = A_->getNodeNumRows ();
00600   for (size_t i = 0; i < numMyRowsA; ++i) {
00601     size_t Nnz;
00602     // Use this class's getrow to make the below code simpler
00603     A_->getLocalRowCopy (i, Indices_ (),Values_ (), Nnz);
00604     if (mode == Teuchos::NO_TRANS) {
00605       for (size_t j = 0; j < Nnz; ++j)
00606         for (size_t k = 0; k < NumVectors; ++k)
00607           y_ptr[k][i] += as<RangeScalar> (Values_[j]) *
00608             as<RangeScalar> (x_ptr[k][Indices_[j]]);
00609     }
00610     else if (mode == Teuchos::TRANS){
00611       for (size_t j = 0; j < Nnz; ++j)
00612         for (size_t k = 0; k < NumVectors; ++k)
00613           y_ptr[k][Indices_[j]] += as<RangeScalar> (Values_[j]) *
00614             as<RangeScalar> (x_ptr[k][i]);
00615     }
00616     else { // mode == Teuchos::CONJ_TRANS
00617       for (size_t j = 0; j < Nnz; ++j)
00618         for (size_t k = 0; k < NumVectors; ++k)
00619           y_ptr[k][Indices_[j]] +=
00620             STRS::conjugate (as<RangeScalar> (Values_[j])) *
00621             as<RangeScalar> (x_ptr[k][i]);
00622     }
00623   }
00624 
00625   const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
00626   for (size_t i = 0 ; i < numMyRowsB ; ++i) {
00627     size_t Nnz;
00628     // Use this class's getrow to make the below code simpler
00629     ExtMatrix_->getLocalRowCopy (i, Indices_ (), Values_ (), Nnz);
00630     if (mode == Teuchos::NO_TRANS) {
00631       for (size_t j = 0; j < Nnz; ++j)
00632         for (size_t k = 0; k < NumVectors; ++k)
00633           y_ptr[k][numMyRowsA+i] += as<RangeScalar> (Values_[j]) *
00634             as<RangeScalar> (x_ptr[k][Indices_[j]]);
00635     }
00636     else if (mode == Teuchos::TRANS) {
00637       for (size_t j = 0; j < Nnz; ++j)
00638         for (size_t k = 0; k < NumVectors; ++k)
00639           y_ptr[k][numMyRowsA+Indices_[j]] += as<RangeScalar> (Values_[j]) *
00640             as<RangeScalar> (x_ptr[k][i]);
00641     }
00642     else { // mode == Teuchos::CONJ_TRANS
00643       for (size_t j = 0; j < Nnz; ++j)
00644         for (size_t k = 0; k < NumVectors; ++k)
00645           y_ptr[k][numMyRowsA+Indices_[j]] +=
00646             STRS::conjugate (as<RangeScalar> (Values_[j])) *
00647             as<RangeScalar> (x_ptr[k][i]);
00648     }
00649   }
00650 }
00651 
00652 
00653 template<class MatrixType>
00654 void
00655 OverlappingRowMatrix<MatrixType>::
00656 importMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
00657                    Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
00658                    Tpetra::CombineMode CM)
00659 {
00660   OvX.doImport (X, *Importer_, CM);
00661 }
00662 
00663 
00664 template<class MatrixType>
00665 void
00666 OverlappingRowMatrix<MatrixType>::
00667 exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
00668                    Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
00669                    Tpetra::CombineMode CM)
00670 {
00671   X.doExport (OvX, *Importer_, CM);
00672 }
00673 
00674 
00675 template<class MatrixType>
00676 bool OverlappingRowMatrix<MatrixType>::hasTransposeApply () const
00677 {
00678   return true;
00679 }
00680 
00681 
00682 template<class MatrixType>
00683 bool OverlappingRowMatrix<MatrixType>::supportsRowViews () const
00684 {
00685   return false;
00686 }
00687 
00688 template<class MatrixType>
00689 std::string OverlappingRowMatrix<MatrixType>::description() const
00690 {
00691   std::ostringstream oss;
00692   if (isFillComplete()) {
00693     oss << "{ isFillComplete: true"
00694         << ", global rows: " << getGlobalNumRows()
00695         << ", global columns: " << getGlobalNumCols()
00696         << ", global entries: " << getGlobalNumEntries()
00697         << " }";
00698   }
00699   else {
00700     oss << "{ isFillComplete: false"
00701         << ", global rows: " << getGlobalNumRows()
00702         << " }";
00703   }
00704   return oss.str();
00705 }
00706 
00707 template<class MatrixType>
00708 void OverlappingRowMatrix<MatrixType>::describe(Teuchos::FancyOStream &out,
00709             const Teuchos::EVerbosityLevel verbLevel) const
00710 {
00711     using std::endl;
00712     using std::setw;
00713     using Teuchos::as;
00714     using Teuchos::VERB_DEFAULT;
00715     using Teuchos::VERB_NONE;
00716     using Teuchos::VERB_LOW;
00717     using Teuchos::VERB_MEDIUM;
00718     using Teuchos::VERB_HIGH;
00719     using Teuchos::VERB_EXTREME;
00720     using Teuchos::RCP;
00721     using Teuchos::null;
00722     using Teuchos::ArrayView;
00723 
00724     Teuchos::EVerbosityLevel vl = verbLevel;
00725     if (vl == VERB_DEFAULT) {
00726       vl = VERB_LOW;
00727     }
00728     RCP<const Teuchos::Comm<int> > comm = this->getComm();
00729     const int myRank = comm->getRank();
00730     const int numProcs = comm->getSize();
00731     size_t width = 1;
00732     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
00733       ++width;
00734     }
00735     width = std::max<size_t> (width, as<size_t> (11)) + 2;
00736     Teuchos::OSTab tab(out);
00737     //    none: print nothing
00738     //     low: print O(1) info from node 0
00739     //  medium: print O(P) info, num entries per process
00740     //    high: print O(N) info, num entries per row
00741     // extreme: print O(NNZ) info: print indices and values
00742     //
00743     // for medium and higher, print constituent objects at specified verbLevel
00744     if (vl != VERB_NONE) {
00745       if (myRank == 0) {
00746         out << this->description() << std::endl;
00747       }
00748       // O(1) globals, minus what was already printed by description()
00749       //if (isFillComplete() && myRank == 0) {
00750       //  out << "Global number of diagonal entries: " << getGlobalNumDiags() << std::endl;
00751       //  out << "Global max number of entries in a row: " << getGlobalMaxNumRowEntries() << std::endl;
00752       //}
00753       // constituent objects
00754       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
00755         if (myRank == 0) {
00756           out << endl << "Row map:" << endl;
00757         }
00758         getRowMap()->describe(out,vl);
00759         //
00760         if (getColMap() != null) {
00761           if (getColMap() == getRowMap()) {
00762             if (myRank == 0) {
00763               out << endl << "Column map is row map.";
00764             }
00765           }
00766           else {
00767             if (myRank == 0) {
00768               out << endl << "Column map:" << endl;
00769             }
00770             getColMap()->describe(out,vl);
00771           }
00772         }
00773         if (getDomainMap() != null) {
00774           if (getDomainMap() == getRowMap()) {
00775             if (myRank == 0) {
00776               out << endl << "Domain map is row map.";
00777             }
00778           }
00779           else if (getDomainMap() == getColMap()) {
00780             if (myRank == 0) {
00781               out << endl << "Domain map is column map.";
00782             }
00783           }
00784           else {
00785             if (myRank == 0) {
00786               out << endl << "Domain map:" << endl;
00787             }
00788             getDomainMap()->describe(out,vl);
00789           }
00790         }
00791         if (getRangeMap() != null) {
00792           if (getRangeMap() == getDomainMap()) {
00793             if (myRank == 0) {
00794               out << endl << "Range map is domain map." << endl;
00795             }
00796           }
00797           else if (getRangeMap() == getRowMap()) {
00798             if (myRank == 0) {
00799               out << endl << "Range map is row map." << endl;
00800             }
00801           }
00802           else {
00803             if (myRank == 0) {
00804               out << endl << "Range map: " << endl;
00805             }
00806             getRangeMap()->describe(out,vl);
00807           }
00808         }
00809         if (myRank == 0) {
00810           out << endl;
00811         }
00812       }
00813       // O(P) data
00814       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
00815         for (int curRank = 0; curRank < numProcs; ++curRank) {
00816           if (myRank == curRank) {
00817             out << "Process rank: " << curRank << std::endl;
00818             out << "  Number of entries: " << getNodeNumEntries() << std::endl;
00819             if (isFillComplete()) {
00820               out << "  Number of diagonal entries: " << getNodeNumDiags() << std::endl;
00821             }
00822             out << "  Max number of entries per row: " << getNodeMaxNumRowEntries() << std::endl;
00823           }
00824           comm->barrier();
00825           comm->barrier();
00826           comm->barrier();
00827         }
00828       }
00829       // O(N) and O(NNZ) data
00830       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
00831         for (int curRank = 0; curRank < numProcs; ++curRank) {
00832           if (myRank == curRank) {
00833             out << std::setw(width) << "Proc Rank"
00834                 << std::setw(width) << "Global Row"
00835                 << std::setw(width) << "Num Entries";
00836             if (vl == VERB_EXTREME) {
00837               out << std::setw(width) << "(Index,Value)";
00838             }
00839             out << endl;
00840             for (size_t r = 0; r < getNodeNumRows (); ++r) {
00841               const size_t nE = getNumEntriesInLocalRow(r);
00842               typename MatrixType::global_ordinal_type gid = getRowMap()->getGlobalElement(r);
00843               out << std::setw(width) << myRank
00844                   << std::setw(width) << gid
00845                   << std::setw(width) << nE;
00846               if (vl == VERB_EXTREME) {
00847                 if (isGloballyIndexed()) {
00848                   ArrayView<const typename MatrixType::global_ordinal_type> rowinds;
00849                   ArrayView<const typename MatrixType::scalar_type> rowvals;
00850                   getGlobalRowView (gid, rowinds, rowvals);
00851                   for (size_t j = 0; j < nE; ++j) {
00852                     out << " (" << rowinds[j]
00853                         << ", " << rowvals[j]
00854                         << ") ";
00855                   }
00856                 }
00857                 else if (isLocallyIndexed()) {
00858                   ArrayView<const typename MatrixType::local_ordinal_type> rowinds;
00859                   ArrayView<const typename MatrixType::scalar_type> rowvals;
00860                   getLocalRowView (r, rowinds, rowvals);
00861                   for (size_t j=0; j < nE; ++j) {
00862                     out << " (" << getColMap()->getGlobalElement(rowinds[j])
00863                         << ", " << rowvals[j]
00864                         << ") ";
00865                   }
00866                 } // globally or locally indexed
00867               } // vl == VERB_EXTREME
00868               out << endl;
00869             } // for each row r on this process
00870 
00871           } // if (myRank == curRank)
00872           comm->barrier();
00873           comm->barrier();
00874           comm->barrier();
00875         }
00876 
00877         out.setOutputToRootOnly(0);
00878         out << "===========\nlocal matrix\n=================" << std::endl;
00879         out.setOutputToRootOnly(-1);
00880         A_->describe(out,Teuchos::VERB_EXTREME);
00881         out.setOutputToRootOnly(0);
00882         out << "===========\nend of local matrix\n=================" << std::endl;
00883         comm->barrier();
00884         out.setOutputToRootOnly(0);
00885         out << "=================\nghost matrix\n=================" << std::endl;
00886         out.setOutputToRootOnly(-1);
00887         ExtMatrix_->describe(out,Teuchos::VERB_EXTREME);
00888         out.setOutputToRootOnly(0);
00889         out << "===========\nend of ghost matrix\n=================" << std::endl;
00890       }
00891     }
00892 }
00893 
00894 
00895 } // namespace Ifpack2
00896 
00897 #define IFPACK2_OVERLAPPINGROWMATRIX_INSTANT(S,LO,GO,N)                 \
00898   template class Ifpack2::OverlappingRowMatrix< Tpetra::CrsMatrix<S, LO, GO, N> >;
00899 
00900 #endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends