C++ Epetra_Comm::Barrier方法代码示例

int fevec4(Epetra_Comm& Comm, bool verbose)
{
  int NumElements = 4;
  Epetra_Map     Map(NumElements, 0, Comm);
  Epetra_FEVector x1(Map);
  const double value = 1.;
  x1.PutScalar (value);
				// replace one element by itself. processor 0
				// does not own this element
  const int GID = 3;
  x1.ReplaceGlobalValues(1, &GID, &value);
  x1.GlobalAssemble (Insert);

  if (Map.MyGID(3)) {
    //insist that the value for GID==3 is 1:
    if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
  }

  std::cout << x1;

  Comm.Barrier();

				// re-apply GlobalAssemble. Nothing should
				// happen
  x1.GlobalAssemble (Insert);
  std::cout << x1;
  if (Map.MyGID(3)) {
    //insist that the value for GID==3 is 1:
    if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
  }

  return 0;
}

//============================================================================
void Ifpack_BreakForDebugger(Epetra_Comm& Comm)
{
  char hostname[80];
  char buf[80];
  if (Comm.MyPID()  == 0) cout << "Host and Process Ids for tasks" << endl;
  for (int i = 0; i <Comm.NumProc() ; i++) {
    if (i == Comm.MyPID() ) {
#if defined(TFLOP) || defined(JANUS_STLPORT)
      sprintf(buf, "Host: %s   PID: %d", "janus", getpid());
#elif defined(_WIN32)
      sprintf(buf,"Windows compiler, unknown hostname and PID!");
#else
      gethostname(hostname, sizeof(hostname));
      sprintf(buf, "Host: %s\tComm.MyPID(): %d\tPID: %d",
              hostname, Comm.MyPID(), getpid());
#endif
      printf("%s\n",buf);
      fflush(stdout);
#if !( defined(_WIN32) )
      sleep(1);
#endif
    }
  }
  if(Comm.MyPID() == 0) {
    printf("\n");
    printf("** Pausing to attach debugger...\n");
    printf("** You may now attach debugger to the processes listed above.\n");
    printf( "**\n");
    printf( "** Enter a character to continue > "); fflush(stdout);
    char go;
    scanf("%c",&go);
  }

  Comm.Barrier();

}

//
//  Amesos_TestMultiSolver.cpp reads in a matrix in Harwell-Boeing format, 
//  calls one of the sparse direct solvers, using blocked right hand sides
//  and computes the error and residual.  
//
//  TestSolver ignores the Harwell-Boeing right hand sides, creating
//  random right hand sides instead.  
//
//  Amesos_TestMultiSolver can test either A x = b or A^T x = b.
//  This can be a bit confusing because sparse direct solvers 
//  use compressed column storage - the transpose of Trilinos'
//  sparse row storage.
//
//  Matrices:
//    readA - Serial.  As read from the file.
//    transposeA - Serial.  The transpose of readA.
//    serialA - if (transpose) then transposeA else readA 
//    distributedA - readA distributed to all processes
//    passA - if ( distributed ) then distributedA else serialA
//
//
int Amesos_TestMultiSolver( Epetra_Comm &Comm, char *matrix_file, int numsolves, 
		      SparseSolverType SparseSolver, bool transpose,
		      int special, AMESOS_MatrixType matrix_type ) {


  int iam = Comm.MyPID() ;

  
  //  int hatever;
  //  if ( iam == 0 )  std::cin >> hatever ; 
  Comm.Barrier();


  Epetra_Map * readMap;
  Epetra_CrsMatrix * readA; 
  Epetra_Vector * readx; 
  Epetra_Vector * readb;
  Epetra_Vector * readxexact;
   
  std::string FileName = matrix_file ;
  int FN_Size = FileName.size() ; 
  std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
  std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
  bool NonContiguousMap = false; 

  if ( LastFiveBytes == ".triU" ) { 
    NonContiguousMap = true; 
    // Call routine to read in unsymmetric Triplet matrix
    EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx, 
						      readb, readxexact, NonContiguousMap ) );
  } else {
    if ( LastFiveBytes == ".triS" ) { 
      NonContiguousMap = true; 
      // Call routine to read in symmetric Triplet matrix
      EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, 
							readMap, readA, readx, 
							readb, readxexact, NonContiguousMap ) );
    } else {
      if (  LastFourBytes == ".mtx" ) { 
	EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap, 
							       readA, readx, readb, readxexact) );
      } else {
	// Call routine to read in HB problem
	Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx, 
						     readb, readxexact) ;
      }
    }
  }

  Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
  Epetra_CrsMatrix *serialA ; 

  if ( transpose ) {
    assert( CrsMatrixTranspose( readA, &transposeA ) == 0 ); 
    serialA = &transposeA ; 
  } else {
    serialA = readA ; 
  }

  // Create uniform distributed map
  Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
  Epetra_Map* map_;

  if( NonContiguousMap ) {
    //
    //  map gives us NumMyElements and MyFirstElement;
    //
    int NumGlobalElements =  readMap->NumGlobalElements();
    int NumMyElements = map.NumMyElements();
    int MyFirstElement = map.MinMyGID();
    std::vector<int> MapMap_( NumGlobalElements );
    readMap->MyGlobalElements( &MapMap_[0] ) ;
    Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ; 
    map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
  } else {
    map_ = new Epetra_Map( map ) ; 
  }


