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_GraphModel.hpp>
00049 #include <Zoltan2_OrderingSolution.hpp>
00050 #include <Zoltan2_Sort.hpp>
00051 #include <queue>
00052 
00053 
00057 
00058 
00059 namespace Zoltan2{
00060 
00061 template <typename Adapter>
00062 class AlgRCM
00063 {
00064   private:
00065     typedef typename Adapter::lno_t lno_t;
00066     typedef typename Adapter::gno_t gno_t;
00067     typedef typename Adapter::scalar_t scalar_t;
00068   
00069   public:
00070   AlgRCM()
00071   {
00072   }
00073 
00074   int order(
00075     const RCP<GraphModel<Adapter> > &model,
00076     const RCP<OrderingSolution<typename Adapter::gid_t,
00077              typename Adapter::lno_t> > &solution,
00078     const RCP<Teuchos::ParameterList> &pl,
00079     const RCP<Teuchos::Comm<int> > &comm
00080   )
00081   {
00082     int ierr= 0;
00083 
00084     HELLO;
00085   
00086     // Get local graph.
00087     ArrayView<const lno_t> edgeIds;
00088     ArrayView<const lno_t> offsets;
00089     ArrayView<StridedData<lno_t, scalar_t> > wgts;
00090   
00091     const size_t nVtx = model->getLocalNumVertices();
00092     model->getLocalEdgeList(edgeIds, offsets, wgts); 
00093     const int numWeightsPerEdge = model->getNumWeightsPerEdge();
00094     if (numWeightsPerEdge > 1){
00095       throw std::runtime_error("Multiple weights not supported.");
00096     }
00097   
00098 #if 0
00099     // Debug
00100     cout << "Debug: Local graph from getLocalEdgeList" << endl;
00101     cout << "rank " << comm->getRank() << ": nVtx= " << nVtx << endl;
00102     cout << "rank " << comm->getRank() << ": edgeIds: " << edgeIds << endl;
00103     cout << "rank " << comm->getRank() << ": offsets: " << offsets << endl;
00104 #endif
00105   
00106     // RCM constructs invPerm, not perm
00107     ArrayRCP<lno_t> invPerm = solution->getPermutationRCP(true);
00108   
00109     // Check if there are actually edges to reorder.
00110     // If there are not, then just use the natural ordering.
00111     if (offsets[nVtx] == 0) {
00112       for (size_t i = 0; i < nVtx; ++i) {
00113         invPerm[i] = i;
00114       }
00115       solution->setHaveInverse(true);
00116       return 0;
00117     }
00118   
00119     // Set the label of each vertex to invalid.
00120     Tpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
00121     for (size_t i = 0; i < nVtx; ++i) {
00122       invPerm[i] = INVALID;
00123     }
00124   
00125     // Loop over all connected components.
00126     // Do BFS within each component.
00127     lno_t root = 0;
00128     std::queue<lno_t> Q;
00129     size_t count = 0; // CM label, reversed later
00130     size_t next = 0;  // next unmarked vertex
00131     Teuchos::Array<std::pair<lno_t, size_t> >  children; // children and their degrees
00132   
00133     while (count < nVtx) {
00134   
00135       // Find suitable root vertex for this component.
00136       // First find an unmarked vertex, use to find root in next component.
00137       while ((next < nVtx) && (static_cast<Tpetra::global_size_t>(invPerm[next]) != INVALID)) next++;
00138 
00139       // Select root method. Pseudoperipheral usually gives the best
00140       // ordering, but the user may choose a faster method.
00141       std::string root_method = pl->get("root_method", "pseudoperipheral");
00142       if (root_method == std::string("first"))
00143         root = next;
00144       else if (root_method == std::string("smallest_degree"))
00145         root = findSmallestDegree(next, nVtx, edgeIds, offsets);
00146       else if (root_method == std::string("pseudoperipheral"))
00147         root = findPseudoPeripheral(next, nVtx, edgeIds, offsets);
00148       else {
00149         // This should never happen if pl was validated.
00150         throw std::runtime_error("invalid root_method");
00151       }
00152 
00153       // Label connected component starting at root
00154       Q.push(root);
00155       //cout << "Debug: invPerm[" << root << "] = " << count << endl;
00156       invPerm[root] = count++;
00157   
00158       while (Q.size()){
00159         // Get a vertex from the queue
00160         lno_t v = Q.front();
00161         Q.pop();
00162         //cout << "Debug: v= " << v << ", offsets[v] = " << offsets[v] << endl;
00163   
00164         // Add unmarked children to list of pairs, to be added to queue.
00165         children.resize(0);
00166         for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00167           lno_t child = edgeIds[ptr];
00168           if (static_cast<Tpetra::global_size_t>(invPerm[child]) == INVALID){
00169             // Not visited yet; add child to list of pairs.
00170             std::pair<lno_t,size_t> newchild;
00171             newchild.first = child;
00172             newchild.second = offsets[child+1] - offsets[child];
00173             children.push_back(newchild); 
00174           }
00175         }
00176         // Sort children by increasing degree
00177         // TODO: If edge weights, sort children by decreasing weight,
00178         SortPairs<lno_t,size_t> zort;
00179         zort.sort(children);
00180 
00181         typename Teuchos::Array<std::pair<lno_t,size_t> >::iterator it = children.begin ();
00182         for ( ; it != children.end(); ++it){
00183           // Push children on the queue in sorted order.
00184           lno_t child = it->first;
00185           invPerm[child] = count++; // Label as we push on Q
00186           Q.push(child);
00187           //cout << "Debug: invPerm[" << child << "] = " << count << endl;
00188         }
00189       }
00190     }
00191   
00192     // Reverse labels for RCM
00193     bool reverse = true; // TODO: Make parameter
00194     if (reverse) {
00195       lno_t temp;
00196       for (size_t i=0; i < nVtx/2; ++i) {
00197         // Swap (invPerm[i], invPerm[nVtx-i])
00198         temp = invPerm[i];
00199         invPerm[i] = invPerm[nVtx-1-i];
00200         invPerm[nVtx-1-i] = temp;
00201       }
00202     }
00203   
00204     solution->setHaveInverse(true);
00205     return ierr;
00206   }
00207 
00208   private:
00209   // Find a smallest degree vertex in component containing v
00210   lno_t findSmallestDegree(
00211     lno_t v,
00212     lno_t nVtx,
00213     ArrayView<const lno_t> edgeIds,
00214     ArrayView<const lno_t> offsets)
00215   {
00216     std::queue<lno_t> Q;
00217     Teuchos::Array<bool> mark(nVtx);
00218 
00219     // Do BFS and compute smallest degree as we go
00220     lno_t smallestDegree = nVtx;
00221     lno_t smallestVertex = 0;
00222 
00223     // Clear mark array - nothing marked yet
00224     for (int i=0; i<nVtx; i++)
00225       mark[i] = false;
00226 
00227     // Start from v
00228     Q.push(v);
00229     while (Q.size()){
00230       // Get first vertex from the queue
00231       v = Q.front();
00232       Q.pop();
00233       // Check degree of v
00234       lno_t deg = offsets[v+1] - offsets[v];
00235       if (deg < smallestDegree){
00236         smallestDegree = deg;
00237         smallestVertex = v;
00238       }
00239       // Add unmarked children to queue
00240       for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00241         lno_t child = edgeIds[ptr];
00242         if (!mark[child]){
00243           mark[child] = true; 
00244           Q.push(child);
00245         }
00246       }
00247     }
00248     return smallestVertex;
00249   }
00250 
00251   // Find a pseudoperipheral vertex in component containing v
00252   lno_t findPseudoPeripheral(
00253     lno_t v,
00254     lno_t nVtx,
00255     ArrayView<const lno_t> edgeIds,
00256     ArrayView<const lno_t> offsets)
00257   {
00258     std::queue<lno_t> Q;
00259     Teuchos::Array<bool> mark(nVtx);
00260 
00261     // Do BFS a couple times, pick vertex last visited (furthest away)
00262     const int numBFS = 2;
00263     for (int bfs=0; bfs<numBFS; bfs++){
00264       // Clear mark array - nothing marked yet
00265       for (int i=0; i<nVtx; i++)
00266         mark[i] = false;
00267       // Start from v
00268       Q.push(v);
00269       while (Q.size()){
00270         // Get first vertex from the queue
00271         v = Q.front();
00272         Q.pop();
00273         // Add unmarked children to queue
00274         for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
00275           lno_t child = edgeIds[ptr];
00276           if (!mark[child]){
00277             mark[child] = true; 
00278             Q.push(child);
00279           }
00280         }
00281       }
00282     }
00283     return v;
00284   }
00285   
00286 };
00287 }
00288 #endif