C++ Epetra_RowMatrix::Comm方法代码示例

// ======================================================================
int GetGlobalAggregates(Epetra_RowMatrix& A, Teuchos::ParameterList& List,
                        double* thisns, Epetra_IntVector& aggrinfo)
{
  int naggregates = GetAggregates(A,List,thisns,aggrinfo);
  const Epetra_Comm& comm = A.Comm();
  std::vector<int> local(comm.NumProc());
  std::vector<int> global(comm.NumProc());
  for (int i=0; i<comm.NumProc(); ++i) local[i] = 0;
  local[comm.MyPID()] = naggregates;
  comm.SumAll(&local[0],&global[0],comm.NumProc());
  int offset = 0;
  for (int i=0; i<comm.MyPID(); ++i) offset += global[i];
  for (int i=0; i<aggrinfo.MyLength(); ++i)
    if (aggrinfo[i]<naggregates) aggrinfo[i] += offset;
    else                         aggrinfo[i] = -1;
  return naggregates;
}

double Ifpack_FrobeniusNorm(const Epetra_RowMatrix& A)
{
  double MyNorm = 0.0, GlobalNorm;

  std::vector<int> colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  for (int i = 0 ; i < A.NumMyRows() ; ++i) {

    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    for (int j = 0 ; j < Nnz ; ++j) {
      MyNorm += colVal[j] * colVal[j];
    }
  }

  A.Comm().SumAll(&MyNorm,&GlobalNorm,1);

  return(sqrt(GlobalNorm));
}

int TestMultiLevelPreconditioner(char ProblemType[],
				 Teuchos::ParameterList & MLList,
				 Epetra_LinearProblem & Problem, double & TotalErrorResidual,
				 double & TotalErrorExactSol,bool cg=false)
{

  Epetra_MultiVector* lhs = Problem.GetLHS();
  Epetra_MultiVector* rhs = Problem.GetRHS();
  Epetra_RowMatrix* A = Problem.GetMatrix();

  // ======================================== //
  // create a rhs corresponding to lhs or 1's //
  // ======================================== //

  lhs->PutScalar(1.0);
  A->Multiply(false,*lhs,*rhs);

  lhs->PutScalar(0.0);

  Epetra_Time Time(A->Comm());

  // =================== //
  // call ML and AztecOO //
  // =================== //

  AztecOO solver(Problem);

  MLList.set("ML output", 10);

  ML_Epetra::MultiLevelPreconditioner * MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList, true);

  // tell AztecOO to use this preconditioner, then solve
  solver.SetPrecOperator(MLPrec);

  if(cg) solver.SetAztecOption(AZ_solver, AZ_cg);
  else solver.SetAztecOption(AZ_solver, AZ_gmres);
  solver.SetAztecOption(AZ_output, 32);
  solver.SetAztecOption(AZ_kspace, 160);

  solver.Iterate(1550, 1e-12);

  delete MLPrec;

  // ==================================================== //
  // compute difference between exact solution and ML one //
  // ==================================================== //

  double d = 0.0, d_tot = 0.0;

  for( int i=0 ; i<lhs->Map().NumMyElements() ; ++i )
    d += ((*lhs)[0][i] - 1.0) * ((*lhs)[0][i] - 1.0);

  A->Comm().SumAll(&d,&d_tot,1);

  // ================== //
  // compute ||Ax - b|| //
  // ================== //

  double Norm;
  Epetra_Vector Ax(rhs->Map());
  A->Multiply(false, *lhs, Ax);
  Ax.Update(1.0, *rhs, -1.0);
  Ax.Norm2(&Norm);

  string msg = ProblemType;

  if (A->Comm().MyPID() == 0) {
    cout << msg << "......Using " << A->Comm().NumProc() << " processes" << endl;
    cout << msg << "......||A x - b||_2 = " << Norm << endl;
    cout << msg << "......||x_exact - x||_2 = " << sqrt(d_tot) << endl;
    cout << msg << "......Total Time = " << Time.ElapsedTime() << endl;
  }

  TotalErrorExactSol += sqrt(d_tot);
  TotalErrorResidual += Norm;

  return( solver.NumIters() );

}

