Zoltan2
Zoltan2_GraphModel.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 
00050 #ifndef _ZOLTAN2_GRAPHMODEL_HPP_
00051 #define _ZOLTAN2_GRAPHMODEL_HPP_
00052 
00053 #include <Zoltan2_Model.hpp>
00054 #include <Zoltan2_InputTraits.hpp>
00055 #include <Zoltan2_MatrixAdapter.hpp>
00056 #include <Zoltan2_GraphAdapter.hpp>
00057 #include <Zoltan2_IdentifierAdapter.hpp>
00058 #include <Zoltan2_VectorAdapter.hpp>
00059 #include <Zoltan2_StridedData.hpp>
00060 
00061 #include <vector>
00062 #include <Teuchos_Hashtable.hpp>
00063 
00064 namespace Zoltan2 {
00065 
00066 
00068 
00104 template <typename User>
00105 size_t removeUndesiredEdges(
00106   const RCP<const Environment> &env,
00107   int myRank,
00108   bool removeSelfEdges,
00109   bool removeOffProcessEdges,
00110   bool removeOffGroupEdges,
00111   ArrayView<const typename InputTraits<User>::gid_t> &gids,
00112   ArrayView<const typename InputTraits<User>::gid_t> &gidNbors,
00113   ArrayView<const int> &procIds,
00114   ArrayView<StridedData<typename InputTraits<User>::lno_t,
00115     typename InputTraits<User>::scalar_t> > &edgeWeights,
00116   ArrayView<const typename InputTraits<User>::lno_t> &offsets,
00117   ArrayRCP<const typename InputTraits<User>::gid_t> &newGidNbors, // out
00118   typename InputTraits<User>::scalar_t **&newWeights,       // out
00119   ArrayRCP<const typename InputTraits<User>::lno_t> &newOffsets)  // out
00120 {
00121   typedef typename InputTraits<User>::gid_t gid_t;
00122   typedef typename InputTraits<User>::scalar_t scalar_t;
00123   typedef typename InputTraits<User>::lno_t lno_t;
00124   size_t numKeep = 0;
00125 
00126   size_t numVtx = offsets.size() - 1;
00127   size_t numNbors = gidNbors.size();
00128 
00129   env->localInputAssertion(__FILE__, __LINE__, "need more input",
00130     (!removeSelfEdges ||
00131       gids.size() >=
00132        static_cast<typename ArrayView<const gid_t>::size_type>(numVtx))
00133       &&
00134     (!removeOffProcessEdges ||
00135       procIds.size() >=
00136        static_cast<typename ArrayView<const int>::size_type>(numNbors)) &&
00137     (!removeOffGroupEdges ||
00138       procIds.size() >=
00139        static_cast<typename ArrayView<const int>::size_type>(numNbors)),
00140     BASIC_ASSERTION);
00141 
00142   // initialize edge weight array
00143 
00144   newWeights = NULL;
00145   int eDim = edgeWeights.size();
00146 
00147   // count desired edges
00148 
00149   lno_t *offs = new lno_t [numVtx + 1];
00150   env->localMemoryAssertion(__FILE__, __LINE__, numVtx+1, offs);
00151   for (size_t i = 0; i < numVtx+1; i++) offs[i] = 0;
00152   ArrayRCP<const lno_t> offArray = arcp(offs, 0, numVtx+1, true);
00153 
00154   const lno_t *allOffs = offsets.getRawPtr();
00155   const gid_t *allIds = gidNbors.getRawPtr();
00156 
00157   const gid_t *vtx = NULL;
00158   const int *proc = NULL;
00159 
00160   if (gids.size() > 0)
00161     vtx = gids.getRawPtr();
00162 
00163   if (procIds.size() > 0)
00164     proc = procIds.getRawPtr();
00165 
00166   offs[0] = 0;
00167   for (size_t i=0; i < numVtx; i++){
00168     offs[i+1] = 0;
00169     int vid = vtx ? vtx[i] : 0;
00170     for (lno_t j=allOffs[i]; j < allOffs[i+1]; j++){
00171       int owner = proc ? proc[j] : 0;
00172       bool keep = (!removeSelfEdges || vid != allIds[j]) &&
00173                (!removeOffProcessEdges || owner == myRank) &&
00174                (!removeOffGroupEdges || owner >= 0);
00175 
00176       if (keep)
00177         offs[i+1]++;
00178     }
00179   }
00180 
00181   // from counters to offsets
00182 
00183   for (size_t i=1; i < numVtx; i++)
00184     offs[i+1] += offs[i];
00185 
00186   numKeep = offs[numVtx];
00187 
00188   // do we need a new neighbor list?
00189 
00190   if (numNbors == numKeep){
00191     newGidNbors = Teuchos::arcpFromArrayView(gidNbors);
00192     newOffsets = Teuchos::arcpFromArrayView(offsets);
00193     return numNbors;
00194   }
00195   else if (numKeep == 0){
00196     newGidNbors = ArrayRCP<const gid_t>(Teuchos::null);
00197     newOffsets = offArray;
00198     return 0;
00199   }
00200 
00201   // Build the subset neighbor lists (id, weight, and offset).
00202 
00203   gid_t *newGids = new gid_t [numKeep];
00204   env->localMemoryAssertion(__FILE__, __LINE__, numKeep, newGids);
00205 
00206   newGidNbors = arcp(newGids, 0, numKeep, true);
00207   newOffsets = offArray;
00208 
00209   if (eDim > 0){
00210     newWeights = new scalar_t * [eDim];
00211     env->localMemoryAssertion(__FILE__, __LINE__, eDim, newWeights);
00212 
00213     if (numKeep) {
00214       for (int w=0; w < eDim; w++){
00215         newWeights[w] = new scalar_t [numKeep];
00216         env->localMemoryAssertion(__FILE__, __LINE__, numKeep, newWeights[w]);
00217       }
00218     }
00219     else {
00220       for (int w=0; w < eDim; w++)
00221         newWeights[w] = NULL;
00222     }
00223   }
00224 
00225   size_t next = 0;
00226   for (size_t i=0; i < numVtx && next < numKeep; i++){
00227     int vid = vtx ? vtx[i] : 0;
00228     for (lno_t j=allOffs[i]; j < allOffs[i+1]; j++){
00229       int owner = proc ? proc[j] : 0;
00230       bool keep = (!removeSelfEdges || vid != allIds[j]) &&
00231                (!removeOffProcessEdges || owner == myRank) &&
00232                (!removeOffGroupEdges || owner >= 0);
00233 
00234       if (keep){
00235         newGids[next] = allIds[j];
00236         for (int w=0; w < eDim; w++){
00237           newWeights[w][next] = edgeWeights[w][j];
00238         }
00239         next++;
00240         if (next == numKeep)
00241           break;
00242 
00243       }  // if (keep)
00244     }
00245   }
00246 
00247   return numKeep;
00248 }
00249 
00251 
00255 template <typename User>
00256 size_t computeLocalEdgeList(
00257   const RCP<const Environment> &env, int rank,
00258   size_t numLocalEdges,           // local edges
00259   size_t numLocalGraphEdges,      // edges in "local" graph
00260   RCP<const IdentifierMap<User> > &idMap,
00261   ArrayRCP<const typename InputTraits<User>::gid_t> &allEdgeIds, // in
00262   ArrayRCP<int> &allProcs,                                 // in
00263   ArrayRCP<const typename InputTraits<User>::lno_t> &allOffs,    // in
00264   ArrayRCP<StridedData<typename InputTraits<User>::lno_t,
00265     typename InputTraits<User>::scalar_t> > &allWeights,         // in
00266   ArrayRCP<const typename InputTraits<User>::lno_t> &edgeLocalIds, //
00267   ArrayRCP<const typename InputTraits<User>::lno_t> &offsets,      // out
00268   ArrayRCP<StridedData<typename InputTraits<User>::lno_t,
00269     typename InputTraits<User>::scalar_t> > &eWeights)             // out
00270 {
00271   typedef typename InputTraits<User>::gid_t gid_t;
00272   typedef typename InputTraits<User>::gno_t gno_t;
00273   typedef typename InputTraits<User>::scalar_t scalar_t;
00274   typedef typename InputTraits<User>::lno_t lno_t;
00275   typedef StridedData<lno_t, scalar_t> input_t;
00276 
00277   bool gnosAreGids = idMap->gnosAreGids();
00278 
00279   edgeLocalIds = ArrayRCP<const lno_t>(Teuchos::null);
00280   eWeights = ArrayRCP<input_t>(Teuchos::null);
00281   offsets = ArrayRCP<const lno_t>(Teuchos::null);
00282 
00283   if (numLocalGraphEdges == 0) {
00284     // Set the offsets array and return
00285     size_t allOffsSize = allOffs.size();
00286     lno_t *offs = new lno_t [allOffsSize];
00287     env->localMemoryAssertion(__FILE__, __LINE__, allOffsSize, offs);
00288     for (size_t i = 0; i < allOffsSize; i++) offs[i] = 0;
00289     offsets = arcp(offs, 0, allOffsSize, true);
00290     return 0;
00291   }
00292 
00293   if (numLocalGraphEdges == numLocalEdges){
00294 
00295     // Entire graph is local.
00296 
00297     lno_t *lnos = new lno_t [numLocalEdges];
00298     env->localMemoryAssertion(__FILE__, __LINE__,numLocalEdges, lnos);
00299     for (size_t i=0; i < numLocalEdges; i++)
00300       lnos[i] = allEdgeIds[i];
00301     edgeLocalIds = arcp(lnos, 0, numLocalEdges, true);
00302     offsets = allOffs;
00303     eWeights = allWeights;
00304   }
00305   else{
00306 
00307     // Create subset list of local graph edges, offsets and weights.
00308 
00309     int nWeightsPerEdge = allWeights.size();
00310 
00311     ArrayRCP<const gid_t> newEgids;
00312     scalar_t **newWeights = NULL;
00313 
00314     ArrayView<const gid_t> dummyVtx;
00315     ArrayView<const gid_t> nborView= allEdgeIds.view(0, numLocalEdges);
00316     ArrayView<const int> nborOwner = allProcs.view(0, numLocalEdges);
00317     ArrayView<input_t> eWgts = allWeights.view(0, nWeightsPerEdge);
00318     ArrayView<const lno_t> offView = allOffs.view(0, allOffs.size());
00319 
00320     try{
00321       numLocalEdges = removeUndesiredEdges<User>(env, rank, false, true, false,
00322                                                  dummyVtx, nborView, nborOwner,
00323                                                  eWgts, offView, newEgids,
00324                                                  newWeights, offsets);
00325     }
00326     Z2_FORWARD_EXCEPTIONS;
00327 
00328     env->localBugAssertion(__FILE__, __LINE__, "local graph miscalculation",
00329       numLocalEdges == numLocalGraphEdges, BASIC_ASSERTION);
00330 
00331     // offsets array was set by removeUndesiredEdges.  Create weight array.
00332 
00333     if (nWeightsPerEdge > 0){
00334       input_t *wgts = new input_t [nWeightsPerEdge];
00335       for (int w=0; w < nWeightsPerEdge; w++){
00336         ArrayRCP<const scalar_t> wgtArray(newWeights[w], 0, numLocalGraphEdges,true);
00337         wgts[w] = input_t(wgtArray, 1);
00338       }
00339       eWeights = arcp(wgts, 0, nWeightsPerEdge);
00340       delete [] newWeights;
00341     }
00342 
00343     // Create local ID array.  First translate gid to gno.
00344 
00345     ArrayRCP<gno_t> gnoArray;
00346 
00347     if (gnosAreGids){
00348       ArrayRCP<const gno_t> gnosConst =
00349         arcp_reinterpret_cast<const gno_t>(newEgids);
00350       gnoArray = arcp_const_cast<gno_t>(gnosConst);
00351     }
00352     else{
00353 
00354       ArrayRCP<gid_t> gidArray = arcp_const_cast<gid_t>(newEgids);
00355       gno_t *gnoList= new gno_t [numLocalGraphEdges];
00356       env->localMemoryAssertion(__FILE__, __LINE__, numLocalGraphEdges,
00357         gnoList);
00358       gnoArray = arcp(gnoList, 0, numLocalGraphEdges, true);
00359 
00360       try {
00361         idMap->gidTranslate(
00362           gidArray.view(0,numLocalGraphEdges),
00363           gnoArray.view(0,numLocalGraphEdges),
00364           TRANSLATE_APP_TO_LIB);
00365       }
00366       Z2_FORWARD_EXCEPTIONS;
00367     }
00368 
00369     // translate gno to lno
00370 
00371     lno_t *lnoList = new lno_t [numLocalGraphEdges];
00372     env->localMemoryAssertion(__FILE__, __LINE__, numLocalGraphEdges,
00373       lnoList);
00374     ArrayView<lno_t> lnoView(lnoList, numLocalGraphEdges);
00375 
00376     try {
00377       idMap->lnoTranslate(
00378         lnoView,
00379         gnoArray.view(0,numLocalGraphEdges),
00380         TRANSLATE_LIB_TO_APP);
00381     }
00382     Z2_FORWARD_EXCEPTIONS;
00383     edgeLocalIds = arcp<const lno_t>(lnoList, 0, numLocalGraphEdges, true);
00384   }
00385 
00386   return numLocalGraphEdges;
00387 }
00388 
00390 
00399 template <typename Adapter>
00400 class GraphModel : public Model<Adapter>
00401 {
00402 public:
00403 
00404 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00405   typedef typename Adapter::scalar_t    scalar_t;
00406   typedef typename Adapter::gno_t       gno_t;
00407   typedef typename Adapter::lno_t       lno_t;
00408   typedef typename Adapter::gid_t       gid_t;
00409   typedef typename Adapter::node_t      node_t;
00410   typedef typename Adapter::user_t      user_t;
00411   typedef typename Adapter::userCoord_t userCoord_t;
00412   typedef IdentifierMap<user_t>         idmap_t;
00413   typedef StridedData<lno_t, scalar_t>  input_t;
00414 #endif
00415 
00417   ~GraphModel() { }
00418 
00430   GraphModel(const MatrixAdapter<user_t,userCoord_t> *ia,
00431     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm,
00432     modelFlag_t &modelFlags);
00433 
00434   GraphModel(const GraphAdapter<user_t,userCoord_t> *ia,
00435     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm,
00436     modelFlag_t &modelFlags);
00437 
00438   GraphModel(const VectorAdapter<userCoord_t> *ia,
00439     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm,
00440     modelFlag_t &flags)
00441   {
00442     throw std::runtime_error("cannot build GraphModel from VectorAdapter");
00443   }
00444 
00445   GraphModel(const IdentifierAdapter<user_t> *ia,
00446     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm,
00447     modelFlag_t &flags)
00448   {
00449     throw std::runtime_error("cannot build GraphModel from IdentifierAdapter");
00450   }
00451 
00454   size_t getLocalNumVertices() const { return numLocalVertices_; }
00455 
00458   size_t getGlobalNumVertices() const { return numGlobalVertices_; }
00459 
00463   size_t getLocalNumGlobalEdges() const { return numLocalEdges_; }
00464 
00468   size_t getLocalNumLocalEdges() const { return numLocalGraphEdges_; }
00469 
00472   size_t getGlobalNumEdges() const { return numGlobalEdges_; }
00473 
00476   int getNumWeightsPerVertex() const { return numWeightsPerVertex_; }
00477 
00480   int getNumWeightsPerEdge() const { return nWeightsPerEdge_; }
00481 
00484   int getCoordinateDim() const { return vCoordDim_; }
00485 
00496   size_t getVertexList(
00497     ArrayView<const gno_t> &Ids,
00498     ArrayView<input_t> &xyz,
00499     ArrayView<input_t> &wgts) const 
00500   {
00501     size_t nv = gids_.size();
00502     if (gnosAreGids_) Ids = gids_.view(0, nv);
00503     else              Ids = gnosConst_.view(0, nv);
00504     xyz = vCoords_.view(0, vCoordDim_);
00505     wgts = vWeights_.view(0, numWeightsPerVertex_);
00506     return nv;
00507   }
00508 
00522   // Implied Vertex LNOs from getVertexList are used as indices to offsets
00523   // array.
00524   // Vertex GNOs are returned as neighbors in edgeIds.
00525 
00526   size_t getEdgeList( ArrayView<const gno_t> &edgeIds,
00527     ArrayView<const int> &procIds, ArrayView<const lno_t> &offsets,
00528     ArrayView<input_t> &wgts) const
00529   {
00530     if (gnosAreGids_) edgeIds = edgeGids_.view(0, numLocalEdges_);
00531     else              edgeIds = edgeGnosConst_.view(0, numLocalEdges_);
00532     procIds = procIdsConst_.view(0, numLocalEdges_);
00533     offsets = offsets_.view(0, numLocalVertices_+1);
00534     wgts = eWeights_.view(0, nWeightsPerEdge_);
00535     return numLocalEdges_;
00536   }
00537 
00561   size_t getLocalEdgeList(
00562     ArrayView<const lno_t> &edgeIds,
00563     ArrayView<const lno_t> &offsets,
00564     ArrayView<input_t> &wgts)
00565   {
00566     if (localGraphEdgeOffsets_.size() == 0) {
00567       // Local graph not created yet
00568       RCP<const IdentifierMap<user_t> > idmap = this->getIdentifierMap();
00569       computeLocalEdgeList(env_, comm_->getRank(),
00570         numLocalEdges_, numLocalGraphEdges_,
00571         idmap, edgeGids_, procIds_, offsets_, eWeights_,
00572         localGraphEdgeLnos_, localGraphEdgeOffsets_, localGraphEdgeWeights_);
00573     }
00574     edgeIds = localGraphEdgeLnos_();
00575     offsets = localGraphEdgeOffsets_();
00576     wgts = localGraphEdgeWeights_();
00577     return numLocalGraphEdges_;
00578   }
00579 
00581   // The Model interface.
00583 
00584   size_t getLocalNumObjects() const { return numLocalVertices_; }
00585 
00586   size_t getGlobalNumObjects() const { return numGlobalVertices_; }
00587 
00588 private:
00589   void shared_constructor(const Adapter *ia, modelFlag_t &modelFlags);
00590 
00591   const RCP<const Environment > env_;
00592   const RCP<const Comm<int> > comm_;
00593 
00594   ArrayRCP<const gid_t> gids_;        // vertices of input graph
00595   ArrayRCP<gno_t> gnos_;
00596 
00597   int numWeightsPerVertex_;
00598   ArrayRCP<input_t> vWeights_;
00599 
00600   int vCoordDim_;
00601   ArrayRCP<input_t> vCoords_;
00602 
00603   // Note: in case of graph subsetting, size of these arrays
00604   // may be larger than numLocalEdges_.  So do not use .size().
00605 
00606   ArrayRCP<const gid_t> edgeGids_;
00607   ArrayRCP<gno_t> edgeGnos_;
00608   ArrayRCP<int> procIds_;
00609   ArrayRCP<const lno_t> offsets_;
00610 
00611   int nWeightsPerEdge_;
00612   ArrayRCP<input_t> eWeights_;
00613 
00614   ArrayRCP<const gno_t> gnosConst_;
00615   ArrayRCP<const gno_t> edgeGnosConst_;
00616   ArrayRCP<const int> procIdsConst_;
00617 
00618   bool gnosAreGids_;
00619 
00620   // For local graphs (graph restricted to local process).  We
00621   // create these arrays only if required by the algorithm.
00622 
00623   ArrayRCP<const lno_t> localGraphEdgeLnos_;
00624   ArrayRCP<const lno_t> localGraphEdgeOffsets_;
00625   ArrayRCP<input_t> localGraphEdgeWeights_;
00626 
00627   // For convenience
00628 
00629   size_t numLocalVertices_;
00630   size_t numGlobalVertices_;
00631   size_t numLocalEdges_;
00632   size_t numGlobalEdges_;
00633   size_t numLocalGraphEdges_;
00634 };
00635 
00636 
00638 template <typename Adapter>
00639 GraphModel<Adapter>::GraphModel(
00640   const MatrixAdapter<user_t,userCoord_t> *ia,
00641   const RCP<const Environment> &env,
00642   const RCP<const Comm<int> > &comm,
00643   modelFlag_t &modelFlags):
00644        env_(env),
00645        comm_(comm),
00646        gids_(),
00647        gnos_(),
00648        numWeightsPerVertex_(0),
00649        vWeights_(),
00650        vCoordDim_(0),
00651        vCoords_(),
00652        edgeGids_(),
00653        edgeGnos_(),
00654        procIds_(),
00655        offsets_(),
00656        nWeightsPerEdge_(0),
00657        eWeights_(),
00658        gnosConst_(),
00659        edgeGnosConst_(),
00660        procIdsConst_(),
00661        gnosAreGids_(false),
00662        localGraphEdgeLnos_(),
00663        localGraphEdgeOffsets_(),
00664        localGraphEdgeWeights_(),
00665        numLocalVertices_(0),
00666        numGlobalVertices_(0),
00667        numLocalEdges_(0),
00668        numGlobalEdges_(0),
00669        numLocalGraphEdges_(0)
00670 {
00671   // Model creation flags
00672   bool symTranspose = modelFlags.test(SYMMETRIZE_INPUT_TRANSPOSE);
00673   bool symBipartite = modelFlags.test(SYMMETRIZE_INPUT_BIPARTITE);
00674   bool vertexCols = modelFlags.test(VERTICES_ARE_MATRIX_COLUMNS);
00675   bool vertexNz = modelFlags.test(VERTICES_ARE_MATRIX_NONZEROS);
00676 
00677   if (symTranspose || symBipartite || vertexCols || vertexNz){
00678     throw std::runtime_error("graph build option not yet implemented");
00679   }
00680 
00681   // Get the matrix from the input adapter
00682   gid_t const *vtxIds=NULL, *nborIds=NULL;
00683   lno_t const  *offsets=NULL;
00684   try{
00685     numLocalVertices_ = ia->getLocalNumIDs();
00686     ia->getIDsView(vtxIds);
00687   }
00688   Z2_FORWARD_EXCEPTIONS;
00689   try{
00690     if (ia->CRSViewAvailable()) {
00691       ia->getCRSView(offsets, nborIds);
00692     }
00693     else {
00694       // TODO:  Add support for CCS matrix layout
00695       throw std::runtime_error("Only MatrixAdapter::getCRSView is supported "
00696                                "in graph model");
00697     }
00698   }
00699   Z2_FORWARD_EXCEPTIONS;
00700 
00701   numLocalEdges_ = offsets[numLocalVertices_];
00702 
00703   gids_ = arcp<const gid_t>(vtxIds, 0, numLocalVertices_, false);
00704   edgeGids_ = arcp<const gid_t>(nborIds, 0, numLocalEdges_, false);
00705   offsets_ = arcp<const lno_t>(offsets, 0, numLocalVertices_ + 1, false);
00706 
00707   nWeightsPerEdge_ = 0;   // no edge weights from a matrix yet.
00708                      // TODO:  use matrix values as edge weights
00709 
00710   shared_constructor(ia, modelFlags);
00711 }
00712 
00713 
00715 template <typename Adapter>
00716 GraphModel<Adapter>::GraphModel(
00717   const GraphAdapter<user_t,userCoord_t> *ia,
00718   const RCP<const Environment> &env,
00719   const RCP<const Comm<int> > &comm,
00720   modelFlag_t &modelFlags):
00721        env_(env),
00722        comm_(comm),
00723        gids_(),
00724        gnos_(),
00725        numWeightsPerVertex_(0),
00726        vWeights_(),
00727        vCoordDim_(0),
00728        vCoords_(),
00729        edgeGids_(),
00730        edgeGnos_(),
00731        procIds_(),
00732        offsets_(),
00733        nWeightsPerEdge_(0),
00734        eWeights_(),
00735        gnosConst_(),
00736        edgeGnosConst_(),
00737        procIdsConst_(),
00738        gnosAreGids_(false),
00739        localGraphEdgeLnos_(),
00740        localGraphEdgeOffsets_(),
00741        localGraphEdgeWeights_(),
00742        numLocalVertices_(0),
00743        numGlobalVertices_(0),
00744        numLocalEdges_(0),
00745        numGlobalEdges_(0),
00746        numLocalGraphEdges_(0)
00747 {
00748 
00749   // This GraphModel is built with vertices == GRAPH_VERTEX from GraphAdapter.
00750   // It is not ready to use vertices == GRAPH_EDGE from GraphAdapter.
00751   env_->localInputAssertion(__FILE__, __LINE__,
00752     "GraphModel from GraphAdapter is implemented only for "
00753     "Graph Vertices as primary object, not for Graph Edges",
00754     ia->getPrimaryEntityType() == Zoltan2::GRAPH_VERTEX, BASIC_ASSERTION);
00755 
00756   // Get the graph from the input adapter
00757 
00758   gid_t const *vtxIds=NULL, *nborIds=NULL;
00759   lno_t const  *offsets=NULL;
00760   try{
00761     numLocalVertices_ = ia->getLocalNumVertices();
00762     ia->getVertexIDsView(vtxIds);
00763     ia->getEdgesView(offsets, nborIds);
00764   }
00765   Z2_FORWARD_EXCEPTIONS;
00766 
00767   numLocalEdges_ = offsets[numLocalVertices_];
00768 
00769   gids_ = arcp<const gid_t>(vtxIds, 0, numLocalVertices_, false);
00770   edgeGids_ = arcp<const gid_t>(nborIds, 0, numLocalEdges_, false);
00771   offsets_ = arcp<const lno_t>(offsets, 0, numLocalVertices_ + 1, false);
00772 
00773   nWeightsPerEdge_ = ia->getNumWeightsPerEdge();
00774 
00775   if (nWeightsPerEdge_ > 0){
00776     input_t *wgts = new input_t [nWeightsPerEdge_];
00777     eWeights_ = arcp(wgts, 0, nWeightsPerEdge_, true);
00778   }
00779 
00780   for (int w=0; w < nWeightsPerEdge_; w++){
00781     const scalar_t *ewgts=NULL;
00782     int stride=0;
00783 
00784     ia->getEdgeWeightsView(ewgts, stride, w);
00785 
00786     ArrayRCP<const scalar_t> wgtArray(ewgts, 0, numLocalEdges_, false);
00787     eWeights_[w] = input_t(wgtArray, stride);
00788   }
00789 
00790   shared_constructor(ia, modelFlags);
00791 }
00792 
00794 template <typename Adapter>
00795 void GraphModel<Adapter>::shared_constructor(
00796   const Adapter *ia,
00797   modelFlag_t &modelFlags)
00798 {
00799   // Model creation flags
00800   bool consecutiveIdsRequired =
00801     modelFlags.test(IDS_MUST_BE_GLOBALLY_CONSECUTIVE);
00802   bool removeSelfEdges = modelFlags.test(SELF_EDGES_MUST_BE_REMOVED);
00803   bool subsetGraph = modelFlags.test(GRAPH_IS_A_SUBSET_GRAPH);
00804 
00805   // A subset graph is special in that it may contain neighbor
00806   // vertices that are not owned by processes in this communicator.
00807   // We remove these.
00808 
00809   ArrayRCP<const int> nborProcs;
00810 
00811   if (subsetGraph){
00812     RCP<const idmap_t> idMap;
00813     try{
00814       idMap = rcp(new idmap_t(env_, comm_, gids_, false));
00815     }
00816     Z2_FORWARD_EXCEPTIONS;
00817 
00818     ArrayRCP<int> procArray;
00819 
00820     if (numLocalEdges_ > 0){
00821       int *pids = new int [numLocalEdges_];
00822       env_->localMemoryAssertion(__FILE__, __LINE__, numLocalEdges_, pids);
00823       procArray = arcp(pids, 0, numLocalEdges_, true);
00824     }
00825 
00826     ArrayView<gno_t> dummyGno;
00827 
00828     try{
00829       // All processes must make this call.
00830       // procOwner will be -1 if edge Id is not in our communicator.
00831 
00832       idMap->gidGlobalTranslate( edgeGids_.view(0, numLocalEdges_),
00833         dummyGno, procArray.view(0, numLocalEdges_));
00834     }
00835     Z2_FORWARD_EXCEPTIONS;
00836 
00837     int outOfSubset = 0;
00838     for (size_t i=0; i < numLocalEdges_; i++){
00839       if (procArray[i] < 0){
00840         outOfSubset++;
00841         break;
00842       }
00843     }
00844 
00845     if (outOfSubset == 0){
00846       procArray.clear();
00847       subsetGraph = false;
00848     }
00849     else{
00850       nborProcs = arcp_const_cast<const int>(procArray);
00851     }
00852   }
00853 
00854   // Now remove undesired edges.
00855 
00856   if (subsetGraph || removeSelfEdges){
00857 
00858     ArrayRCP<const gid_t> newEdges;
00859     ArrayRCP<const lno_t> newOffs;
00860     scalar_t **newWeights = NULL;
00861     size_t numNewEdges = 0;
00862 
00863     ArrayView<const gid_t> vtxView = gids_.view(0, numLocalVertices_);
00864     ArrayView<const gid_t> nborView= edgeGids_.view(0, numLocalEdges_);
00865     ArrayView<const int> nborOwner = nborProcs.view(0, nborProcs.size());
00866     ArrayView<input_t> eWgts = eWeights_.view(0, nWeightsPerEdge_);
00867     ArrayView<const lno_t> offView = offsets_.view(0, numLocalVertices_ + 1);
00868 
00869     try{
00870       numNewEdges = removeUndesiredEdges<user_t>(env_, comm_->getRank(),
00871         removeSelfEdges,
00872         false,
00873         subsetGraph,
00874         vtxView,
00875         nborView,
00876         nborOwner,
00877         eWgts,
00878         offView,
00879         newEdges,
00880         newWeights,
00881         newOffs);
00882     }
00883     Z2_FORWARD_EXCEPTIONS;
00884 
00885     nborProcs.clear();
00886 
00887     if (numNewEdges < numLocalEdges_){
00888       edgeGids_ = newEdges;
00889       offsets_ = newOffs;
00890       numLocalEdges_ = numNewEdges;
00891 
00892       for (int w=0; w < nWeightsPerEdge_; w++){
00893         ArrayRCP<const scalar_t> wgtArray(newWeights[w], 0, numNewEdges, true);
00894         eWeights_[w] = input_t(wgtArray, 1);
00895       }
00896     }
00897     delete [] newWeights;
00898   }
00899 
00900   // Create an IdentifierMap, which maps the user's global IDs to
00901   //   Zoltan2 internal global numbers if necessary.
00902   //   The map can also give us owners of our vertex neighbors.
00903 
00904   RCP<const idmap_t> idMap;
00905 
00906   try{
00907     idMap = rcp(new idmap_t(env_, comm_, gids_, consecutiveIdsRequired));
00908   }
00909   Z2_FORWARD_EXCEPTIONS;
00910 
00911   // Model base class needs to have IdentifierMap.
00912 
00913   this->setIdentifierMap(idMap);
00914 
00915   numGlobalVertices_ = idMap->getGlobalNumberOfIds();
00916   gnosAreGids_ = idMap->gnosAreGids();
00917 
00918   // Compute internal global numbers if we can not use the
00919   // user's global Ids.  Also find the process owning each
00920   // neighboring vertex.
00921 
00922   ArrayView<const gid_t> gidArray(Teuchos::null);  // edge gid
00923   ArrayView<gno_t> gnoArray(Teuchos::null);        // edge gno
00924   ArrayView<int> procArray(Teuchos::null);         // edge owner
00925 
00926   if (numLocalVertices_){
00927 
00928     if (!gnosAreGids_){   // need vertex global numbers, edge global numbers
00929       gno_t *tmp = new gno_t [numLocalVertices_];
00930       env_->localMemoryAssertion(__FILE__, __LINE__, numLocalVertices_, tmp);
00931       gnos_ = arcp(tmp, 0, numLocalVertices_);
00932 
00933       try{
00934         ArrayRCP<gid_t> tmpGids = arcp_const_cast<gid_t>(gids_);
00935 
00936         idMap->gidTranslate(tmpGids(0,numLocalVertices_),
00937           gnos_(0,numLocalVertices_), TRANSLATE_APP_TO_LIB);
00938       }
00939       Z2_FORWARD_EXCEPTIONS;
00940 
00941       if (numLocalEdges_){
00942         tmp = new gno_t [numLocalEdges_];
00943         env_->localMemoryAssertion(__FILE__, __LINE__, numLocalEdges_, tmp);
00944         edgeGnos_ = arcp(tmp, 0, numLocalEdges_);
00945         gnoArray = edgeGnos_.view(0, numLocalEdges_);
00946       }
00947     }
00948 
00949     if (numLocalEdges_){
00950       gidArray = edgeGids_.view(0, numLocalEdges_);
00951 
00952       int *p = new int [numLocalEdges_];
00953       env_->localMemoryAssertion(__FILE__, __LINE__, numLocalEdges_, p);
00954       procIds_ = arcp(p, 0, numLocalEdges_);
00955       procArray = procIds_.view(0, numLocalEdges_);
00956     }
00957   }
00958 
00959   try{
00960     // All processes must make this call.
00961     idMap->gidGlobalTranslate(gidArray, gnoArray, procArray);
00962   }
00963   Z2_FORWARD_EXCEPTIONS;
00964 
00965   gnosConst_ = arcp_const_cast<const gno_t>(gnos_);
00966   edgeGnosConst_ = arcp_const_cast<const gno_t>(edgeGnos_);
00967   procIdsConst_ = arcp_const_cast<const int>(procIds_);
00968 
00969   // Number of edges such that neighbor is on the local process.
00970   // We only create the list of local graph edges if the user
00971   // calls getLocalEdgeList().
00972 
00973   numLocalGraphEdges_ = 0;
00974   int *pids = procArray.getRawPtr();
00975   int me = comm_->getRank();
00976   for (size_t i=0; i < numLocalEdges_; i++)
00977     if (pids[i] == me) numLocalGraphEdges_++;
00978 
00979   // Vertex weights
00980 
00981   numWeightsPerVertex_ = ia->getNumWeightsPerID();
00982 
00983   if (numWeightsPerVertex_ > 0){
00984     input_t *weightInfo = new input_t [numWeightsPerVertex_];
00985     env_->localMemoryAssertion(__FILE__, __LINE__, numWeightsPerVertex_, weightInfo);
00986 
00987     for (int idx=0; idx < numWeightsPerVertex_; idx++){
00988       bool useNumNZ = ia->useDegreeAsWeight(idx); 
00989       if (useNumNZ){
00990         scalar_t *wgts = new scalar_t [numLocalVertices_];
00991         env_->localMemoryAssertion(__FILE__, __LINE__, numLocalVertices_, wgts);
00992         ArrayRCP<const scalar_t> wgtArray =
00993           arcp(wgts, 0, numLocalVertices_, true);
00994         for (size_t i=0; i < numLocalVertices_; i++){
00995           wgts[i] = offsets_[i+1] - offsets_[i];
00996         }
00997         weightInfo[idx] = input_t(wgtArray, 1);
00998       }
00999       else{
01000         const scalar_t *weights=NULL;
01001         int stride=0;
01002         ia->getWeightsView(weights, stride, idx);
01003         ArrayRCP<const scalar_t> wgtArray = arcp(weights, 0,
01004                                                  stride*numLocalVertices_,
01005                                                  false);
01006         weightInfo[idx] = input_t(wgtArray, stride);
01007       }
01008     }
01009 
01010     vWeights_ = arcp<input_t>(weightInfo, 0, numWeightsPerVertex_, true);
01011   }
01012 
01013 
01014   // Vertex coordinates
01015 
01016   if (ia->coordinatesAvailable()) {
01017     vCoordDim_ = ia->getCoordinateInput()->getNumEntriesPerID();
01018 
01019     if (vCoordDim_ > 0){
01020       input_t *coordInfo = new input_t [vCoordDim_];
01021       env_->localMemoryAssertion(__FILE__, __LINE__, vCoordDim_, coordInfo);
01022 
01023       for (int dim=0; dim < vCoordDim_; dim++){
01024         const scalar_t *coords=NULL;
01025         int stride=0;
01026         ia->getCoordinateInput()->getEntriesView(coords, stride, dim);
01027         ArrayRCP<const scalar_t> coordArray = arcp(coords, 0,
01028                                                    stride*numLocalVertices_,
01029                                                    false);
01030         coordInfo[dim] = input_t(coordArray, stride);
01031       }
01032 
01033       vCoords_ = arcp<input_t>(coordInfo, 0, vCoordDim_, true);
01034     }
01035   }
01036 
01037   reduceAll<int, size_t>(*comm_, Teuchos::REDUCE_SUM, 1,
01038     &numLocalEdges_, &numGlobalEdges_);
01039 
01040   env_->memory("After construction of graph model");
01041 }
01042 
01043 
01044 }   // namespace Zoltan2
01045 
01046 #endif
01047