//.........这里部分代码省略.........

int MatrixMarketFileToBlockMap( const char *filename, const Epetra_Comm & comm, Epetra_BlockMap * & map) {

  const int lineLength = 1025;
  char line[lineLength];
  char token[lineLength];
  int M, N, numProc, MaxElementSize, MinElementSize, NumMyElements, IndexBase, NumGlobalElements, firstGid;

  FILE * handle = 0;

  bool inHeader = true;

  handle = fopen(filename,"r");
  if (handle == 0)
    EPETRA_CHK_ERR(-1); // file not found

  while (inHeader) {
    if(fgets(line, lineLength, handle)==0) return(-1);
    if(sscanf(line, "%s", token)==0) return(-1);
    if (!strcmp(token, "%NumProc:")) inHeader = false;
  }

  if(fgets(line, lineLength, handle)==0) return(-1); // numProc value
  if(sscanf(line, "%s %d", token, &numProc)==0) return(-1);

  if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize header line
  if(fgets(line, lineLength, handle)==0) return(-1); // MaxElementSize value
  if(sscanf(line, "%s %d", token, &MaxElementSize)==0) return(-1);

  if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize header line
  if(fgets(line, lineLength, handle)==0) return(-1); // MinElementSize value
  if(sscanf(line, "%s %d", token, &MinElementSize)==0) return(-1);

  if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase header line
  if(fgets(line, lineLength, handle)==0) return(-1); // IndexBase value
  if(sscanf(line, "%s %d", token, &IndexBase)==0) return(-1);

  if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements header line
  if(fgets(line, lineLength, handle)==0) return(-1); // NumGlobalElements value
  if(sscanf(line, "%s %d", token, &NumGlobalElements)==0) return(-1);

  int ierr = 0;
  if (comm.NumProc()==numProc) {
    if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
    firstGid = 0;
    for (int i=0; i<comm.MyPID(); i++) {
      if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value
      if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);
      firstGid += NumMyElements;
    }
 
    if(fgets(line, lineLength, handle)==0) return(-1); // This PE's NumMyElements value
    if(sscanf(line, "%s %d", token, &NumMyElements)==0) return(-1);

    for (int i=comm.MyPID()+1; i<numProc; i++) {
      if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
    }
  }
  else {
    ierr = 1; // Warning error, different number of processors.

    if(fgets(line, lineLength, handle)==0) return(-1); // NumMyElements header line
    for (int i=0; i<numProc; i++) {
      if(fgets(line, lineLength, handle)==0) return(-1); // ith NumMyElements value (dump these)
    }

    NumMyElements = NumGlobalElements/comm.NumProc();
    firstGid = comm.MyPID()*NumMyElements;
    int remainder = NumGlobalElements%comm.NumProc();
    if (comm.MyPID()<remainder) NumMyElements++;
    int extra = remainder;
    if (comm.MyPID()<remainder) extra = comm.MyPID();
    firstGid += extra;
  }
  if(fgets(line, lineLength, handle)==0) return(-1); // Number of rows, columns
  if(sscanf(line, "%d %d", &M, &N)==0) return(-1);

  bool doSizes = (N>1);
  Epetra_IntSerialDenseVector v1(NumMyElements);
  Epetra_IntSerialDenseVector v2(NumMyElements);
  for (int i=0; i<firstGid; i++) {
    if(fgets(line, lineLength, handle)==0) return(-1); // dump these
  }

  if (doSizes) {
    for (int i=0; i<NumMyElements; i++) {
      if(fgets(line, lineLength, handle)==0) return(-1);
      if(sscanf(line, "%d %d", &v1[i], &v2[i])==0) return(-1); // load v1, v2
    }
  }
  else {
    for (int i=0; i<NumMyElements; i++) {
      if(fgets(line, lineLength, handle)==0) return(-1);
      if(sscanf(line, "%d", &v1[i])==0) return(-1); // load v1
      v2[i] = MinElementSize; // Fill with constant size
    }
  }
  if (fclose(handle)) return(-1);

  comm.Barrier();

//.........这里部分代码省略.........

//------------------------------------------------------------------------------
int check_rowpermute_multivector_local(Epetra_Comm& Comm,
				       bool verbose)
{
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  Comm.Barrier();
  bool verbose1 = verbose;

  if (verbose) verbose = (MyPID==0);

  if (verbose) {
    cerr << "================check_rowpermute_multivector_local=========="
	 <<endl;
  }

  int NumMyElements = 5;
  int NumGlobalElements = NumMyElements*NumProc;
 
  Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);

  int* p = new int[NumMyElements];
  int firstGlobalRow = MyPID*NumMyElements;

  //Set up a permutation that will reverse the order of all LOCAL rows. (i.e.,
  //this test won't cause any inter-processor data movement.)

  if (verbose) {
    cout << "Permutation P:"<<endl;
  }

  int i;

  for(i=0; i<NumMyElements; ++i) {
    p[i] = firstGlobalRow+NumMyElements-1-i;
    if (verbose1) {
      cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
    }
  }

  Epetra_MultiVector v(Map, 3);

  double* v0 = v[0];
  double* v1 = v[1];
  double* v2 = v[2];

  for(i=0; i<NumMyElements; ++i) {
    v0[i] = 1.0*(firstGlobalRow+i) + 0.1;
    v1[i] = 1.0*(firstGlobalRow+i) + 0.2;
    v2[i] = 1.0*(firstGlobalRow+i) + 0.3;
  }
 
  if (verbose1) {
    cout << "*************** MultiVector v: ********************"<<endl;
    cout << v << endl;
  }

  EpetraExt::Permutation<Epetra_MultiVector> P(Copy, Map, p);

  Epetra_MultiVector& Pv = P(v);

  if (verbose1) {
    cout <<"************* permuted MultiVector Pv: ****************"<<endl;
    cout << Pv << endl;
  }

  return(0);
}

//-------------------------------------------------------------------------------
int check_colpermute_crsmatrix(Epetra_Comm& Comm,
			       bool verbose)
{
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  Comm.Barrier();
  bool verbose1 = verbose;

  if (verbose) verbose = (MyPID==0);

  if (verbose) {
    cerr << "================check_colpermute_crsmatrix=========="
	 <<endl;
  }

  int NumMyElements = 5;
  int NumGlobalElements = NumMyElements*NumProc;
 
  Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);

  int* p = new int[NumMyElements];
  int firstGlobalRow = MyPID*NumMyElements;

  if (verbose) {
    cout << "Permutation P:"<<endl;
  }

  int i;

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    p[i] = NumGlobalElements - row - 1;
    if (verbose1) {
      cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
    }
  }

  Epetra_CrsMatrix A(Copy, Map, 1);

  int col;
  double val;

  //set up a tri-diagonal graph.

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    col = NumGlobalElements - row - 1;
    val = 1.0*col;

    A.InsertGlobalValues(row, 1, &val, &col);

    if (col > 0) {
      int colm1 = col-1;
      val = 1.0*colm1;
      A.InsertGlobalValues(row, 1, &val, &colm1);
    }

    if (col < NumGlobalElements-1) {
      int colp1 = col+1;
      val = 1.0*colp1;
      A.InsertGlobalValues(row, 1, &val, &colp1);
    }
  }
 
  A.FillComplete();

  if (verbose1) {
    cout << "*************** matrix A: ********************"<<endl;
    cout << A << endl;
  }

  EpetraExt::Permutation<Epetra_CrsMatrix> P(Copy, Map, p);

  bool column_permutation = true;
  Epetra_CrsMatrix& B = P(A, column_permutation);

  if (verbose1) {
    cout <<"************* column-permuted matrix B: ****************"<<endl;
    cout << B << endl;
  }

  delete [] p;

  return(0);
}