int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {

  int ierr = 0;
  EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
  EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
  EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
  EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
  EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
  EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
  EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalCols64()==B.NumGlobalCols64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalDiagonals64()==B.NumGlobalDiagonals64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalNonzeros64()==B.NumGlobalNonzeros64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalRows64()==B.NumGlobalRows64(),ierr);
  EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
  EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
  EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
  for (int i=0; i<A.NumMyRows(); i++) {
    int nA, nB;
    A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
    EPETRA_TEST_ERR(!nA==nB,ierr);
  }
  EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
  EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
  EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
  EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
  EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
  EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
  EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);

  int NumVectors = 5;
  { // No transpose case
    Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
    Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
    Epetra_MultiVector YA2(YA1);
    Epetra_MultiVector YB1(YA1);
    Epetra_MultiVector YB2(YA1);
    X.Random();

    bool transA = false;
    A.SetUseTranspose(transA);
    B.SetUseTranspose(transA);
    A.Apply(X,YA1);
    A.Multiply(transA, X, YA2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
    B.Apply(X,YB1);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
    B.Multiply(transA, X, YB2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);

  }
  {// transpose case
    Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
    Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
    Epetra_MultiVector YA2(YA1);
    Epetra_MultiVector YB1(YA1);
    Epetra_MultiVector YB2(YA1);
    X.Random();

    bool transA = true;
    A.SetUseTranspose(transA);
    B.SetUseTranspose(transA);
    A.Apply(X,YA1);
    A.Multiply(transA, X, YA2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
    B.Apply(X,YB1);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
    B.Multiply(transA, X,YB2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);

  }

  Epetra_Vector diagA(A.RowMatrixRowMap());
  EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
  Epetra_Vector diagB(B.RowMatrixRowMap());
  EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
  EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);

  Epetra_Vector rowA(A.RowMatrixRowMap());
  EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
  Epetra_Vector rowB(B.RowMatrixRowMap());
  EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
  EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);

  Epetra_Vector colA(A.RowMatrixColMap());
  EPETRA_TEST_ERR(A.InvColSums(colA),ierr);
  Epetra_Vector colB(B.RowMatrixColMap());
  EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
  EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr);

  EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
  EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);

  EPETRA_TEST_ERR(A.RightScale(colA),ierr);
  EPETRA_TEST_ERR(B.RightScale(colB),ierr);


  EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
  EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);

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

