Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_OverlappingPartitioner_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_OVERLAPPINGPARTITIONER_DEF_HPP
00044 #define IFPACK2_OVERLAPPINGPARTITIONER_DEF_HPP
00045 #include "Ifpack2_ConfigDefs.hpp"
00046 #include "Ifpack2_OverlappingPartitioner_decl.hpp"
00047 #include "Teuchos_Array.hpp"
00048 #include "Teuchos_ArrayRCP.hpp"
00049 #include <vector>
00050 #include <string>
00051 
00052 namespace Ifpack2 {
00053 
00054 template<class GraphType>
00055 OverlappingPartitioner<GraphType>::
00056 OverlappingPartitioner (const Teuchos::RCP<const row_graph_type>& graph) :
00057   NumLocalParts_ (1),
00058   Graph_ (graph),
00059   OverlappingLevel_ (0),
00060   IsComputed_ (false),
00061   verbose_ (false)
00062 {}
00063 
00064 
00065 template<class GraphType>
00066 OverlappingPartitioner<GraphType>::~OverlappingPartitioner() {}
00067 
00068 
00069 template<class GraphType>
00070 int
00071 OverlappingPartitioner<GraphType>::numLocalParts () const 
00072 {
00073   return NumLocalParts_;
00074 }
00075 
00076 
00077 template<class GraphType>
00078 int OverlappingPartitioner<GraphType>::overlappingLevel () const
00079 {
00080   return OverlappingLevel_;
00081 }
00082 
00083 
00084 template<class GraphType>
00085 typename GraphType::local_ordinal_type 
00086 OverlappingPartitioner<GraphType>::
00087 operator () (const local_ordinal_type MyRow) const
00088 {
00089   TEUCHOS_TEST_FOR_EXCEPTION(
00090     MyRow < 0 || Teuchos::as<size_t> (MyRow) > Graph_->getNodeNumRows (), 
00091     std::runtime_error, 
00092     "Ifpack2::OverlappingPartitioner::operator(): "
00093     "Invalid local row index " << MyRow << ".");
00094    
00095   return Partition_[MyRow];
00096 }
00097 
00098 
00099 //==============================================================================
00100 template<class GraphType>
00101 typename GraphType::local_ordinal_type 
00102 OverlappingPartitioner<GraphType>::
00103 operator() (const local_ordinal_type i, const local_ordinal_type j) const
00104 {
00105   TEUCHOS_TEST_FOR_EXCEPTION(
00106     i < 0 || i > Teuchos::as<local_ordinal_type> (NumLocalParts_), 
00107     std::runtime_error, 
00108     "Ifpack2::OverlappingPartitioner::operator(): "
00109     "Invalid local row index i=" << i << ".");
00110   TEUCHOS_TEST_FOR_EXCEPTION(
00111     j < 0 || j > Teuchos::as<local_ordinal_type> (Parts_[i].size ()),
00112     std::runtime_error, 
00113     "Ifpack2::OverlappingPartitioner::operator(): "
00114     "Invalid node index j=" << j << ".");
00115   return Parts_[i][j];
00116 }
00117 
00118 //==============================================================================
00119 template<class GraphType>
00120 size_t
00121 OverlappingPartitioner<GraphType>::
00122 numRowsInPart (const local_ordinal_type Part) const
00123 {
00124   TEUCHOS_TEST_FOR_EXCEPTION(
00125     Part < 0 || Part > Teuchos::as<local_ordinal_type> (NumLocalParts_), 
00126     std::runtime_error, 
00127     "Ifpack2::OverlappingPartitioner::numRowsInPart: "
00128     "Invalid partition index Part=" << Part << ".");
00129   return Parts_[Part].size ();
00130 }
00131 
00132 //==============================================================================
00133 template<class GraphType>
00134 void
00135 OverlappingPartitioner<GraphType>::
00136 rowsInPart (const local_ordinal_type Part,
00137       Teuchos::ArrayRCP<local_ordinal_type>& List) const
00138 {
00139   // Let numRowsInPart do the sanity checking...
00140   const size_t numRows = numRowsInPart (Part); 
00141   for (size_t i = 0; i < numRows; ++i) {
00142     List[i] = Parts_[Part][i];
00143   }
00144 }
00145 
00146 //==============================================================================
00147 template<class GraphType>
00148 Teuchos::ArrayView<const typename GraphType::local_ordinal_type>
00149 OverlappingPartitioner<GraphType>::nonOverlappingPartition () const
00150 {
00151   return Partition_.view (0, Graph_->getNodeNumRows ());
00152 }
00153 
00154 //==============================================================================
00155 template<class GraphType>
00156 void
00157 OverlappingPartitioner<GraphType>::
00158 setParameters (Teuchos::ParameterList& List)
00159 {
00160   NumLocalParts_    = List.get("partitioner: local parts", NumLocalParts_);
00161   OverlappingLevel_ = List.get("partitioner: overlap", OverlappingLevel_);
00162   verbose_          = List.get("partitioner: print level", verbose_);
00163 
00164   if (NumLocalParts_ < 0) {
00165     NumLocalParts_ = Graph_->getNodeNumRows() / (-NumLocalParts_);
00166   }
00167   if (NumLocalParts_ == 0) {
00168     NumLocalParts_ = 1;
00169   }
00170   
00171   // Sanity checking
00172   TEUCHOS_TEST_FOR_EXCEPTION(
00173     NumLocalParts_ < 0 || 
00174     Teuchos::as<size_t> (NumLocalParts_) > Graph_->getNodeNumRows(),
00175     std::runtime_error, 
00176     "Ifpack2::OverlappingPartitioner::setParameters: "
00177     "Invalid NumLocalParts_ = " << NumLocalParts_ << ".");
00178   TEUCHOS_TEST_FOR_EXCEPTION(
00179     OverlappingLevel_ < 0, std::runtime_error,
00180     "Ifpack2::OverlappingPartitioner::setParameters: "
00181     "Invalid OverlappingLevel_ = " << OverlappingLevel_ << ".");
00182 
00183   setPartitionParameters(List);
00184 }
00185 
00186 //==============================================================================
00187 template<class GraphType>
00188 void OverlappingPartitioner<GraphType>::compute()
00189 {
00190   using std::cout;
00191   using std::endl;
00192 
00193   TEUCHOS_TEST_FOR_EXCEPTION(
00194     NumLocalParts_ < 1 || OverlappingLevel_ < 0,
00195     std::runtime_error, 
00196     "Ifpack2::OverlappingPartitioner::compute: "
00197     "Invalid NumLocalParts_ or OverlappingLevel_.");
00198 
00199   // std::string's constructor has some overhead, so it's better to
00200   // use const char[] for local constant strings.
00201   const char printMsg[] = "OverlappingPartitioner: ";
00202 
00203   if (verbose_ && (Graph_->getComm()->getRank() == 0)) {
00204     cout << printMsg << "Number of local parts          = " 
00205    << NumLocalParts_ << endl;
00206     cout << printMsg << "Approx. Number of global parts = " 
00207    << NumLocalParts_ * Graph_->getComm ()->getSize () << endl;
00208     cout << printMsg << "Amount of overlap              = " 
00209    << OverlappingLevel_ << endl;
00210   }
00211 
00212   // 1.- allocate memory 
00213   Partition_.resize (Graph_->getNodeNumRows ());
00214   Parts_.resize (NumLocalParts_);
00215 
00216   // 2.- sanity checks on input graph
00217   TEUCHOS_TEST_FOR_EXCEPTION( 
00218     ! Graph_->isFillComplete (), std::runtime_error, 
00219     "Ifpack2::OverlappingPartitioner::compute: "
00220     "The input graph must be fill complete.");
00221 
00222   TEUCHOS_TEST_FOR_EXCEPTION( 
00223     Graph_->getGlobalNumRows () != Graph_->getGlobalNumCols (),
00224     std::runtime_error, 
00225     "Ifpack2::OverlappingPartitioner::compute: "
00226     "The input graph must be (globally) square.");
00227  
00228   // 3.- perform non-overlapping partition 
00229   computePartitions ();
00230 
00231   // 4.- compute the partitions with overlapping
00232   computeOverlappingPartitions ();
00233 
00234   // 5.- mark as computed
00235   IsComputed_ = true;
00236 }
00237 
00238 //==============================================================================
00239 template<class GraphType>
00240 void OverlappingPartitioner<GraphType>::computeOverlappingPartitions()
00241 {
00242   const local_ordinal_type invalid = 
00243     Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
00244 
00245   // Old FIXME from Ifpack: the first part of this function should be elsewhere
00246  
00247   // start defining the subgraphs for no overlap
00248 
00249   std::vector<size_t> sizes;
00250   sizes.resize (NumLocalParts_);
00251 
00252   // 1.- compute how many rows are in each subgraph
00253   for (int i = 0; i < NumLocalParts_; ++i) {
00254     sizes[i] = 0;
00255   }
00256 
00257   for (size_t i = 0; i < Graph_->getNodeNumRows (); ++i) {
00258     TEUCHOS_TEST_FOR_EXCEPTION(
00259       Partition_[i] >= NumLocalParts_, std::runtime_error, 
00260       "Ifpack2::OverlappingPartitioner::computeOverlappingPartitions: "
00261       "Partition_[i] > NumLocalParts_.");
00262     // invalid indicates that this unknown is not in a nonoverlapping
00263     // partition
00264     if (Partition_[i] != invalid) {
00265       sizes[Partition_[i]]++;
00266     }
00267   }
00268 
00269   // 2.- allocate space for each subgraph
00270   for (int i = 0; i < NumLocalParts_; ++i) {
00271     Parts_[i].resize (sizes[i]);
00272   }
00273 
00274   // 3.- cycle over all rows and populate the vectors
00275   for (int i = 0; i < NumLocalParts_; ++i) {
00276     sizes[i] = 0;
00277   }
00278 
00279   for (size_t i = 0; i < Graph_->getNodeNumRows (); ++i) {
00280     const local_ordinal_type part = Partition_[i];
00281     if (part != invalid) {
00282       const size_t count = sizes[part];
00283       Parts_[part][count] = i;
00284       sizes[part]++;
00285     }
00286   }
00287 
00288   // If there is no overlap, we're done, so return
00289   if (OverlappingLevel_ == 0) {
00290     return;
00291   }
00292 
00293   // wider overlap requires further computations
00294   for (int level = 1; level <= OverlappingLevel_; ++level) {
00295     std::vector<std::vector<size_t> > tmp;
00296     tmp.resize (NumLocalParts_);
00297 
00298     // cycle over all rows in the local graph (that is the overlapping
00299     // graph). For each row, all columns will belong to the subgraph
00300     // of row `i'.
00301 
00302     int MaxNumEntries_tmp = Graph_->getNodeMaxNumRowEntries();
00303     Teuchos::Array<local_ordinal_type> Indices;
00304     Indices.resize (MaxNumEntries_tmp);
00305 
00306     for (int part = 0; part < NumLocalParts_ ; ++part) {
00307       for (size_t i = 0; i < Teuchos::as<size_t> (Parts_[part].size ()); ++i) {  
00308   const local_ordinal_type LRID = Parts_[part][i];
00309         
00310   size_t NumIndices;
00311   Graph_->getLocalRowCopy (LRID, Indices (), NumIndices);
00312 
00313   for (size_t j = 0; j < NumIndices; ++j) {
00314     // use *local* indices only
00315     const local_ordinal_type col = Indices[j];
00316           if (Teuchos::as<size_t> (col) >= Graph_->getNodeNumRows ()) {
00317             continue;
00318     }
00319 
00320     // has this column already been inserted?
00321     std::vector<size_t>::iterator where = 
00322       std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (col));
00323 
00324     if (where == tmp[part].end()) {
00325       tmp[part].push_back (col);
00326     }
00327   }
00328 
00329         // has this column already been inserted?
00330   std::vector<size_t>::iterator where =
00331           std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (LRID));
00332         
00333         // This happens here b/c Vanka on Stokes with Stabilized elements will have
00334         // a zero pivot entry if this gets pushed back first. So... Last.
00335         if (where == tmp[part].end ()) {
00336           tmp[part].push_back (LRID);
00337         }
00338       }
00339     }
00340 
00341     // now I convert the STL vectors into Teuchos Array RCP's
00342     //
00343     // FIXME (mfh 12 July 2013) You could have started with ArrayRCP
00344     // in the first place (which implements push_back and iterators)
00345     // and avoided the copy.
00346     for (int i = 0; i < NumLocalParts_; ++i) {
00347       Parts_[i].resize (tmp[i].size ());
00348       for (size_t j = 0; j < tmp[i].size (); ++j) {
00349   Parts_[i][j] = tmp[i][j];
00350       }
00351     }
00352   }
00353 }
00354 
00355 //==============================================================================
00356 template<class GraphType>
00357 bool OverlappingPartitioner<GraphType>::isComputed() const
00358 {
00359   return IsComputed_;
00360 }
00361 
00362 //==============================================================================
00363 template<class GraphType>
00364 std::ostream& 
00365 OverlappingPartitioner<GraphType>::print (std::ostream& os) const
00366 {
00367   Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
00368   fos.setOutputToRootOnly (0);
00369   describe (fos);
00370   return os;
00371 }
00372 
00373 //==============================================================================
00374 template<class GraphType>
00375 std::string OverlappingPartitioner<GraphType>::description() const
00376 {
00377   std::ostringstream oss;
00378   oss << Teuchos::Describable::description();
00379   if (isComputed()) {
00380     oss << "{status = computed";
00381   }
00382   else {
00383     oss << "{status = is not computed";
00384   }
00385   oss <<"}";
00386   return oss.str();
00387 }
00388 
00389 //==============================================================================
00390 template<class GraphType>
00391 void  OverlappingPartitioner<GraphType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
00392 {
00393   using std::endl;
00394   if (verbLevel == Teuchos::VERB_NONE) {
00395     return;
00396   }
00397 
00398   os << "================================================================================" << endl;
00399   os << "Ifpack2::OverlappingPartitioner" << endl;
00400   os << "Number of local rows  = " << Graph_->getNodeNumRows() << endl;
00401   os << "Number of global rows = " << Graph_->getGlobalNumRows() << endl;
00402   os << "Number of local parts = " << NumLocalParts_ << endl;
00403   os << "Overlapping level     = " << OverlappingLevel_ << endl;
00404   os << "Is computed           = " << IsComputed_ << endl;
00405   os << "================================================================================" << endl;
00406 }
00407 
00408 
00409 }// namespace Ifpack2
00410 
00411 #define IFPACK2_OVERLAPPINGPARTITIONER_INSTANT(LO,GO,N) \
00412   template class Ifpack2::OverlappingPartitioner<Tpetra::CrsGraph< LO, GO, N > >; \
00413   template class Ifpack2::OverlappingPartitioner<Tpetra::RowGraph< LO, GO, N > >;
00414 
00415 #endif // IFPACK2_OVERLAPPINGPARTITIONER_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends