C++ Epetra_BlockMap::GID方法代码示例

int
LOCA::Epetra::AugmentedOp::blockMap2PointMap(const Epetra_BlockMap& BlockMap,
				    Epetra_Map*& PointMap) const
{
  // Generate an Epetra_Map that has the same number and distribution of points
  // as the input Epetra_BlockMap object.  The global IDs for the output PointMap
  // are computed by using the MaxElementSize of the BlockMap.  For variable block
  // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.

  int MaxElementSize = BlockMap.MaxElementSize();
  int PtNumMyElements = BlockMap.NumMyPoints();
  int * PtMyGlobalElements = 0;
  if (PtNumMyElements>0) PtMyGlobalElements = new int[PtNumMyElements];

  int NumMyElements = BlockMap.NumMyElements();

  int curID = 0;
  for (int i=0; i<NumMyElements; i++) {
    int StartID = BlockMap.GID(i)*MaxElementSize;
    int ElementSize = BlockMap.ElementSize(i);
    for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
  }
  assert(curID==PtNumMyElements); // Sanity test

  PointMap = new Epetra_Map(-1, PtNumMyElements, PtMyGlobalElements, BlockMap.IndexBase(), BlockMap.Comm());

  if (PtNumMyElements>0) delete [] PtMyGlobalElements;

  if (!BlockMap.PointSameAs(*PointMap)) {EPETRA_CHK_ERR(-1);} // Maps not compatible
  return(0);
}

//=========================================================================
int Ifpack_CrsRiluk::BlockMap2PointMap(const Epetra_BlockMap & BlockMap, Teuchos::RefCountPtr<Epetra_Map>* PointMap) {
	// Generate an Epetra_Map that has the same number and distribution of points
	// as the input Epetra_BlockMap object.  The global IDs for the output PointMap
	// are computed by using the MaxElementSize of the BlockMap.  For variable block
	// sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.

	int MaxElementSize = BlockMap.MaxElementSize();
	int PtNumMyElements = BlockMap.NumMyPoints();
	vector<int> PtMyGlobalElements;
	if (PtNumMyElements>0) PtMyGlobalElements.resize(PtNumMyElements);

	int NumMyElements = BlockMap.NumMyElements();

	int curID = 0;
	for (int i=0; i<NumMyElements; i++) {
		int StartID = BlockMap.GID(i)*MaxElementSize;
		int ElementSize = BlockMap.ElementSize(i);
		for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
	}
	assert(curID==PtNumMyElements); // Sanity test

	(*PointMap) = Teuchos::rcp( new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements[0], BlockMap.IndexBase(), BlockMap.Comm()) );

	if (!BlockMap.PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);} // Maps not compatible
  return(0);
}

Teuchos::Array<int>
Albany::NodeSetSolutionCullingStrategy::
selectedGIDs(const Epetra_BlockMap &sourceMap) const
{
  Teuchos::Array<int> result;
  {
    Teuchos::Array<int> mySelectedGIDs;
    {
      const NodeSetList &nodeSets = disc_->getNodeSets();

      const NodeSetList::const_iterator it = nodeSets.find(nodeSetLabel_);
      if (it != nodeSets.end()) {
        typedef NodeSetList::mapped_type NodeSetEntryList;
        const NodeSetEntryList &sampleNodeEntries = it->second;

        for (NodeSetEntryList::const_iterator jt = sampleNodeEntries.begin(); jt != sampleNodeEntries.end(); ++jt) {
          typedef NodeSetEntryList::value_type NodeEntryList;
          const NodeEntryList &sampleEntries = *jt;
          for (NodeEntryList::const_iterator kt = sampleEntries.begin(); kt != sampleEntries.end(); ++kt) {
            mySelectedGIDs.push_back(sourceMap.GID(*kt));
          }
        }
      }
    }

    const Epetra_Comm &comm = sourceMap.Comm();

    {
      int selectedGIDCount;
      {
        int mySelectedGIDCount = mySelectedGIDs.size();
        comm.SumAll(&mySelectedGIDCount, &selectedGIDCount, 1);
      }
      result.resize(selectedGIDCount);
    }

    const int ierr = Epetra::GatherAllV(
        comm,
        mySelectedGIDs.getRawPtr(), mySelectedGIDs.size(),
        result.getRawPtr(), result.size());
    TEUCHOS_ASSERT(ierr == 0);
  }

  std::sort(result.begin(), result.end());

  return result;
}

int BlockMapToHandle(FILE * handle, const Epetra_BlockMap & map) {

  const Epetra_Comm & comm = map.Comm();
  int numProc = comm.NumProc();
  bool doSizes = !map.ConstantElementSize();

  if (numProc==1) {
    int * myElements = map.MyGlobalElements();
    int * elementSizeList = 0;
    if (doSizes) elementSizeList = map.ElementSizeList();
    return(writeBlockMap(handle, map.NumGlobalElements(), myElements, elementSizeList, doSizes));
  }

  int numRows = map.NumMyElements();
  
  Epetra_Map allGidsMap(-1, numRows, 0,comm);
  
  Epetra_IntVector allGids(allGidsMap);
  for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
  
  Epetra_IntVector allSizes(allGidsMap);
  for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
  
  // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
  int numChunks = numProc;
  int stripSize = allGids.GlobalLength()/numChunks;
  int remainder = allGids.GlobalLength()%numChunks;
  int curStart = 0;
  int curStripSize = 0;
  Epetra_IntSerialDenseVector importGidList;
  Epetra_IntSerialDenseVector importSizeList;
  if (comm.MyPID()==0) {
    importGidList.Size(stripSize+1); // Set size of vector to max needed
    if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
  }
  for (int i=0; i<numChunks; i++) {
    if (comm.MyPID()==0) { // Only PE 0 does this part
      curStripSize = stripSize;
      if (i<remainder) curStripSize++; // handle leftovers
      for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
      curStart += curStripSize;
    }
    // The following import map will be non-trivial only on PE 0.
    Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
    Epetra_Import gidImporter(importGidMap, allGidsMap);
    
    Epetra_IntVector importGids(importGidMap);
    if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
    Epetra_IntVector importSizes(importGidMap);
    if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1); 
    
    // importGids (and importSizes, if non-trivial block map)
    // now have a list of GIDs (and sizes, respectively) for the current strip of map.
    
    int * myElements = importGids.Values();
    int * elementSizeList = 0;
    if (doSizes) elementSizeList = importSizes.Values();
    // Finally we are ready to write this strip of the map to file
    writeBlockMap(handle, importGids.MyLength(), myElements, elementSizeList, doSizes);
  }
  return(0);
}

static int compute_hypergraph_metrics(const Epetra_BlockMap &rowmap, 
            const Epetra_BlockMap &colmap,
            int numGlobalColumns,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, double &cutn, double &cutl)  // output
{
  const Epetra_Comm &comm  = rowmap.Comm();
#ifdef HAVE_MPI
  const Epetra_MpiComm* mpiComm =
    dynamic_cast<const Epetra_MpiComm*>(&comm);

  MPI_Comm mcomm = mpiComm->Comm();
#endif
  int nProcs = comm.NumProc();
  int myProc = comm.MyPID();
  double min, avg;
  std::map<int, float> vertexWeights;
  std::map<int, std::map<int, float > > graphEdgeWeights;
  std::map<int, float> hyperEdgeWeights;

  costs.getCosts(vertexWeights,  // vertex global ID -> weight
           graphEdgeWeights,    //  vertex global ID -> map from neighbor global ID to edge weight
           hyperEdgeWeights);   // hyperedge global ID -> weight

  Epetra_Vector vwgt(rowmap);
  int numVWgts = vertexWeights.size();

  if (numVWgts > 0){
    double *wvals = new double [numVWgts];
    int *gids = new int [numVWgts];

    std::map<int, float>::iterator vnext = vertexWeights.begin();
    int i=0;
    while (vnext != vertexWeights.end()){
      wvals[i] =  vnext->second;
      gids[i] = vnext->first;
      vnext++;
      i++;
    }

    vwgt.ReplaceGlobalValues(i, wvals, gids);

    delete [] wvals;
    delete [] gids;
  }
  else{
    vwgt.PutScalar(1.0);   // default to unit weights
  }

  compute_balance(vwgt, myGoalWeight, min, balance, avg);

  if (balance < 0){
    return 1;
  }

  /* Compute cutl and cutn. 
   */

  int totalHEWeights = 0; 

  int numHEWeights = hyperEdgeWeights.size();

  comm.SumAll(&numHEWeights, &totalHEWeights, 1);
 
  if ((totalHEWeights > 0) && (totalHEWeights <  numGlobalColumns)){
    if (myProc == 0)
      std::cerr << "Must supply either no h.e. weights or else supply at least one for each column" << std::endl;
      return -1;
  }

  std::map<int, float>::iterator heWgtIter;

  // Create a set containing all the columns in my rows.  We assume all
  // the rows are in the same partition.

  int numMyCols = colmap.NumMyElements();

  std::set<int> colGIDS;
  std::set<int>::iterator gidIter;

  for (int j=0; j<numMyCols; j++){
    colGIDS.insert(colmap.GID(j));
  }
  
  /* Divide columns among processes, then each process computes its
   * assigned columns' cutl and cutn.
   *  TODO - numGlobalColumns can be less than nprocs

   * Fix this when a process is assigned no columns. TODO
   */
  int ncols = numGlobalColumns / nProcs;
  int leftover = numGlobalColumns - (nProcs * ncols);
  std::vector<int> colCount(nProcs, 0);
  for (int i=0; i<nProcs; i++){
    colCount[i] = ncols;
    if (i < leftover) colCount[i]++;
  }
  int *colTotals = NULL;
  double *colWeights = NULL;
//.........这里部分代码省略.........

static int compute_graph_metrics(const Epetra_BlockMap &rowmap, const Epetra_BlockMap &colmap,
            std::vector<std::vector<int> > &rows, 
            Isorropia::Epetra::CostDescriber &costs, double &myGoalWeight,
            double &balance, int &numCuts, double &cutWgt, double &cutn, double &cutl)
{
  const Epetra_Comm &comm  = rowmap.Comm();
  int myProc = comm.MyPID();
  int myCols = colmap.NumMyElements();
  double min, avg;
  std::map<int, float> vertexWeights;
  std::map<int, std::map<int, float > > graphEdgeWeights;
  std::map<int, float> hyperEdgeWeights;

  costs.getCosts(vertexWeights,  // vertex global ID -> weight
           graphEdgeWeights,    //  vertex global ID -> map from neighbor global ID to edge weight
           hyperEdgeWeights);   // hyperedge global ID -> weight

  // Compute the balance

  Epetra_Vector vwgt(rowmap);
  int numVWgts = vertexWeights.size();

  if (numVWgts > 0){
    double *wvals = new double [numVWgts];
    int *gids = new int [numVWgts];

    std::map<int, float>::iterator vnext = vertexWeights.begin();
    int i=0;
    while (vnext != vertexWeights.end()){
      wvals[i] =  vnext->second;
      gids[i] = vnext->first;
      vnext++;
      i++;
    }

    vwgt.ReplaceGlobalValues(i, wvals, gids);

    delete [] wvals;
    delete [] gids;
  }
  else{
    vwgt.PutScalar(1.0);   // default to unit weights
  }

  compute_balance(vwgt, myGoalWeight, min, balance, avg);

  if (balance < 0){
    return 1;
  }

  // Compute the measures based on cut edges

  int *procID = new int [myCols];
  int *GID = new int [myCols];
  int *tmp = new int [myCols];

  for (int i=0; i < myCols; i++){
    GID[i] = colmap.GID(i);
  }

  rowmap.RemoteIDList(myCols, GID, procID, tmp);  // matrix is square

  delete [] tmp;

  int haveEdgeWeights = graphEdgeWeights.size();

  int localNumCuts = 0;
  double localCutWgt = 0.0;
  double localCutn = 0.0;
  double localCutl = 0.0;

  for (int i=0; i < rowmap.NumMyElements(); i++){
    int vtxGID = rowmap.GID(i);

    int numEdges = rows[i].size();

    if (numEdges > 0){

      std::map<int, std::map<int, float> >::iterator wnext;
      if (haveEdgeWeights){
        wnext = graphEdgeWeights.find(vtxGID);
        if (wnext == graphEdgeWeights.end()){
          std::cerr << "Graph edge weights are missing for vertex " << vtxGID;
          std::cerr << std::endl;
          return -1;
        }
      }

      double heWeight = 0.0;
      std::set<int> nbors;

      for (int j=0; j < numEdges; j++){

        int colGID = GID[rows[i][j]];

        int nborProc = procID[rows[i][j]];

        if (colGID == vtxGID) continue;  // skip self edges

        float wgt = 1.0;
//.........这里部分代码省略.........

int checkmap(Epetra_BlockMap & Map, int NumGlobalElements, int NumMyElements, 
	     int *MyGlobalElements, int ElementSize, int * ElementSizeList,
	     int NumGlobalPoints, int NumMyPoints,
	     int IndexBase, Epetra_Comm& Comm,
	     bool DistributedGlobal,
	     bool IsOneToOne)
{

  int i, ierr=0, forierr=0;// forierr is used in for loops, then is tested
  // after for loop completes to see if it is non zero - potentially prevents
  // thousands of error messages

  if (ElementSizeList==0)
    {
      EPETRA_TEST_ERR(!Map.ConstantElementSize(),ierr);
    }
  else
    EPETRA_TEST_ERR(Map.ConstantElementSize(),ierr);
  
  EPETRA_TEST_ERR(DistributedGlobal!=Map.DistributedGlobal(),ierr);

  EPETRA_TEST_ERR(IsOneToOne!=Map.IsOneToOne(),ierr);

  int *MyElementSizeList;

  if (ElementSizeList==0)
    {
      EPETRA_TEST_ERR(Map.ElementSize()!=ElementSize,ierr);
      
      MyElementSizeList = new int[NumMyElements];
      
      EPETRA_TEST_ERR(Map.ElementSizeList(MyElementSizeList)!=0,ierr);
      forierr = 0;
      for (i=0; i<NumMyElements; i++) 
        forierr += MyElementSizeList[i]!=ElementSize;
      EPETRA_TEST_ERR(forierr,ierr);

      EPETRA_TEST_ERR(Map.MaxMyElementSize() != ElementSize,ierr);
      EPETRA_TEST_ERR(Map.MinMyElementSize() != ElementSize,ierr);
    }
  else
    {
      MyElementSizeList = new int[NumMyElements];
      EPETRA_TEST_ERR(Map.ElementSizeList(MyElementSizeList)!=0,ierr);
      int MaxSize = MyElementSizeList[0];
      int MinSize = MyElementSizeList[0];
      forierr=0;
      for (i=0; i<NumMyElements; i++) {
        forierr += MyElementSizeList[i]!=ElementSizeList[i];
	if (MyElementSizeList[i] > MaxSize)
	  MaxSize = MyElementSizeList[i];
	if (MyElementSizeList[i] < MinSize)
	  MinSize = MyElementSizeList[i];

	// Test ElementSize(int LID) method	

	forierr += Map.ElementSize(Map.LID(MyGlobalElements[i])) != ElementSizeList[i];
      }
      EPETRA_TEST_ERR(forierr,ierr);
   
      EPETRA_TEST_ERR(MaxSize !=Map.MaxMyElementSize(),ierr);
      EPETRA_TEST_ERR(MinSize !=Map.MinMyElementSize(),ierr);
    }

  const Epetra_Comm & Comm1 = Map.Comm();

  EPETRA_TEST_ERR(Comm1.NumProc()!=Comm.NumProc(),ierr);

  EPETRA_TEST_ERR(Comm1.MyPID()!=Comm.MyPID(),ierr);

  EPETRA_TEST_ERR(Map.IndexBase()!=IndexBase,ierr);

  EPETRA_TEST_ERR(!Map.LinearMap() && MyGlobalElements==0,ierr);

  EPETRA_TEST_ERR(Map.LinearMap() && MyGlobalElements!=0,ierr);

  EPETRA_TEST_ERR(Map.MaxAllGID()!=NumGlobalElements-1+IndexBase,ierr);

  EPETRA_TEST_ERR(Map.MaxElementSize()!=ElementSize,ierr);

  int MaxLID = Map.MaxLID();
  EPETRA_TEST_ERR(MaxLID!=NumMyElements-1,ierr);

  int MaxMyGID = (Comm.MyPID()+1)*NumMyElements-1+IndexBase;
  if (Comm.MyPID()>2) MaxMyGID+=3;
  if (!DistributedGlobal) MaxMyGID = NumMyElements-1+IndexBase;
  EPETRA_TEST_ERR(Map.MaxMyGID()!=MaxMyGID,ierr);

  EPETRA_TEST_ERR(Map.MinAllGID()!=IndexBase,ierr);

  if (ElementSizeList==0)
    {
      EPETRA_TEST_ERR(Map.MinElementSize()!=ElementSize,ierr);
    }
  else EPETRA_TEST_ERR(Map.MinElementSize()!=2,ierr);

  int MinLID = Map.MinLID();
  EPETRA_TEST_ERR(MinLID!=0,ierr);

  int MinMyGID = Comm.MyPID()*NumMyElements+IndexBase;
//.........这里部分代码省略.........