//-------------------------------------------------------------------------------
int check_rowpermute_crsmatrix_global_diagonal(Epetra_Comm& Comm,
			bool verbose)
{
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  Comm.Barrier();
  bool verbose1 = verbose;

  if (verbose) verbose = (MyPID==0);

  if (verbose) {
    cerr << "================check_rowpermute_crsmatrix_global_diagonal=========="
	 <<endl;
  }

  int NumMyElements = 5;
  int NumGlobalElements = NumMyElements*NumProc;
 
  Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);

  int* p = new int[NumMyElements];
  int firstGlobalRow = MyPID*NumMyElements;

  //Now set up a permutation that will GLOBALLY reverse the order of all rows.
  //(i.e., if there are multiple processors, there will be inter-processor
  //data movement as rows are migrated.)

  int i;

  Epetra_CrsMatrix A(Copy, Map, 1);

  int col;
  double val;

  //set up a diagonal matrix A. It's diagonal because that's the easiest
  //to fill and to examine output before and after permutation...

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    val = 1.0*row;
    col = row;

    A.InsertGlobalValues(row, 1, &val, &col);
  }
 
  A.FillComplete();

  if (verbose1) {
    cout << "******************* matrix A: ****************************"<<endl;
    cout << A << endl;
  }

  if (verbose) {
    cout << "Permutation P:"<<endl;
  }

  for(i=0; i<NumMyElements; ++i) {
    int globalrow = NumGlobalElements-(firstGlobalRow+i)-1;
    p[i] = globalrow;
    if (verbose1) {
      cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
    }
  }

  EpetraExt::Permutation<Epetra_CrsMatrix> Pglobal(Copy, Map, p);

  Epetra_CrsMatrix& Bglobal = Pglobal(A);

  if (verbose1) {
    cout << "******************* permuted matrix Bglobal: *******************" <<endl;
    cout << Bglobal << endl;
  }

  return(0);
}

//------------------------------------------------------------------------------
int check_colpermute_crsgraph(Epetra_Comm& Comm,
			      bool verbose)
{
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  Comm.Barrier();
  bool verbose1 = verbose;

  if (verbose) verbose = (MyPID==0);

  if (verbose) {
    cerr << "================check_colpermute_crsgraph=========="
	 <<endl;
  }

  int NumMyElements = 5;
  int NumGlobalElements = NumMyElements*NumProc;
 
  Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);

  int* p = new int[NumMyElements];
  int firstGlobalRow = MyPID*NumMyElements;

  if (verbose) {
    cout << "Permutation P:"<<endl;
  }

  int i;

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    p[i] = NumGlobalElements - row - 1;
    if (verbose1) {
      cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
    }
  }

  Epetra_CrsGraph Agrph(Copy, Map, 1);

  int col;

  //set up a tri-diagonal graph.

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    col = NumGlobalElements - row - 1;

    Agrph.InsertGlobalIndices(row, 1, &col);

    if (col > 0) {
      int colm1 = col-1;
      Agrph.InsertGlobalIndices(row, 1, &colm1);
    }

    if (col < NumGlobalElements-1) {
      int colp1 = col+1;
      Agrph.InsertGlobalIndices(row, 1, &colp1);
    }
  }
 
  Agrph.FillComplete();

  if (verbose1) {
    cout << "*************** graph Agrph: ********************"<<endl;
    cout << Agrph << endl;
  }

  EpetraExt::Permutation<Epetra_CrsGraph> P(Copy, Map, p);

  bool column_permutation = true;
  Epetra_CrsGraph& Bgrph = P(Agrph, column_permutation);

  if (verbose1) {
    cout <<"************* column-permuted graph Bgrph: ****************"<<endl;
    cout << Bgrph << endl;
  }

  delete [] p;

  return(0);
}

