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

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

void Epetra_Export::Construct( const Epetra_BlockMap &  sourceMap, const Epetra_BlockMap & targetMap)
{

  int i;

  // Build three ID lists:
  // NumSameIDs - Number of IDs in TargetMap and SourceMap that are identical, up to the first
  //              nonidentical ID.
  // NumPermuteIDs - Number of IDs in SourceMap that must be indirectly loaded but are on this processor.
  // NumExportIDs - Number of IDs that are in SourceMap but not in TargetMap, and thus must be exported.

  int NumSourceIDs = sourceMap.NumMyElements();
  int NumTargetIDs = targetMap.NumMyElements();

  int_type *TargetGIDs = 0;
  if (NumTargetIDs>0) {
    TargetGIDs = new int_type[NumTargetIDs];
    targetMap.MyGlobalElements(TargetGIDs);
  }

  int_type * SourceGIDs = 0;
  if (NumSourceIDs>0) {
    SourceGIDs = new int_type[NumSourceIDs];
    sourceMap.MyGlobalElements(SourceGIDs);
  }

  int MinIDs = EPETRA_MIN(NumSourceIDs, NumTargetIDs);

  NumSameIDs_ = 0;
  for (i=0; i< MinIDs; i++) if (TargetGIDs[i]==SourceGIDs[i]) NumSameIDs_++; else break;

  // Find count of Source IDs that are truly remote and those that are local but permuted

  NumPermuteIDs_ = 0;
  NumExportIDs_ = 0;
  for (i=NumSameIDs_; i< NumSourceIDs; i++)
    if (targetMap.MyGID(SourceGIDs[i])) NumPermuteIDs_++; // Check if Source GID is a local Target GID
    else NumExportIDs_++; // If not, then it is remote

  // Define remote and permutation lists

  int_type * ExportGIDs = 0;
  if (NumExportIDs_>0) {
    ExportLIDs_ = new int[NumExportIDs_];
    ExportGIDs = new int_type[NumExportIDs_];
  }
  if (NumPermuteIDs_>0)  {
    PermuteToLIDs_ = new int[NumPermuteIDs_];
    PermuteFromLIDs_ = new int[NumPermuteIDs_];
  }

  NumPermuteIDs_ = 0;
  NumExportIDs_ = 0;
  for (i=NumSameIDs_; i< NumSourceIDs; i++) {
    if (targetMap.MyGID(SourceGIDs[i])) {
      PermuteFromLIDs_[NumPermuteIDs_] = i;
      PermuteToLIDs_[NumPermuteIDs_++] = targetMap.LID(SourceGIDs[i]);
    }
    else {
      //NumSend_ +=sourceMap.ElementSize(i); // Count total number of entries to send
      NumSend_ +=sourceMap.MaxElementSize(); // Count total number of entries to send (currently need max)
      ExportGIDs[NumExportIDs_] = SourceGIDs[i];
      ExportLIDs_[NumExportIDs_++] = i;
    }
  }

  if ( NumExportIDs_>0 && !sourceMap.DistributedGlobal())
    ReportError("Warning in Epetra_Export: Serial Export has remote IDs. (Exporting from Subset of Source Map)", 1);

  // Test for distributed cases
  int ierr = 0;

  if (sourceMap.DistributedGlobal()) {

    if (NumExportIDs_>0) ExportPIDs_ = new int[NumExportIDs_];
    ierr = targetMap.RemoteIDList(NumExportIDs_, ExportGIDs, ExportPIDs_, 0); // Get remote PIDs
    if( ierr ) throw ReportError("Error in Epetra_BlockMap::RemoteIDList", ierr);

    //Get rid of IDs not in Target Map
    if(NumExportIDs_>0) {
      int cnt = 0;
      for( i = 0; i < NumExportIDs_; ++i )
  if( ExportPIDs_[i] == -1 ) ++cnt;
      if( cnt ) {
  int_type * NewExportGIDs = 0;
  int * NewExportPIDs = 0;
  int * NewExportLIDs = 0;
  int cnt1 = NumExportIDs_-cnt;
  if (cnt1) {
    NewExportGIDs = new int_type[cnt1];
    NewExportPIDs = new int[cnt1];
    NewExportLIDs = new int[cnt1];
  }
  cnt = 0;
  for( i = 0; i < NumExportIDs_; ++i )
    if( ExportPIDs_[i] != -1 ) {
      NewExportGIDs[cnt] = ExportGIDs[i];
      NewExportPIDs[cnt] = ExportPIDs_[i];
      NewExportLIDs[cnt] = ExportLIDs_[i];
      ++cnt;
//.........这里部分代码省略.........

//==============================================================================
// Epetra_Export constructor for a Epetra_BlockMap object
Epetra_Export::Epetra_Export( const Epetra_BlockMap &  SourceMap, const Epetra_BlockMap & TargetMap)
  : Epetra_Object("Epetra::Export"), 
    TargetMap_(TargetMap),
    SourceMap_(SourceMap),
    NumSameIDs_(0),
    NumPermuteIDs_(0),
    PermuteToLIDs_(0),
    PermuteFromLIDs_(0),
    NumRemoteIDs_(0),
    RemoteLIDs_(0),
    NumExportIDs_(0),
    ExportLIDs_(0),
    ExportPIDs_(0),
    NumSend_(0),
    NumRecv_(0),
    Distor_(0)
{

  int i;

  // Build three ID lists:
  // NumSameIDs - Number of IDs in TargetMap and SourceMap that are identical, up to the first
  //              nonidentical ID.
  // NumPermuteIDs - Number of IDs in SourceMap that must be indirectly loaded but are on this processor.
  // NumExportIDs - Number of IDs that are in SourceMap but not in TargetMap, and thus must be exported.

  int NumSourceIDs = SourceMap.NumMyElements();
  int NumTargetIDs = TargetMap.NumMyElements();

  int *TargetGIDs = 0;
  if (NumTargetIDs>0) {
    TargetGIDs = new int[NumTargetIDs];
    TargetMap.MyGlobalElements(TargetGIDs);
  }

  int * SourceGIDs = 0;
  if (NumSourceIDs>0) {
    SourceGIDs = new int[NumSourceIDs];
    SourceMap.MyGlobalElements(SourceGIDs);
  }

  int MinIDs = EPETRA_MIN(NumSourceIDs, NumTargetIDs);

  NumSameIDs_ = 0;
  for (i=0; i< MinIDs; i++) if (TargetGIDs[i]==SourceGIDs[i]) NumSameIDs_++; else break;

  // Find count of Source IDs that are truly remote and those that are local but permuted

  NumPermuteIDs_ = 0;
  NumExportIDs_ = 0;
  for (i=NumSameIDs_; i< NumSourceIDs; i++) 
    if (TargetMap.MyGID(SourceGIDs[i])) NumPermuteIDs_++; // Check if Source GID is a local Target GID
    else NumExportIDs_++; // If not, then it is remote

  // Define remote and permutation lists

  int * ExportGIDs = 0;
  if (NumExportIDs_>0) {
    ExportLIDs_ = new int[NumExportIDs_];
    ExportGIDs = new int[NumExportIDs_];
  }
  if (NumPermuteIDs_>0)  {
    PermuteToLIDs_ = new int[NumPermuteIDs_];
    PermuteFromLIDs_ = new int[NumPermuteIDs_];
  }

  NumPermuteIDs_ = 0;
  NumExportIDs_ = 0;
  for (i=NumSameIDs_; i< NumSourceIDs; i++) {
    if (TargetMap.MyGID(SourceGIDs[i])) {
      PermuteFromLIDs_[NumPermuteIDs_] = i;
      PermuteToLIDs_[NumPermuteIDs_++] = TargetMap.LID(SourceGIDs[i]);
    }
    else {
      //NumSend_ +=SourceMap.ElementSize(i); // Count total number of entries to send
      NumSend_ +=SourceMap.MaxElementSize(); // Count total number of entries to send (currently need max)
      ExportGIDs[NumExportIDs_] = SourceGIDs[i];
      ExportLIDs_[NumExportIDs_++] = i;
    }
  }
     
  if ( NumExportIDs_>0 && !SourceMap.DistributedGlobal()) 
    ReportError("Warning in Epetra_Export: Serial Export has remote IDs. (Exporting from Subset of Source Map)", 1);

  // Test for distributed cases
  int ierr = 0;

  if (SourceMap.DistributedGlobal()) {

    if (NumExportIDs_>0) ExportPIDs_ = new int[NumExportIDs_];
    ierr = TargetMap.RemoteIDList(NumExportIDs_, ExportGIDs, ExportPIDs_, 0); // Get remote PIDs
    if( ierr ) throw ReportError("Error in Epetra_BlockMap::RemoteIDList", ierr);

    //Get rid of IDs not in Target Map
    if(NumExportIDs_>0) {
      int cnt = 0;
      for( i = 0; i < NumExportIDs_; ++i )
	if( ExportPIDs_[i] == -1 ) ++cnt;
//.........这里部分代码省略.........

void Epetra_Import::Construct( const Epetra_BlockMap &  targetMap, const Epetra_BlockMap & sourceMap)
{
  int i;
  
  // Build three ID lists:
  // NumSameIDs - Number of IDs in TargetMap and SourceMap that are identical, up to the first
  //              nonidentical ID.
  // NumPermuteIDs - Number of IDs in SourceMap that must be indirectly loaded but are on this processor.
  // NumRemoteIDs - Number of IDs that are in SourceMap but not in TargetMap, and thus must be imported.
  
  int NumSourceIDs = sourceMap.NumMyElements();
  int NumTargetIDs = targetMap.NumMyElements();
  
  int_type *TargetGIDs = 0;
  if (NumTargetIDs>0) {
    TargetGIDs = new int_type[NumTargetIDs];
    targetMap.MyGlobalElements(TargetGIDs);
  }
  
  int_type * SourceGIDs = 0;
  if (NumSourceIDs>0) {
    SourceGIDs = new int_type[NumSourceIDs];
    sourceMap.MyGlobalElements(SourceGIDs);
  }
  
  int MinIDs = EPETRA_MIN(NumSourceIDs, NumTargetIDs);
  
  
  NumSameIDs_ = 0;
  for (i=0; i< MinIDs; i++) if (TargetGIDs[i]==SourceGIDs[i]) NumSameIDs_++; else break;
  
  
  // Find count of Target IDs that are truly remote and those that are local but permuted

  NumPermuteIDs_ = 0;
  NumRemoteIDs_ = 0;
  for (i=NumSameIDs_; i< NumTargetIDs; i++) 
    if (sourceMap.MyGID(TargetGIDs[i])) NumPermuteIDs_++; // Check if Target GID is a local Source GID
    else NumRemoteIDs_++; // If not, then it is remote
  
  
  
  // Define remote and permutation lists
  
  int_type * RemoteGIDs=0;
  RemoteLIDs_ = 0;
  if (NumRemoteIDs_>0) {
    RemoteLIDs_ = new int[NumRemoteIDs_];
    RemoteGIDs = new int_type[NumRemoteIDs_];
  }
  if (NumPermuteIDs_>0)  {
    PermuteToLIDs_ = new int[NumPermuteIDs_];
    PermuteFromLIDs_ = new int[NumPermuteIDs_];
  }
  
  NumPermuteIDs_ = 0;
  NumRemoteIDs_ = 0;
  for (i=NumSameIDs_; i< NumTargetIDs; i++) {
    if (sourceMap.MyGID(TargetGIDs[i])) {
      PermuteToLIDs_[NumPermuteIDs_] = i;
      PermuteFromLIDs_[NumPermuteIDs_++] = sourceMap.LID(TargetGIDs[i]);
    }
    else {
      //NumRecv_ +=TargetMap.ElementSize(i); // Count total number of entries to receive
      NumRecv_ +=targetMap.MaxElementSize(); // Count total number of entries to receive (currently need max)
      RemoteGIDs[NumRemoteIDs_] = TargetGIDs[i];
      RemoteLIDs_[NumRemoteIDs_++] = i;
    }
  }

  if( NumRemoteIDs_>0 && !sourceMap.DistributedGlobal() )
    ReportError("Warning in Epetra_Import: Serial Import has remote IDs. (Importing to Subset of Target Map)", 1);
  
  // Test for distributed cases
  
  int * RemotePIDs = 0;

  if (sourceMap.DistributedGlobal()) {
    
    if (NumRemoteIDs_>0)  RemotePIDs = new int[NumRemoteIDs_];
    int ierr = sourceMap.RemoteIDList(NumRemoteIDs_, RemoteGIDs, RemotePIDs, 0); // Get remote PIDs
    if (ierr) throw ReportError("Error in sourceMap.RemoteIDList call", ierr);

    //Get rid of IDs that don't exist in SourceMap
    if(NumRemoteIDs_>0) {
      int cnt = 0;
      for( i = 0; i < NumRemoteIDs_; ++i )
        if( RemotePIDs[i] == -1 ) ++cnt;
      if( cnt ) {
        if( NumRemoteIDs_-cnt ) {
          int_type * NewRemoteGIDs = new int_type[NumRemoteIDs_-cnt];
          int * NewRemotePIDs = new int[NumRemoteIDs_-cnt];
          int * NewRemoteLIDs = new int[NumRemoteIDs_-cnt];
          cnt = 0;
          for( i = 0; i < NumRemoteIDs_; ++i )
            if( RemotePIDs[i] != -1 ) {
              NewRemoteGIDs[cnt] = RemoteGIDs[i];
              NewRemotePIDs[cnt] = RemotePIDs[i];
              NewRemoteLIDs[cnt] = targetMap.LID(RemoteGIDs[i]);
              ++cnt;
//.........这里部分代码省略.........

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

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;
//.........这里部分代码省略.........