Zoltan2
Zoltan2_AlgRCM.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //   Zoltan2: A package of combinatorial algorithms for scientific computing
00006 //                  Copyright 2012 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Karen Devine      (kddevin@sandia.gov)
00039 //                    Erik Boman        (egboman@sandia.gov)
00040 //                    Siva Rajamanickam (srajama@sandia.gov)
00041 //
00042 // ***********************************************************************
00043 //
00044 // @HEADER
00045 #ifndef _ZOLTAN2_ALGRCM_HPP_
00046 #define _ZOLTAN2_ALGRCM_HPP_
00047 
00048 #include <Zoltan2_Algorithm.hpp>
00049 #include <Zoltan2_GraphModel.hpp>
00050 #include <Zoltan2_OrderingSolution.hpp>
00051 #include <Zoltan2_Sort.hpp>
00052 #include <queue>
00053 
00054 
00058 
00059 
00060 namespace Zoltan2{
00061 
00062 template <typename Adapter>
00063 class AlgRCM : public Algorithm<Adapter>
00064 {
00065   private:
00066 
00067   const RCP<GraphModel<Adapter> > model;
00068   const RCP<Teuchos::ParameterList> &pl;
00069   const RCP<Teuchos::Comm<int> > &comm;
00070 
00071   public:
00072 
00073   typedef typename Adapter::lno_t lno_t;
00074   typedef typename Adapter::gid_t gid_t;
00075   typedef typename Adapter::scalar_t scalar_t;
00076 
00077   AlgRCM(
00078     const RCP<GraphModel<Adapter> > &model__,
00079     const RCP<Teuchos::ParameterList> &pl__,
00080     const RCP<Teuchos::Comm<int> > &comm__
00081   ) : model(model__), pl(pl__), comm(comm__)
00082   {
00083   }
00084 
00085   int order(const RCP<OrderingSolution<gid_t, lno_t> > &solution)
00086   {
00087     int ierr= 0;
00088 
00089     HELLO;
00090   
00091     // Get local graph.
00092     ArrayView<const lno_t> edgeIds;
00093     ArrayView<const lno_t> offsets;
00094     ArrayView<StridedData<lno_t, scalar_t> > wgts;
00095   
00096     const size_t nVtx = model->getLocalNumVertices();
00097     model->getLocalEdgeList(edgeIds, offsets, wgts); 
00098     const int numWeightsPerEdge = model->getNumWeightsPerEdge();
00099     if (numWeightsPerEdge > 1){
00100       throw std::runtime_error("Multiple weights not supported.");
00101     }
00102   
00103 #if 0
00104     // Debug
00105     cout << "Debug: Local graph from getLocalEdgeList" << endl;
00106     cout << "rank " << comm->getRank() << ": nVtx= " << nVtx << endl;
00107     cout << "rank " << comm->getRank() << ": edgeIds: " << edgeIds << endl;
00108     cout << "rank " << comm->getRank() << ": offsets: " << offsets << endl;
00109 #endif
00110   
00111     // RCM constructs invPerm, not perm
00112     ArrayRCP<lno_t> invPerm = solution->getPermutationRCP(true);
00113   
00114     // Check if there are actually edges to reorder.
00115     // If there are not, then just use the natural ordering.
00116     if (offsets[nVtx] == 0) {
00117       for (size_t i = 0; i < nVtx; ++i) {
00118         invPerm[i] = i;
00119       }
00120       solution->setHaveInverse(true);
00121       return 0;
00122     }
00123   
00124     // Set the label of each vertex to invalid.
00125     Tpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
00126     for (size_t i = 0; i < nVtx; ++i) {
00127       invPerm[i] = INVALID;
00128     }
00129   
00130     // Loop over all connected components.
00131     // Do BFS within each component.
00132     lno_t root = 0;
00133     std::queue<lno_t> Q;
00134     size_t count = 0; // CM label, reversed later
00135     size_t next = 0;  // next unmarked vertex
00136     Teuchos::Array<std::pair<lno_t, size_t> >  children; // children and their degrees
00137   
00138     while (count < nVtx) {
00139   
00140       // Find suitable root vertex for this component.
00141       // First find an unmarked vertex, use to find root in next component.
00142       while ((next < nVtx) && (static_cast<Tpetra::global_size_t>(invPerm[next]) != INVALID)) next++;
00143 
00144       // Select root method. Pseudoperipheral usually gives the best
00145       // ordering, but the user may choose a faster method.
00146       std::string root_method = pl->get("root_method", "pseudoperipheral");
00147       if (root_method == std::string("first"))
00148         root = next;
00149       else if (root_method == std::string("smallest_degree"))
00150         root = findSmallestDegree(next, nVtx, edgeIds, offsets);
00151       else if (root_method == std::string("pseudoperipheral"))
00152         root = findPseudoPeripheral(next, nVtx, edgeIds, offsets);
00153       else {
00154         // This should never happen if pl was validated.
00155         throw std::runtime_error("invalid root_method");
00156       }
00157 
00158       // Label connected component starting at root
00159       Q.push(root);
00160       //cout << "Debug: invPerm[" << root << "] = " << count << endl;
00161       invPerm[root] = count++;
00162   
00163       while (Q.size()){
00164         // Get a vertex from the queue
00165         lno_t v = Q.front();
00166         Q.pop();
00167         //cout << "Debug: v= " << v << ", offsets[v] = " << offsets[v] << endl;
00168   
00169         // Add unmarked children to list of pairs, to be added to queue.
00170         children.resize(0);
00171         for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00172           lno_t child = edgeIds[ptr];
00173           if (static_cast<Tpetra::global_size_t>(invPerm[child]) == INVALID){
00174             // Not visited yet; add child to list of pairs.
00175             std::pair<lno_t,size_t> newchild;
00176             newchild.first = child;
00177             newchild.second = offsets[child+1] - offsets[child];
00178             children.push_back(newchild); 
00179           }
00180         }
00181         // Sort children by increasing degree
00182         // TODO: If edge weights, sort children by decreasing weight,
00183         SortPairs<lno_t,size_t> zort;
00184         zort.sort(children);
00185 
00186         typename Teuchos::Array<std::pair<lno_t,size_t> >::iterator it = children.begin ();
00187         for ( ; it != children.end(); ++it){
00188           // Push children on the queue in sorted order.
00189           lno_t child = it->first;
00190           invPerm[child] = count++; // Label as we push on Q
00191           Q.push(child);
00192           //cout << "Debug: invPerm[" << child << "] = " << count << endl;
00193         }
00194       }
00195     }
00196   
00197     // Reverse labels for RCM
00198     bool reverse = true; // TODO: Make parameter
00199     if (reverse) {
00200       lno_t temp;
00201       for (size_t i=0; i < nVtx/2; ++i) {
00202         // Swap (invPerm[i], invPerm[nVtx-i])
00203         temp = invPerm[i];
00204         invPerm[i] = invPerm[nVtx-1-i];
00205         invPerm[nVtx-1-i] = temp;
00206       }
00207     }
00208   
00209     solution->setHaveInverse(true);
00210     return ierr;
00211   }
00212 
00213   private:
00214   // Find a smallest degree vertex in component containing v
00215   lno_t findSmallestDegree(
00216     lno_t v,
00217     lno_t nVtx,
00218     ArrayView<const lno_t> edgeIds,
00219     ArrayView<const lno_t> offsets)
00220   {
00221     std::queue<lno_t> Q;
00222     Teuchos::Array<bool> mark(nVtx);
00223 
00224     // Do BFS and compute smallest degree as we go
00225     lno_t smallestDegree = nVtx;
00226     lno_t smallestVertex = 0;
00227 
00228     // Clear mark array - nothing marked yet
00229     for (int i=0; i<nVtx; i++)
00230       mark[i] = false;
00231 
00232     // Start from v
00233     Q.push(v);
00234     while (Q.size()){
00235       // Get first vertex from the queue
00236       v = Q.front();
00237       Q.pop();
00238       // Check degree of v
00239       lno_t deg = offsets[v+1] - offsets[v];
00240       if (deg < smallestDegree){
00241         smallestDegree = deg;
00242         smallestVertex = v;
00243       }
00244       // Add unmarked children to queue
00245       for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00246         lno_t child = edgeIds[ptr];
00247         if (!mark[child]){
00248           mark[child] = true; 
00249           Q.push(child);
00250         }
00251       }
00252     }
00253     return smallestVertex;
00254   }
00255 
00256   // Find a pseudoperipheral vertex in component containing v
00257   lno_t findPseudoPeripheral(
00258     lno_t v,
00259     lno_t nVtx,
00260     ArrayView<const lno_t> edgeIds,
00261     ArrayView<const lno_t> offsets)
00262   {
00263     std::queue<lno_t> Q;
00264     Teuchos::Array<bool> mark(nVtx);
00265 
00266     // Do BFS a couple times, pick vertex last visited (furthest away)
00267     const int numBFS = 2;
00268     for (int bfs=0; bfs<numBFS; bfs++){
00269       // Clear mark array - nothing marked yet
00270       for (int i=0; i<nVtx; i++)
00271         mark[i] = false;
00272       // Start from v
00273       Q.push(v);
00274       while (Q.size()){
00275         // Get first vertex from the queue
00276         v = Q.front();
00277         Q.pop();
00278         // Add unmarked children to queue
00279         for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00280           lno_t child = edgeIds[ptr];
00281           if (!mark[child]){
00282             mark[child] = true; 
00283             Q.push(child);
00284           }
00285         }
00286       }
00287     }
00288     return v;
00289   }
00290   
00291 };
00292 }
00293 #endif