int checkResults( bool trans, 
		  Epetra_LinearProblemRedistor * redistor, 
		  Epetra_LinearProblem * A,
		  Epetra_LinearProblem * R,
		  bool verbose) {
  
  int m = A->GetRHS()->MyLength();
  int n = A->GetLHS()->MyLength();
  assert( m == n ) ; 

  Epetra_MultiVector *x = A->GetLHS() ; 
  Epetra_MultiVector x1( *x ) ; 
  //  Epetra_MultiVector Difference( x1 ) ; 
  Epetra_MultiVector *b = A->GetRHS();
  Epetra_RowMatrix *matrixA = A->GetMatrix();
  assert( matrixA != 0 ) ; 
  int iam = matrixA->Comm().MyPID();
  
  //  Epetra_Time timer(A->Comm());
  //  double start = timer.ElapsedTime();

  matrixA->Multiply(trans, *b, x1) ;   // x = Ab 

  int M,N,nz;
  int *ptr, *ind;
  double *val, *rhs, *lhs;
  int Nrhs, ldrhs, ldlhs;

  redistor->ExtractHbData( M, N, nz, ptr, ind, 
		    val, Nrhs, rhs, ldrhs, 
		    lhs, ldlhs);

  assert( M == N ) ; 
  if ( verbose ) {
    cout << " iam = " << iam 
	 << " m = " << m  << " n = " << n  << " M = " << M << endl ;  
    
    cout << " iam = " << iam << " ptr = " << ptr[0] << " "   << ptr[1] << " "  << ptr[2] << " "  << ptr[3] << " "  << ptr[4] << " " << ptr[5] << endl ;
    
    cout << " iam = " << iam << " ind = " << ind[0] << " "   << ind[1] << " "  << ind[2] << " "  << ind[3] << " "  << ind[4] << " " << ind[5] << endl ;
    
    cout << " iam = " << iam << " val = " << val[0] << " "   << val[1] << " "  << val[2] << " "  << val[3] << " "  << val[4] << " " << val[5] << endl ;
  }
  //  Create a serial map in case we end up needing it 
  //  If it is created inside the else block below it would have to
  //  be with a call to new().
  int NumMyElements_ = 0 ;
  if (matrixA->Comm().MyPID()==0) NumMyElements_ = n;
  Epetra_Map SerialMap( n, NumMyElements_, 0, matrixA->Comm() );

  //  These are unnecessary and useless
  //  Epetra_Vector serial_A_rhs( SerialMap ) ; 
  //  Epetra_Vector serial_A_lhs( SerialMap ) ;

  //  Epetra_Export exporter( matrixA->BlockRowMap(), SerialMap) ;

  //
  //  In each process, we will compute Rb putting the answer into LHS
  //


  for ( int k = 0 ; k < Nrhs; k ++ ) { 
    for ( int i = 0 ; i < M ; i ++ ) { 
      lhs[ i + k * ldlhs ] = 0.0; 
    }
    for ( int i = 0 ; i < M ; i++ ) { 
      for ( int l = ptr[i]; l < ptr[i+1]; l++ ) {
	int j = ind[l] ; 
	if ( verbose && N < 40 ) {
	  cout << " i = " << i << " j = " << j ;
	  cout << " l = " << l << " val[l] = " << val[l] ;
	  cout << " rhs = " << rhs[ j + k * ldrhs ] << endl ;
	}
	lhs[ i + k * ldrhs ] += val[l] * rhs[ j + k * ldrhs ] ;
      }
    }

    if ( verbose && N < 40 ) { 
      cout << " lhs = " ; 
      for ( int j = 0 ; j < N ; j++ ) cout << " " << lhs[j] ; 
      cout << endl ; 
      cout << " rhs = " ; 
      for ( int j = 0 ; j < N ; j++ ) cout << " " << rhs[j] ; 
      cout << endl ; 
    }

    const Epetra_Comm &comm = matrixA->Comm() ; 
#ifdef HAVE_COMM_ASSERT_EQUAL
    //
    //  Here we double check to make sure that lhs and rhs are 
    //  replicated.  
    //
    for ( int j = 0 ; j < N ; j++ ) { 
      assert( Comm_assert_equal( &comm, lhs[ j + k * ldrhs ] ) ) ; 
      assert( Comm_assert_equal( &comm, rhs[ j + k * ldrhs ] ) ) ; 
    } 
#endif
  }

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

// ======================================================================
void GetPtent(const Epetra_RowMatrix& A, Teuchos::ParameterList& List,
              double* thisns, Teuchos::RCP<Epetra_CrsMatrix>& Ptent,
              Teuchos::RCP<Epetra_MultiVector>& NextNS, const int domainoffset)
{
  const int nsdim = List.get<int>("null space: dimension",-1);
  if (nsdim<=0) ML_THROW("null space dimension not given",-1);
  const Epetra_Map& rowmap = A.RowMatrixRowMap();
  const int mylength = rowmap.NumMyElements();

  // wrap nullspace into something more handy
  Epetra_MultiVector ns(View,rowmap,thisns,mylength,nsdim);

  // vector to hold aggregation information
  Epetra_IntVector aggs(rowmap,false);

  // aggregation with global aggregate numbering
  int naggregates = GetGlobalAggregates(const_cast<Epetra_RowMatrix&>(A),List,thisns,aggs);

  // build a domain map for Ptent
  // find first aggregate on proc
  int firstagg = -1;
  int offset = -1;
  for (int i=0; i<mylength; ++i)
    if (aggs[i]>=0)
    {
      offset = firstagg = aggs[i];
      break;
    }
  offset *= nsdim;
  if (offset<0) ML_THROW("could not find any aggregate on proc",-2);
  std::vector<int> coarsegids(naggregates*nsdim);
  for (int i=0; i<naggregates; ++i)
    for (int j=0; j<nsdim; ++j)
    {
      coarsegids[i*nsdim+j] = offset + domainoffset;
      ++offset;
    }
  Epetra_Map pdomainmap(-1,naggregates*nsdim,&coarsegids[0],0,A.Comm());

  // loop aggregates and build ids for dofs
  std::map<int,std::vector<int> > aggdofs;
  std::map<int,std::vector<int> >::iterator fool;
  for (int i=0; i<naggregates; ++i)
  {
    std::vector<int> gids(0);
    aggdofs.insert(std::pair<int,std::vector<int> >(firstagg+i,gids));
  }
  for (int i=0; i<mylength; ++i)
  {
    if (aggs[i]<0) continue;
    std::vector<int>& gids = aggdofs[aggs[i]];
    gids.push_back(aggs.Map().GID(i));
  }

#if 0 // debugging output
  for (int proc=0; proc<A.Comm().NumProc(); ++proc)
  {
    if (proc==A.Comm().MyPID())
    {
      for (fool=aggdofs.begin(); fool!=aggdofs.end(); ++fool)
      {
        std::cout << "Proc " << proc << " Aggregate " << fool->first << " Dofs ";
        std::vector<int>& gids = fool->second;
        for (int i=0; i<(int)gids.size(); ++i) std::cout << gids[i] << " ";
        std::cout << std::endl;
      }
    }
    fflush(stdout);
    A.Comm().Barrier();
  }
#endif

  // coarse level nullspace to be filled
  NextNS = Teuchos::rcp(new Epetra_MultiVector(pdomainmap,nsdim,true));
  Epetra_MultiVector& nextns = *NextNS;

  // matrix Ptent
  Ptent = Teuchos::rcp(new Epetra_CrsMatrix(Copy,rowmap,nsdim));

  // loop aggregates and extract the appropriate slices of the nullspace.
  // do QR and assemble Q and R to Ptent and NextNS
  for (fool=aggdofs.begin(); fool!=aggdofs.end(); ++fool)
  {
    // extract aggregate-local junk of nullspace
    const int aggsize = (int)fool->second.size();
    Epetra_SerialDenseMatrix Bagg(aggsize,nsdim);
    for (int i=0; i<aggsize; ++i)
      for (int j=0; j<nsdim; ++j)
        Bagg(i,j) = (*ns(j))[ns.Map().LID(fool->second[i])];

    // Bagg = Q*R
    int m = Bagg.M();
    int n = Bagg.N();
    int lwork = n*10;
    int info = 0;
    int k = std::min(m,n);
    if (k!=n) ML_THROW("Aggregate too small, fatal!",-1);
    std::vector<double> work(lwork);
    std::vector<double> tau(k);
//.........这里部分代码省略.........

// ======================================================================
int Ifpack_PrintSparsity(const Epetra_RowMatrix& A, const char* InputFileName,
                         const int NumPDEEqns)
{

  int ltit;
  long long m,nc,nr,maxdim;
  double lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz = 0.0;
  double xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx;
  bool square = false;
  /*change square to .true. if you prefer a square frame around
    a rectangular matrix */
  double conv = 2.54;
  char munt = 'E'; /* put 'E' for centimeters, 'U' for inches */
  int ptitle = 0; /* position of the title, 0 under the drawing,
                     else above */
  FILE* fp = NULL;
  int NumMyRows;
  //int NumMyCols;
  long long NumGlobalRows;
  long long NumGlobalCols;
  int MyPID;
  int NumProc;
  char FileName[1024];
  char title[1024];

  const Epetra_Comm& Comm = A.Comm();

  /* --------------------- execution begins ---------------------- */

  if (strlen(A.Label()) != 0)
    strcpy(title, A.Label());
  else
    sprintf(title, "%s", "matrix");

  if (InputFileName == 0)
    sprintf(FileName, "%s.ps", title);
  else
    strcpy(FileName, InputFileName);

  MyPID = Comm.MyPID();
  NumProc = Comm.NumProc();

  NumMyRows = A.NumMyRows();
  //NumMyCols = A.NumMyCols();

  NumGlobalRows = A.NumGlobalRows64();
  NumGlobalCols = A.NumGlobalCols64();

  if (NumGlobalRows != NumGlobalCols)
    IFPACK_CHK_ERR(-1); // never tested

  /* to be changed for rect matrices */
  maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
  maxdim /= NumPDEEqns;

  m = 1 + maxdim;
  nr = NumGlobalRows / NumPDEEqns + 1;
  nc = NumGlobalCols / NumPDEEqns + 1;

  if (munt == 'E') {
    u2dot = 72.0/conv;
    paperx = 21.0;
    siz = 10.0;
  }
  else {
    u2dot = 72.0;
    paperx = 8.5*conv;
    siz = siz*conv;
  }

  /* left and right margins (drawing is centered) */

  lrmrgn = (paperx-siz)/2.0;

  /* bottom margin : 2 cm */

  botmrgn = 2.0;
  /* c scaling factor */
  scfct = siz*u2dot/m;
  /* matrix frame line witdh */
  frlw = 0.25;
  /* font size for title (cm) */
  fnstit = 0.5;
  /* mfh 23 Jan 2013: title is always nonnull, since it's an array of
     fixed nonzero length.  The 'if' test thus results in a compiler
     warning. */
  /*if (title) ltit = strlen(title);*/
  /*else       ltit = 0;*/
  ltit = strlen(title);

  /* position of title : centered horizontally */
  /*                     at 1.0 cm vertically over the drawing */
  ytitof = 1.0;
  xtit = paperx/2.0;
  ytit = botmrgn+siz*nr/m + ytitof;
  /* almost exact bounding box */
  xl = lrmrgn*u2dot - scfct*frlw/2;
  xr = (lrmrgn+siz)*u2dot + scfct*frlw/2;
  yb = botmrgn*u2dot - scfct*frlw/2;
//.........这里部分代码省略.........

int Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix& A,
                                 const bool abs, const int steps)
{

  bool verbose = (A.Comm().MyPID() == 0);
  double min_val =  DBL_MAX;
  double max_val = -DBL_MAX;

  std::vector<int>    colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  for (int i = 0 ; i < A.NumMyRows() ; ++i) {

    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    for (int j = 0 ; j < Nnz ; ++j) {
      double v = colVal[j];
      if (abs)
        if (v < 0) v = -v;
      if (v < min_val)
        min_val = v;
      if (v > max_val)
        max_val = v;
    }
  }

  if (verbose) {
    cout << endl;
    Ifpack_PrintLine();
    cout << "Label of matrix = " << A.Label() << endl;
    cout << endl;
  }

  double delta = (max_val - min_val) / steps;
  for (int k = 0 ; k < steps ; ++k) {

    double below = delta * k + min_val;
    double above = below + delta;
    int MyBelow = 0, GlobalBelow;

    for (int i = 0 ; i < A.NumMyRows() ; ++i) {

      int Nnz;
      IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                        &colVal[0],&colInd[0]));

      for (int j = 0 ; j < Nnz ; ++j) {
        double v = colVal[j];
        if (abs)
          if (v < 0) v = -v;
        if (v >= below && v < above) MyBelow++;
      }

    }
    A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
    if (verbose) {
      printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
             below, above, GlobalBelow,
             100.0 * GlobalBelow / A.NumGlobalNonzeros64());
    }
  }

  if (verbose) {
    Ifpack_PrintLine();
    cout << endl;
  }

  return(0);
}

int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
                   const int NumPDEEqns)
{

  int NumMyRows = A.NumMyRows();
  long long NumGlobalRows = A.NumGlobalRows64();
  long long NumGlobalCols = A.NumGlobalCols64();
  long long MyBandwidth = 0, GlobalBandwidth;
  long long MyLowerNonzeros = 0, MyUpperNonzeros = 0;
  long long GlobalLowerNonzeros, GlobalUpperNonzeros;
  long long MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
  long long MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
  double MyMin, MyAvg, MyMax;
  double GlobalMin, GlobalAvg, GlobalMax;
  long long GlobalStorage;

  bool verbose = (A.Comm().MyPID() == 0);

  GlobalStorage = sizeof(int*) * NumGlobalRows +
    sizeof(int) * A.NumGlobalNonzeros64() +
    sizeof(double) * A.NumGlobalNonzeros64();

  if (verbose) {
    print();
    Ifpack_PrintLine();
    print<const char*>("Label", A.Label());
    print<long long>("Global rows", NumGlobalRows);
    print<long long>("Global columns", NumGlobalCols);
    print<long long>("Stored nonzeros", A.NumGlobalNonzeros64());
    print<long long>("Nonzeros / row", A.NumGlobalNonzeros64() / NumGlobalRows);
    print<double>("Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
  }

  long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
  long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
  long long NumMyDirichletRows = 0, NumGlobalDirichletRows;

  std::vector<int> colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  Epetra_Vector Diag(A.RowMatrixRowMap());
  Epetra_Vector RowSum(A.RowMatrixRowMap());
  Diag.PutScalar(0.0);
  RowSum.PutScalar(0.0);

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

    long long GRID = A.RowMatrixRowMap().GID64(i);
    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    if (Nnz == 0)
      NumMyEmptyRows++;

    if (Nnz == 1)
      NumMyDirichletRows++;

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

      double v = colVal[j];
      if (v < 0) v = -v;
      if (colVal[j] != 0.0)
        NumMyActualNonzeros++;

      long long GCID = A.RowMatrixColMap().GID64(colInd[j]);

      if (GCID != GRID)
        RowSum[i] += v;
      else
        Diag[i] = v;

      if (GCID < GRID)
        MyLowerNonzeros++;
      else if (GCID > GRID)
        MyUpperNonzeros++;
      long long b = GCID - GRID;
      if (b < 0) b = -b;
      if (b > MyBandwidth)
        MyBandwidth = b;
    }

    if (Diag[i] > RowSum[i])
      MyDiagonallyDominant++;

    if (Diag[i] >= RowSum[i])
      MyWeaklyDiagonallyDominant++;

    RowSum[i] += Diag[i];
  }

  // ======================== //
  // summing up global values //
  // ======================== //

  A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
  A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
  A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
  A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
  A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