//------------------------------------------------------------------------------
int check_rowpermute_crsgraph_local_diagonal(Epetra_Comm& Comm,
					  bool verbose)
{
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  Comm.Barrier();
  bool verbose1 = verbose;

  if (verbose) verbose = (MyPID==0);

  if (verbose) {
    cerr << "================check_rowpermute_crsgraph_local_diagonal=========="
	 <<endl;
  }

  int NumMyElements = 5;
  int NumGlobalElements = NumMyElements*NumProc;
 
  Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);

  int* p = new int[NumMyElements];
  int firstGlobalRow = MyPID*NumMyElements;

  //Set up a permutation that will reverse the order of all LOCAL rows. (i.e.,
  //this test won't cause any inter-processor data movement.)

  if (verbose) {
    cout << "Permutation P:"<<endl;
  }

  int i;

  for(i=0; i<NumMyElements; ++i) {
    p[i] = firstGlobalRow+NumMyElements-1-i;
    if (verbose1) {
      cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
    }
  }

  Epetra_CrsGraph Agrph(Copy, Map, 1);

  int col;

  //set up a diagonal graph. It's diagonal because that's the easiest
  //to fill and to examine output before and after permutation...

  for(i=0; i<NumMyElements; ++i) {
    int row = firstGlobalRow+i;
    col = row;

    Agrph.InsertGlobalIndices(row, 1, &col);
  }
 
  Agrph.FillComplete();

  if (verbose1) {
    cout << "*************** graph Agrph: ********************"<<endl;
    cout << Agrph << endl;
  }

  EpetraExt::Permutation<Epetra_CrsGraph> P(Copy, Map, p);

  Epetra_CrsGraph& Bgrph = P(Agrph);

  if (verbose1) {
    cout <<"************* permuted graph Bgrph: ****************"<<endl;
    cout << Bgrph << endl;
  }

  return(0);
}

//.........这里部分代码省略.........
  } else {
    map_ = new Epetra_Map( map ) ; 
  }


  Epetra_CrsMatrix A(Copy, *map_, 0);


  const Epetra_Map &OriginalMap = serialA->RowMatrixRowMap() ; 
  assert( OriginalMap.SameAs(*readMap) ); 
  Epetra_Export exporter(OriginalMap, *map_);
  Epetra_Export exporter2(OriginalMap, *map_);
  Epetra_Export MatrixExporter(OriginalMap, *map_);
  Epetra_CrsMatrix AwithDiag(Copy, *map_, 0);

  Epetra_Vector x(*map_);
  Epetra_Vector b(*map_);
  Epetra_Vector xexact(*map_);
  Epetra_Vector resid(*map_);
  Epetra_Vector readresid(*readMap);
  Epetra_Vector tmp(*map_);
  Epetra_Vector readtmp(*readMap);

  //  Epetra_Vector xcomp(*map_);      // X as computed by the solver
  bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ; 
  if ( distribute_matrix ) { 
    // Create Exporter to distribute read-in matrix and vectors
    //
    //  Initialize x, b and xexact to the values read in from the file
    //
    x.Export(*readx, exporter, Add);
    b.Export(*readb, exporter, Add);
    xexact.Export(*readxexact, exporter, Add);
    Comm.Barrier();
    
    A.Export(*serialA, exporter, Add);
    assert(A.FillComplete()==0);    
    
    Comm.Barrier();

    passA = &A; 

    passx = &x; 
    passb = &b;
    passxexact = &xexact;
    passresid = &resid;
    passtmp = &tmp;

  } else { 

    passA = serialA; 
    passx = readx; 
    passb = readb;
    passxexact = readxexact;
    passresid = &readresid;
    passtmp = &readtmp;
  }

  Epetra_MultiVector CopyB( *passb ) ;


  double Anorm = passA->NormInf() ; 
  SparseDirectTimingVars::SS_Result.Set_Anorm(Anorm) ;

  Epetra_LinearProblem Problem(  (Epetra_RowMatrix *) passA, 
				 (Epetra_MultiVector *) passx,