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

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);
}

// ============================================================================
std::shared_ptr<const Tpetra::Map<int,int>>
BorderingHelpers::
extendMapBy1(const Epetra_BlockMap & map)
{
  const Teuchos::Comm<int> & comm = map.Comm();
  // Create a new map that hosts one more entry.
  const int numGlobalElements = map.NumGlobalElements() + 1;
  const int numMyElements = map.NumMyElements();
  int * myGlobalElements = map.MyGlobalElements();
  // The following if-else construction just makes sure that
  // the Tpetra::Map<int,int> constructor is called with an extended
  // map on proc 0, and with the regular old stuff on all
  // other procs.
  std::shared_ptr<Tpetra::Map<int,int>> extendedMap;
  if (comm.MyPID() == 0) {
    // Copy over the global indices.
    std::vector<int> a(numMyElements+1);
    for (int k = 0; k < numMyElements; k++)
      a[k] = myGlobalElements[k];
    // Append one more.
    a[numMyElements] = map.NumGlobalElements();

    extendedMap = std::make_shared<Tpetra::Map<int,int>>(
        numGlobalElements,
        numMyElements+1,
        &a[0],
        map.IndexBase(),
        comm
        );
  } else {
    extendedMap = std::make_shared<Tpetra::Map<int,int>>(
        numGlobalElements,
        numMyElements,
        myGlobalElements,
        map.IndexBase(),
        comm
        );
  }

  return extendedMap;
}

// FIXME long long
Epetra_BlockMap
Epetra_Util::Create_OneToOne_BlockMap(const Epetra_BlockMap& usermap,
              bool high_rank_proc_owns_shared)
{
// FIXME long long

  //if usermap is already 1-to-1 then we'll just return a copy of it.
  if (usermap.IsOneToOne()) {
    Epetra_BlockMap newmap(usermap);
    return(newmap);
  }

  int myPID = usermap.Comm().MyPID();
  Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);

  int numMyElems = usermap.NumMyElements();
  const int* myElems = usermap.MyGlobalElements();

  int* owner_procs = new int[numMyElems*2];
  int* sizes = owner_procs+numMyElems;

  directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
         0, sizes, high_rank_proc_owns_shared);

  //we'll fill a list of map-elements which belong on this processor

  int* myOwnedElems = new int[numMyElems*2];
  int* ownedSizes = myOwnedElems+numMyElems;
  int numMyOwnedElems = 0;

  for(int i=0; i<numMyElems; ++i) {
    int GID = myElems[i];
    int owner = owner_procs[i];

    if (myPID == owner) {
      ownedSizes[numMyOwnedElems] = sizes[i];
      myOwnedElems[numMyOwnedElems++] = GID;
    }
  }

  Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
         sizes, usermap.IndexBase(), usermap.Comm());

  delete [] myOwnedElems;
  delete [] owner_procs;
  delete directory;

  return(one_to_one_map);
}

//EpetraMap_To_TpetraMap: takes in Epetra_Map object, converts it to its equivalent Tpetra::Map object,
//and returns an RCP pointer to this Tpetra::Map
Teuchos::RCP<const Tpetra_Map> Petra::EpetraMap_To_TpetraMap(const Epetra_BlockMap& epetraMap_,
                                                      const Teuchos::RCP<const Teuchos::Comm<int> >& commT_)
{
  const std::size_t numElements = Teuchos::as<std::size_t>(epetraMap_.NumMyElements());
  const auto indexBase = Teuchos::as<GO>(epetraMap_.IndexBase());
  if (epetraMap_.DistributedGlobal() || epetraMap_.Comm().NumProc() == Teuchos::OrdinalTraits<int>::one()) {
    Teuchos::Array<Tpetra_GO> indices(numElements);
    int *epetra_indices = epetraMap_.MyGlobalElements();
    for(LO i=0; i < numElements; i++)
       indices[i] = epetra_indices[i];
    const Tpetra::global_size_t computeGlobalElements = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
    return Teuchos::rcp(new Tpetra_Map(computeGlobalElements, indices, indexBase, commT_));
  } else {
    return Teuchos::rcp(new Tpetra_Map(numElements, indexBase, commT_, Tpetra::LocallyReplicated));
  }
}

//=========================================================================
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);
}

int BlockMapToMatrixMarketFile( const char *filename, const Epetra_BlockMap & map, 
				 const char * mapName,
				 const char *mapDescription, 
				 bool writeHeader) {
  int M = map.NumGlobalElements();
  int N = 1;
  if (map.MaxElementSize()>1) N = 2; // Non-trivial block map, store element sizes in second column

  FILE * handle = 0;

  if (map.Comm().MyPID()==0) { // Only PE 0 does this section

    handle = fopen(filename,"w");
    if (!handle) return(-1);
    MM_typecode matcode;
    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_array(&matcode);
    mm_set_integer(&matcode);

    if (writeHeader==true) { // Only write header if requested (true by default)
    
      if (mm_write_banner(handle, matcode)) return(-1);
      
      if (mapName!=0) fprintf(handle, "%% \n%% %s\n", mapName);
      if (mapDescription!=0) fprintf(handle, "%% %s\n%% \n", mapDescription);

    }
  }
    
  if (writeHeader==true) { // Only write header if requested (true by default)

    // Make an Epetra_IntVector of length numProc such that all elements are on PE 0 and
    // the ith element is NumMyElements from the ith PE

    Epetra_Map map1(-1, 1, 0, map.Comm()); // map with one element on each processor
    int length = 0;
    if (map.Comm().MyPID()==0) length = map.Comm().NumProc();
    Epetra_Map map2(-1, length, 0, map.Comm());
    Epetra_Import lengthImporter(map2, map1);
    Epetra_IntVector v1(map1);
    Epetra_IntVector v2(map2);
    v1[0] = map.NumMyElements();
    if (v2.Import(v1, lengthImporter, Insert)) return(-1);
    if (map.Comm().MyPID()==0) { 
      fprintf(handle, "%s", "%Format Version:\n");
      //int version = 1; // We may change the format scheme at a later date.
      fprintf(handle, "%% %d \n", map.Comm().NumProc());
      fprintf(handle, "%s", "%NumProc: Number of processors:\n");
      fprintf(handle, "%% %d \n", map.Comm().NumProc());
      fprintf(handle, "%s", "%MaxElementSize: Maximum element size:\n");
      fprintf(handle, "%% %d \n", map.MaxElementSize());
      fprintf(handle, "%s", "%MinElementSize: Minimum element size:\n");
      fprintf(handle, "%% %d \n", map.MinElementSize());
      fprintf(handle, "%s", "%IndexBase: Index base of map:\n");
      fprintf(handle, "%% %d \n", map.IndexBase());
      fprintf(handle, "%s", "%NumGlobalElements: Total number of GIDs in map:\n");
      fprintf(handle, "%% %d \n", map.NumGlobalElements());
      fprintf(handle, "%s", "%NumMyElements: BlockMap lengths per processor:\n");
      for ( int i=0; i< v2.MyLength(); i++) fprintf(handle, "%% %d\n", v2[i]);
      
      if (mm_write_mtx_array_size(handle, M, N)) return(-1);
    }
  }
  if (BlockMapToHandle(handle, map)) return(-1); // Everybody calls this routine
  
  if (map.Comm().MyPID()==0) // Only PE 0 opened a file
    if (fclose(handle)) return(-1);
  return(0);
}

//.........这里部分代码省略.........
  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;
  if (colCount[myProc] > 0){
    colTotals = new int [colCount[myProc]];
    if (totalHEWeights > 0){
      colWeights = new double [colCount[myProc]];
    } 
  }
  int *colLocal= new int [ncols + 1];
  double *localWeights = NULL;
  if (totalHEWeights > 0){
    localWeights = new double [ncols + 1];
  }

  int base = colmap.IndexBase();
  int colStart = base;

  for (int i=0; i<nProcs; i++){

    // All processes send info to the process reponsible
    // for the next group of columns

    int ncols = colCount[i];
    int colEnd = colStart + ncols;
    for (int j=colStart,k=0; j < colEnd; j++,k++){
      gidIter = colGIDS.find(j);
      if (gidIter != colGIDS.end()){
        colLocal[k] = 1;     // column j has rows in my partition
      }
      else{
        colLocal[k] = 0;
      }
      if (totalHEWeights > 0){
        std::map<int, float>::iterator heWgtIter = hyperEdgeWeights.find(j);
        if (heWgtIter != hyperEdgeWeights.end()){
          // I have the edge weight for column j
          localWeights[k] = heWgtIter->second;
        }
        else{
          localWeights[k] = 0.0;
        }
      }
      
    }
#ifdef HAVE_MPI
    int rc = MPI_Reduce(colLocal, colTotals, ncols, MPI_INT, MPI_SUM, i, mcomm);
    if (totalHEWeights > 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;
//.........这里部分代码省略.........

//==============================================================================
bool Epetra_BlockMap::SameAs(const Epetra_BlockMap & Map) const {

  // Quickest test: See if both maps share an inner data class
  if (this->BlockMapData_ == Map.BlockMapData_) 
    return(true);

  if(!GlobalIndicesTypeMatch(Map))
    return(false);

  // Next check other global properties that are easy global attributes
  if (BlockMapData_->MinAllGID_ != Map.MinAllGID64() ||
      BlockMapData_->MaxAllGID_ != Map.MaxAllGID64() ||
      BlockMapData_->NumGlobalElements_ != Map.NumGlobalElements64() ||
      BlockMapData_->IndexBase_ != Map.IndexBase()) 
    return(false);
  
  // Last possible global check for constant element sizes
  if (BlockMapData_->ConstantElementSize_ && BlockMapData_->ElementSize_!=Map.ElementSize()) 
    return(false);

  // If we get this far, we need to check local properties and then check across
  // all processors to see if local properties are all true
 
  int numMyElements = BlockMapData_->NumMyElements_;

  int MySameMap = 1; // Assume not needed
  
  // First check if number of element is the same in each map
  if (numMyElements != Map.NumMyElements()) MySameMap = 0;
  
  // If numMyElements is the same, check to see that list of GIDs is the same
  if (MySameMap==1) {
    if (LinearMap() && Map.LinearMap() ) {
      // For linear maps, just need to check whether lower bound is the same
      if (MinMyGID64() != Map.MinMyGID64() )
        MySameMap = 0;
    }
    else {
      for (int i = 0; i < numMyElements; i++) {
        if (GID64(i) != Map.GID64(i)) {
          MySameMap = 0;
          break;
        }
      }
    }
  }
//    for (int i = 0; i < numMyElements; i++)
//      if (GID64(i) != Map.GID64(i)) MySameMap = 0;

  // If GIDs are the same, check to see element sizes are the same
  if (MySameMap==1 && !BlockMapData_->ConstantElementSize_) {
    int * sizeList1 = ElementSizeList();
    int * sizeList2 = Map.ElementSizeList();
    for (int i = 0; i < numMyElements; i++) if (sizeList1[i] != sizeList2[i]) MySameMap=0;
  }
  // Now get min of MySameMap across all processors

  int GlobalSameMap = 0;
  int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
  assert(err==0);

  return(GlobalSameMap==1);
}