//.........这里部分代码省略.........

int TestMultiLevelPreconditioner(char ProblemType[],
				 Teuchos::ParameterList & MLList,
				 Epetra_LinearProblem & Problem, double & TotalErrorResidual,
				 double & TotalErrorExactSol)
{
  
  Epetra_MultiVector* lhs = Problem.GetLHS();
  Epetra_MultiVector* rhs = Problem.GetRHS();
  Epetra_RowMatrix* A = Problem.GetMatrix();
  
  // ======================================== //
  // create a rhs corresponding to lhs or 1's //
  // ======================================== //
  
  lhs->PutScalar(1.0);
  A->Multiply(false,*lhs,*rhs);

  lhs->PutScalar(0.0);
  
  Epetra_Time Time(A->Comm());

  Epetra_MultiVector lhs2(*lhs);
  Epetra_MultiVector rhs2(*rhs);
  
  // =================== //
  // call ML and AztecOO //
  // =================== //

  AztecOO solver(Problem);
  
  MLList.set("ML output", 0);
  ML_set_random_seed(24601);
  ML_Epetra::MultiLevelPreconditioner * MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList, true);
  
  // tell AztecOO to use this preconditioner, then solve
  solver.SetPrecOperator(MLPrec);
  
  solver.SetAztecOption(AZ_solver, AZ_gmres);
  solver.SetAztecOption(AZ_output, 32);
  solver.SetAztecOption(AZ_kspace, 160);
  
  solver.Iterate(1550, 1e-12);
  
  delete MLPrec;



  // ================================= //
  // call ML and AztecOO a second time //
  // ================================= // 
  Epetra_LinearProblem Problem2(A,&lhs2,&rhs2);

  AztecOO solver2(Problem2);
  ML_set_random_seed(24601);  
  ML_Epetra::MultiLevelPreconditioner * MLPrec2 = new ML_Epetra::MultiLevelPreconditioner(*A, MLList, true);
  
  // tell AztecOO to use this preconditioner, then solve
  solver2.SetPrecOperator(MLPrec2);
  
  solver2.SetAztecOption(AZ_solver, AZ_gmres);
  solver2.SetAztecOption(AZ_output, 32);
  solver2.SetAztecOption(AZ_kspace, 160);
  
  solver2.Iterate(1550, 1e-12);
  
  
  
  // ==================================================== //
  // compute difference between the two ML solutions //
  // ==================================================== //
  
  double d = 0.0, d_tot = 0.0;
  
  for( int i=0 ; i<lhs->Map().NumMyElements() ; ++i )
    d += ((*lhs)[0][i] - lhs2[0][i]) * ((*lhs)[0][i] - lhs2[0][i]);
  
  A->Comm().SumAll(&d,&d_tot,1);
  string msg = ProblemType;
  if (A->Comm().MyPID() == 0) {
    cout << msg << "......Using " << A->Comm().NumProc() << " processes" << endl;
    cout << msg << "......||x_1 - x_2||_2 = " << sqrt(d_tot) << endl;
    cout << msg << "......Total Time = " << Time.ElapsedTime() << endl;
  }
  
  TotalErrorExactSol += sqrt(d_tot);
  
  return( solver.NumIters() );
  
}