C++ Epetra_Vector::Update方法代码示例

/*!
 * \brief Solve.
 */
void JacobiSolver::iterate( const int max_iters, const double tolerance )
{
    // Extract the linear problem.
    Epetra_CrsMatrix *A = 
	dynamic_cast<Epetra_CrsMatrix*>( d_linear_problem->GetMatrix() );
    Epetra_Vector *x = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetLHS() );
    const Epetra_Vector *b = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetRHS() );

    // Setup the residual.
    Epetra_Map row_map = A->RowMap();
    Epetra_Vector residual( row_map );

    // Iterate.
    Epetra_CrsMatrix H = buildH( A );
    Epetra_Vector temp_vec( row_map );
    d_num_iters = 0;
    double residual_norm = 1.0;
    double b_norm;
    b->NormInf( &b_norm );
    double conv_crit = b_norm*tolerance;
    while ( residual_norm > conv_crit && d_num_iters < max_iters )
    {
	H.Apply( *x, temp_vec );
	x->Update( 1.0, temp_vec, 1.0, *b, 0.0 );

	A->Apply( *x, temp_vec );
	residual.Update( 1.0, *b, -1.0, temp_vec, 0.0 );

	residual.NormInf( &residual_norm );
	++d_num_iters;
    }
}

int power_method(bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0,
		 Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
{
	
  // Fill z with random Numbers
  Epetra_Vector z(z0);
	
  // variable needed for iteration
  double normz, residual;

  int ierr = 1;
	
  for(int iter = 0; iter < niters; iter++) {
		z.Norm2(&normz); // Compute 2-norm of z
		q.Scale(1.0/normz, z);
		A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
		q.Dot(z, lambda); // Approximate maximum eigenvaluE
		if(iter%100==0 || iter+1==niters) {
			resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
			resid.Norm2(&residual);
			if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda
											 << "  Residual of A*q - lambda*q = " << residual << endl;
		}
		if(residual < tolerance) {
			ierr = 0;
			break;
		}
	}
  return(ierr);
}

/*----------------------------------------------------------------------*
 |  evaluate nonlinear function (public, derived)             m.gee 3/06|
 *----------------------------------------------------------------------*/
bool NLNML::NLNML_CoarseLevelNoxInterface::computeF(
                                 const Epetra_Vector& x, Epetra_Vector& F,
               const FillType fillFlag)
{
  bool err;
  if (!Level())
  {
    err = fineinterface_->computeF(x,F,fillFlag);
    if (err==false)
    {
      cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::computeF:\n"
           << "**ERR**: call to fine-userinterface returned false on level " << level_ << "\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
  }
  else
  {
    RefCountPtr<Epetra_Vector> Ffine = rcp(new Epetra_Vector(fineinterface_->getGraph()->RowMap(),false));
    RefCountPtr<Epetra_Vector> xfine = rcp(prolong_this_to_fine(x));
    err = fineinterface_->computeF(*xfine,*Ffine,fillFlag);
    if (err==false)
    {
      cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::computeF:\n"
           << "**ERR**: call to fine-userinterface returned false on level " << level_ << "\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
    RefCountPtr<Epetra_Vector> Fcoarse = rcp(restrict_fine_to_this(*Ffine));
    F.Scale(1.0,*Fcoarse);
  }
  if (isFAS())
    F.Update(-1.0,*fxbar_,1.0,*fbar_,1.0);
  return err;
}

int
powerMethod (double & lambda, 
             Epetra_CrsMatrix& A, 
             const int niters, 
             const double tolerance,
             const bool verbose)
{
  // In the power iteration, z = A*q.  Thus, q must be in the domain
  // of A, and z must be in the range of A.  The residual vector is of
  // course in the range of A.
  Epetra_Vector q (A.OperatorDomainMap ());
  Epetra_Vector z (A.OperatorRangeMap ());
  Epetra_Vector resid (A.OperatorRangeMap ());

  Epetra_Flops* counter = A.GetFlopCounter();
  if (counter != 0) {
    q.SetFlopCounter(A);
    z.SetFlopCounter(A);
    resid.SetFlopCounter(A);
  }

  // Initialize the starting vector z with random data.
  z.Random();

  double normz, residual;
  int ierr = 1;
  for (int iter = 0; iter < niters; ++iter)
    {
      z.Norm2 (&normz);        // normz  := ||z||_2
      q.Scale (1.0/normz, z);  // q      := z / normz
      A.Multiply(false, q, z); // z      := A * q
      q.Dot(z, &lambda);       // lambda := dot (q, z)

      // Compute the residual vector and display status output every
      // 100 iterations, or if we have reached the maximum number of
      // iterations.
      if (iter % 100 == 0 || iter + 1 == niters)
        {
          resid.Update (1.0, z, -lambda, q, 0.0); // resid := A*q - lambda*q
          resid.Norm2 (&residual);                // residual := ||resid||_2
          if (verbose) 
            cout << "Iter = " << iter << "  Lambda = " << lambda 
                 << "  Residual of A*q - lambda*q = " << residual << endl;
        } 
      if (residual < tolerance) { // We've converged!
        ierr = 0;
        break;
      }
    }
  return ierr;
}

/*----------------------------------------------------------------------*
 |  Constructor (public)                                     m.gee 12/04|
 |  IMPORTANT:                                                          |
 |  No matter on which level we are here, the vector xfine is ALWAYS    |
 |  a fine grid vector here!                                            |
 *----------------------------------------------------------------------*/
ML_NOX::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel(int level, int nlevel, int plevel, 
                                 ML* ml, ML_Aggregate* ag, Epetra_CrsMatrix** P,
                                 ML_NOX::Ml_Nox_Fineinterface& interface, const Epetra_Comm& comm,
                                 const Epetra_Vector& xfine, double fd_alpha, double fd_beta,
                                 bool fd_centered, bool isDiagonalOnly, int bsize) 
: fineinterface_(interface),
comm_(comm)                                                          
{
   level_          = level;
   nlevel_         = nlevel;
   ml_printlevel_  = plevel;
   ml_             = ml;
   ag_             = ag;
   fd_alpha_       = fd_alpha;
   fd_beta_        = fd_beta;
   fd_centered_    = fd_centered;
   isDiagonalOnly_ = isDiagonalOnly;
   A_              = 0;
   coarseinterface_= 0;
   bsize_          = bsize;
   

   // we need the graph of the operator on this level. On the fine grid we can just ask the
   // fineinterface for it, on the coarser levels it has to be extracted from the ML-hierachy
   if (level_==0)
   {
      // the Epetra_CrsGraph-copy-constructor shares data with the original one.
      // We want a really deep copy here so we cannot use it
      // graph_ will be given to the FiniteDifferencing class and will be destroyed by it
      graph_ = ML_NOX::deepcopy_graph(interface.getGraph());
   }
   else
   {
      // Note that ML has no understanding of global indices, so it makes up new GIDs
      // (This also holds for the Prolongators P)
      Epetra_CrsMatrix* tmpMat  = 0;
      int               maxnnz  = 0;
      double            cputime = 0.0;
      ML_Operator2EpetraCrsMatrix(&(ml_->Amat[level_]), tmpMat, maxnnz, false, cputime);
      // copy the graph
      double t0 = GetClock();
      graph_ = ML_NOX::deepcopy_graph(&(tmpMat->Graph()));
      // delete the copy of the Epetra_CrsMatrix
      if (tmpMat) delete tmpMat; tmpMat = 0;
      double t1 = GetClock();
      if (ml_printlevel_ > 0 && 0 == comm_.MyPID())
         cout << "matrixfreeML (level " << level_ << "): extraction/copy of Graph in " << cputime+t1-t0 << " sec\n"
              << "                        max-nonzeros in Graph: " << maxnnz << "\n";
   }
   
   // create this levels coarse interface
   coarseinterface_ = new ML_NOX::Nox_CoarseProblem_Interface(fineinterface_,level_,ml_printlevel_,
                                                              P,&(graph_->RowMap()),nlevel_);

   // restrict the xfine-vector to this level 
   Epetra_Vector* xthis = coarseinterface_->restrict_fine_to_this(xfine);
   if (!xthis)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
          << "**ERR**: ML_Epetra::Nox_CoarseProblem_Interface::restrict_fine_to_this returned NULL on level " 
          << level_ << "\n"
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }

   Epetra_Vector* xc = new Epetra_Vector(graph_->RowMap(),false);
   // FIXME: after intesive testing, this test might be obsolet
#if 0
   bool samemap = xc->Map().PointSameAs(xthis->Map());
   if (samemap)
   {
#endif
      xc->Update(1.0,*xthis,0.0);
#if 0
   }
   else
   {
      cout << "**WRN** Maps are not equal in\n"
           << "**WRN** file/line: " << __FILE__ << "/" << __LINE__ << "\n";
      // import the xthis vector in the Map that ML produced for graph_
      Epetra_Import* importer = new Epetra_Import(graph_->RowMap(),xthis->Map());  
      int ierr = xc->Import(*xthis,*importer,Insert);
      if (ierr)
      {
         cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
              << "**ERR**: export from xthis to xc returned err=" << ierr <<"\n"
              << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
      }
      if (importer) delete importer; importer = 0;
   }
#endif   
   if (xthis) delete xthis; xthis = 0;

   // create the coloring of the graph
   if (ml_printlevel_>0 && comm_.MyPID()==0) 
//.........这里部分代码省略.........

/*----------------------------------------------------------------------*
 |  recreate this level (public)                             m.gee 01/05|
 |  this function assumes, that the graph of the fine level problem has |
 |  not changed since call to the constructor and therefore             |
 | the graph and it's coloring do not have to be recomputed             |  
 |  IMPORTANT:                                                          |
 |  No matter on which level we are here, the vector xfine is ALWAYS    |
 |  a fine grid vector here!                                            |
 *----------------------------------------------------------------------*/
bool ML_NOX::ML_Nox_MatrixfreeLevel::recreateLevel(int level, int nlevel, int plevel, 
                                                   ML* ml, ML_Aggregate* ag, Epetra_CrsMatrix** P,
                                                   ML_NOX::Ml_Nox_Fineinterface& interface, 
                                                   const Epetra_Comm& comm, const Epetra_Vector& xfine) 
{
   // make some tests
   if (level != level_)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::recreateLevel:\n"
          << "**ERR**: level_ " << level_ << " not equal level " << level << "\n"
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }
   if (nlevel != nlevel_)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::recreateLevel:\n"
          << "**ERR**: nlevel_ " << nlevel_ << " not equal nlevel " << nlevel << "\n" 
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }
   // printlevel might have changed
   ml_printlevel_ = plevel;
   ml_            = ml;
   ag_            = ag;
   destroyP();  // safer to use the new Ps
   setP(NULL);

   // we need the graph of the operator on this level. On the fine grid we can just ask the
   // fineinterface for it, on the coarser levels it has to be extracted from the ML-hierachy
   bool same;
   if (level_==0)
   {
      const Epetra_CrsGraph* graph = interface.getGraph();
      // check whether the old graph matches the new one
      same = compare_graphs(graph,graph_);
      destroyFD(); // we are here to recompute the FD-operator (this destroys graph_)
      graph_ = ML_NOX::deepcopy_graph(graph);
   }
   else
   {
      // Note that ML has no understanding of global indices, so it makes up new GIDs
      // (This also holds for the Prolongators P)
      Epetra_CrsMatrix* tmpMat  = 0;
      int               maxnnz  = 0;
      double            cputime = 0.0;
      ML_Operator2EpetraCrsMatrix(&(ml_->Amat[level_]), tmpMat, maxnnz, false, cputime);
      // get a view from the graph
      const Epetra_CrsGraph& graph = tmpMat->Graph();
      // compare the graph to the existing one
      same = compare_graphs(&graph,graph_);
      destroyFD(); // we are here to recompute the FD-operator (this destroys graph_)
      double t0 = GetClock();
      graph_ = ML_NOX::deepcopy_graph(&graph);
      // delete the copy of the Epetra_CrsMatrix
      if (tmpMat) delete tmpMat; tmpMat = 0;
      double t1 = GetClock();
      if (ml_printlevel_ > 0 && 0 == comm_.MyPID())
         cout << "matrixfreeML (level " << level_ << "): extraction/copy of Graph in " << cputime+t1-t0 << " sec\n"
              << "                        max-nonzeros in Graph: " << maxnnz << "\n";
   }
   
   // recreate this levels coarse interface
   if (same)
      coarseinterface_->recreate(ml_printlevel_,P,&(graph_->RowMap()));
   else
   {
      delete coarseinterface_;
      coarseinterface_ = new ML_NOX::Nox_CoarseProblem_Interface(fineinterface_,level_,ml_printlevel_,
                                                                 P,&(graph_->RowMap()),nlevel_);
   }
   
   // restrict the xfine-vector to this level 
   Epetra_Vector* xthis = coarseinterface_->restrict_fine_to_this(xfine);
   if (!xthis)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
          << "**ERR**: ML_Epetra::Nox_CoarseProblem_Interface::restrict_fine_to_this returned NULL on level " 
          << level_ << "\n"
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }

   Epetra_Vector* xc = new Epetra_Vector(graph_->RowMap(),false);
   // FIXME: after intesive testing, this test might be obsolet
#if 0
   bool samemap = xc->Map().PointSameAs(xthis->Map());
   if (samemap)
   {
#endif
      xc->Update(1.0,*xthis,0.0);
#if 0
   }
   else
   {
//.........这里部分代码省略.........

void BiCGSTAB(Epetra_CrsMatrix &A, 
	      Epetra_Vector &x, 
	      Epetra_Vector &b, 
	      Ifpack_CrsRiluk *M, 
	      int Maxiter, 
	      double Tolerance, 
	      double *residual, bool verbose) {

  // Allocate vectors needed for iterations
  Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
  Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
  Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
  Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
  Epetra_Vector r(x.Map()); r.SetFlopCounter(x);
  Epetra_Vector rtilde(x.Map()); rtilde.Random(); rtilde.SetFlopCounter(x);
  Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
  

  A.Multiply(false, x, r); // r = A*x

  r.Update(1.0, b, -1.0); // r = b - A*x

  double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
  double alpha = 1.0, omega = 1.0, sigma;
  double omega_num, omega_den;
  r.Norm2(&r_norm);
  b.Norm2(&b_norm);
  scaled_r_norm = r_norm/b_norm;
  r.Dot(rtilde,&rhon);

  if (verbose) cout << "Initial residual = " << r_norm
		    << " Scaled residual = " << scaled_r_norm << endl;


  for (int i=0; i<Maxiter; i++) { // Main iteration loop   

    double beta = (rhon/rhonm1) * (alpha/omega);
    rhonm1 = rhon;

    /* p    = r + beta*(p - omega*v)       */
    /* phat = M^-1 p                       */
    /* v    = A phat                       */

    double dtemp = - beta*omega;

    p.Update(1.0, r, dtemp, v, beta);
    if (M==0) 
      phat.Scale(1.0, p);
    else
      M->Solve(false, p, phat);
    A.Multiply(false, phat, v);

    
    rtilde.Dot(v,&sigma);
    alpha = rhon/sigma;    

    /* s = r - alpha*v                     */
    /* shat = M^-1 s                       */
    /* r = A shat (r is a tmp here for t ) */

    s.Update(-alpha, v, 1.0, r, 0.0);
    if (M==0) 
      shat.Scale(1.0, s);
    else
      M->Solve(false, s, shat);
    A.Multiply(false, shat, r);

    r.Dot(s, &omega_num);
    r.Dot(r, &omega_den);
    omega = omega_num / omega_den;

    /* x = x + alpha*phat + omega*shat */
    /* r = s - omega*r */

    x.Update(alpha, phat, omega, shat, 1.0);
    r.Update(1.0, s, -omega); 
    
    r.Norm2(&r_norm);
    scaled_r_norm = r_norm/b_norm;
    r.Dot(rtilde,&rhon);

    if (verbose) cout << "Iter "<< i << " residual = " << r_norm
		      << " Scaled residual = " << scaled_r_norm << endl;

    if (scaled_r_norm < Tolerance) break;
  }
  return;
}

  int TestOneMatrix( std::string HBname, std::string MMname, std::string TRIname, Epetra_Comm &Comm, bool verbose ) { 

  if ( Comm.MyPID() != 0 ) verbose = false ; 

  Epetra_Map * readMap = 0;

  Epetra_CrsMatrix * HbA = 0; 
  Epetra_Vector * Hbx = 0; 
  Epetra_Vector * Hbb = 0; 
  Epetra_Vector * Hbxexact = 0;
   
  Epetra_CrsMatrix * TriplesA = 0; 
  Epetra_Vector * Triplesx = 0; 
  Epetra_Vector * Triplesb = 0;
  Epetra_Vector * Triplesxexact = 0;
   
  Epetra_CrsMatrix * MatrixMarketA = 0; 
  Epetra_Vector * MatrixMarketx = 0; 
  Epetra_Vector * MatrixMarketb = 0;
  Epetra_Vector * MatrixMarketxexact = 0;
   
  int TRI_Size = TRIname.size() ; 
  std::string LastFiveBytes = TRIname.substr( EPETRA_MAX(0,TRI_Size-5), TRI_Size );

  if ( LastFiveBytes == ".TimD" ) { 
    // Call routine to read in a file with a Tim Davis header and zero-based indexing
    EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra64( &TRIname[0], false, Comm, 
						      readMap, TriplesA, Triplesx, 
						      Triplesb, Triplesxexact, false, true, true ) );
    delete readMap;
  } else {
    if ( LastFiveBytes == ".triU" ) { 
    // Call routine to read in unsymmetric Triplet matrix
      EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra64( &TRIname[0], false, Comm, 
							readMap, TriplesA, Triplesx, 
							Triplesb, Triplesxexact, false, false ) );
      delete readMap;
    } else {
      if ( LastFiveBytes == ".triS" ) { 
	// Call routine to read in symmetric Triplet matrix
	EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra64( &TRIname[0], true, Comm, 
							  readMap, TriplesA, Triplesx, 
							  Triplesb, Triplesxexact, false, false ) );
        delete readMap;
      } else {
	assert( false ) ; 
      }
    }
  }

  EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra64( &MMname[0], Comm, readMap, 
							 MatrixMarketA, MatrixMarketx, 
							 MatrixMarketb, MatrixMarketxexact) );
  delete readMap;

  // Call routine to read in HB problem
  Trilinos_Util_ReadHb2Epetra64( &HBname[0], Comm, readMap, HbA, Hbx, 
			       Hbb, Hbxexact) ;


#if 0
  std::cout << " HbA " ; 
  HbA->Print( std::cout ) ; 
  std::cout << std::endl ; 

  std::cout << " MatrixMarketA " ; 
  MatrixMarketA->Print( std::cout ) ; 
  std::cout << std::endl ; 

  std::cout << " TriplesA " ; 
  TriplesA->Print( std::cout ) ; 
  std::cout << std::endl ; 
#endif


  int TripleErr = 0 ; 
  int MMerr = 0 ; 
  for ( int i = 0 ; i < 10 ; i++ ) 
    {
      double resid_Hb_Triples;
      double resid_Hb_Matrix_Market;
      double norm_A ;
      Hbx->Random();
      //
      //  Set the output vectors to different values:
      //
      Triplesb->PutScalar(1.1);
      Hbb->PutScalar(1.2);
      MatrixMarketb->PutScalar(1.3);

      HbA->Multiply( false, *Hbx, *Hbb );
      norm_A = HbA->NormOne( ) ; 

      TriplesA->Multiply( false, *Hbx, *Triplesb );
      Triplesb->Update( 1.0, *Hbb, -1.0 ) ; 


      MatrixMarketA->Multiply( false, *Hbx, *MatrixMarketb );
      MatrixMarketb->Update( 1.0, *Hbb, -1.0 ) ; 

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

// ================================================ ====== ==== ==== == =
// the tentative null space is in input because the user
// has to remember to allocate and fill it, and then to delete
// it after calling this method.
int ML_Epetra::MultiLevelPreconditioner::
ComputeAdaptivePreconditioner(int TentativeNullSpaceSize,
                              double* TentativeNullSpace)
{

  if ((TentativeNullSpaceSize == 0) || (TentativeNullSpace == 0))
    ML_CHK_ERR(-1);
   
  // ================================== //
  // get parameters from the input list //
  // ================================== //
  
  // maximum number of relaxation sweeps
  int MaxSweeps = List_.get("adaptive: max sweeps", 10);
  // number of std::vector to be added to the tentative null space
  int NumAdaptiveVectors = List_.get("adaptive: num vectors", 1);

  if (verbose_) {
    std::cout << PrintMsg_ << "*** Adaptive Smoother Aggregation setup ***" << std::endl;
    std::cout << PrintMsg_ << "    Maximum relaxation sweeps     = " << MaxSweeps << std::endl;
    std::cout << PrintMsg_ << "    Additional vectors to compute = " << NumAdaptiveVectors << std::endl;
  }

  // ==================================================== //
  // compute the preconditioner, set null space from user //
  // (who will have to delete std::vector TentativeNullSpace)  //
  // ==================================================== //
  
  double* NewNullSpace = 0;
  double* OldNullSpace = TentativeNullSpace;
  int OldNullSpaceSize = TentativeNullSpaceSize;

  // need some work otherwise matvec() with Epetra_Vbr fails.
  // Also, don't differentiate between range and domain here
  // as ML will not work if range != domain
  const Epetra_VbrMatrix* VbrA = NULL;
  VbrA = dynamic_cast<const Epetra_VbrMatrix*>(RowMatrix_);

  Epetra_Vector* LHS = 0;
  Epetra_Vector* RHS = 0;

  if (VbrA != 0) {
    LHS = new Epetra_Vector(VbrA->DomainMap());
    RHS = new Epetra_Vector(VbrA->DomainMap());
  } else {
    LHS = new Epetra_Vector(RowMatrix_->OperatorDomainMap());
    RHS = new Epetra_Vector(RowMatrix_->OperatorDomainMap());
  }

  // destroy what we may already have
  if (IsComputePreconditionerOK_ == true) {
    DestroyPreconditioner();
  }

  // build the preconditioner for the first time
  List_.set("null space: type", "pre-computed");
  List_.set("null space: dimension", OldNullSpaceSize);
  List_.set("null space: vectors", OldNullSpace);
  ComputePreconditioner();

  // ====================== //
  // add one std::vector at time //
  // ====================== //
  
  for (int istep = 0 ; istep < NumAdaptiveVectors ; ++istep) {

    if (verbose_) {
      std::cout << PrintMsg_ << "\tAdaptation step " << istep << std::endl;
      std::cout << PrintMsg_ << "\t---------------" << std::endl;
    }

    // ==================== //
    // look for "bad" modes //
    // ==================== //

    // note: should an error occur, ML_CHK_ERR will return,
    // and LHS and RHS will *not* be delete'd (--> memory leak).
    // Anyway, this means that something wrong happened in the code
    // and should be fixed by the user.

    LHS->Random();
    double Norm2;

    for (int i = 0 ; i < MaxSweeps ; ++i) {
      // RHS = (I - ML^{-1} A) LHS
      ML_CHK_ERR(RowMatrix_->Multiply(false,*LHS,*RHS));
      // FIXME: can do something slightly better here
      ML_CHK_ERR(ApplyInverse(*RHS,*RHS));
      ML_CHK_ERR(LHS->Update(-1.0,*RHS,1.0));
      LHS->Norm2(&Norm2);
      if (verbose_) {
        std::cout << PrintMsg_ << "\titer " << i << ", ||x||_2 = ";
        std::cout << Norm2 << std::endl;
      }
    }

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

//.........这里部分代码省略.........
    // global entries per process, but will distribute them differently,
    // in round-robin (1-D cyclic) fashion instead of contiguously.

    // We'll use the version of the Map constructor that takes, on each
    // MPI process, a list of the global indices in the Map belonging to
    // that process.  You can use this constructor to construct an
    // overlapping (also called "not 1-to-1") Map, in which one or more
    // entries are owned by multiple processes.  We don't do that here;
    // we make a nonoverlapping (also called "1-to-1") Map.
    const int numGblIndsPerProc = 5;
    global_ordinal_type* gblIndList = new global_ordinal_type [numGblIndsPerProc];

    const int numProcs = comm.NumProc ();
    const int myRank = comm.MyPID ();
    for (int k = 0; k < numGblIndsPerProc; ++k) {
        gblIndList[k] = myRank + k*numProcs;
    }

    Epetra_Map cyclicMap (numGlobalEntries, numGblIndsPerProc,
                          gblIndList, indexBase, comm);
    // The above constructor makes a deep copy of the input index list,
    // so it's safe to deallocate that list after this constructor
    // completes.
    if (gblIndList != NULL) {
        delete [] gblIndList;
        gblIndList = NULL;
    }

    // If there's more than one MPI process in the communicator,
    // then cyclicMap is definitely NOT contiguous.
    if (comm.NumProc () > 1 && cyclicMap.LinearMap ()) {
        throw std::logic_error ("The cyclic Map claims to be contiguous.");
    }

    // contigMap and cyclicMap should always be compatible.  However, if
    // the communicator contains more than 1 process, then contigMap and
    // cyclicMap are NOT the same.
    // if (! contigMap.isCompatible (*cyclicMap)) {
    //   throw std::logic_error ("contigMap should be compatible with cyclicMap, "
    //                           "but it's not.");
    // }
    if (comm.NumProc () > 1 && contigMap.SameAs (cyclicMap)) {
        throw std::logic_error ("contigMap should not be the same as cyclicMap.");
    }

    //////////////////////////////////////////////////////////////////////
    // We have maps now, so we can create vectors.
    //////////////////////////////////////////////////////////////////////

    // Create an Epetra_Vector with the contiguous Map we created above.
    // This version of the constructor will fill the vector with zeros.
    // The Vector constructor takes a Map by value, but that's OK,
    // because Epetra_Map has shallow copy semantics.  It uses reference
    // counting internally to avoid copying data unnecessarily.
    Epetra_Vector x (contigMap);

    // The copy constructor performs a deep copy.
    // x and y have the same Map.
    Epetra_Vector y (x);

    // Create a Vector with the 1-D cyclic Map.  Calling the constructor
    // with false for the second argument leaves the data uninitialized,
    // so that you can fill it later without paying the cost of
    // initially filling it with zeros.
    Epetra_Vector z (cyclicMap, false);

    // Set the entries of z to (pseudo)random numbers.  Please don't
    // consider this a good parallel pseudorandom number generator.
    (void) z.Random ();

    // Set the entries of x to all ones.
    (void) x.PutScalar (1.0);

    // Define some constants for use below.
    const double alpha = 3.14159;
    const double beta = 2.71828;
    const double gamma = -10.0;

    // x = beta*x + alpha*z
    //
    // This is a legal operation!  Even though the Maps of x and z are
    // not the same, their Maps are compatible.  Whether it makes sense
    // or not depends on your application.
    (void) x.Update (alpha, z, beta);

    (void) y.PutScalar (42.0); // Set all entries of y to 42.0
    // y = gamma*y + alpha*x + beta*z
    y.Update (alpha, x, beta, z, gamma);

    // Compute the 2-norm of y.
    //
    // The norm may have a different type than scalar_type.
    // For example, if scalar_type is complex, then the norm is real.
    // The ScalarTraits "traits class" gives us the type of the norm.
    double theNorm = 0.0;
    (void) y.Norm2 (&theNorm);

    // Print the norm of y on Proc 0.
    out << "Norm of y: " << theNorm << endl;
}

//.........这里部分代码省略.........
      EPETRA_CHK_ERR( A_umfpack.SymbolicFactorization(  ) ); 
      EPETRA_CHK_ERR( A_umfpack.NumericFactorization(  ) ); 
      EPETRA_CHK_ERR( A_umfpack.Solve(  ) ); 
#endif
#ifdef HAVE_AMESOS_KLU
    } else if ( SparseSolver == KLU ) {


      using namespace Teuchos;

      Amesos_Time AT; 
      int setupTimePtr = -1, symTimePtr = -1, numTimePtr = -1, refacTimePtr = -1, solveTimePtr = -1;
      AT.CreateTimer(Comm, 2);
      AT.ResetTimer(0);

      Teuchos::ParameterList ParamList ;
      // ParamList.set("OutputLevel",2);
      Amesos_Klu A_klu( Problem ); 
      ParamList.set( "MaxProcs", -3 );
      ParamList.set( "TrustMe", false );
      // ParamList.set( "Refactorize", true );
      EPETRA_CHK_ERR( A_klu.SetParameters( ParamList ) ) ; 
      EPETRA_CHK_ERR( A_klu.SetUseTranspose( transpose ) ); 
      setupTimePtr = AT.AddTime("Setup", setupTimePtr, 0);
      EPETRA_CHK_ERR( A_klu.SymbolicFactorization(  ) ); 
      symTimePtr = AT.AddTime("Symbolic", symTimePtr, 0);
      EPETRA_CHK_ERR( A_klu.NumericFactorization(  ) ); 
      numTimePtr = AT.AddTime("Numeric", numTimePtr, 0);
      EPETRA_CHK_ERR( A_klu.NumericFactorization(  ) ); 
      refacTimePtr = AT.AddTime("Refactor", refacTimePtr, 0);
      // for ( int i=0; i<100000 ; i++ ) 
      EPETRA_CHK_ERR( A_klu.Solve(  ) ); 
      solveTimePtr = AT.AddTime("Solve", solveTimePtr, 0);

      double SetupTime = AT.GetTime(setupTimePtr);
      double SymbolicTime = AT.GetTime(symTimePtr);
      double NumericTime = AT.GetTime(numTimePtr);
      double RefactorTime = AT.GetTime(refacTimePtr);
      double SolveTime = AT.GetTime(solveTimePtr);

      std::cout << __FILE__ << "::"  << __LINE__ << " SetupTime = " << SetupTime << std::endl ; 
      std::cout << __FILE__ << "::"  << __LINE__ << " SymbolicTime = " << SymbolicTime - SetupTime << std::endl ; 
      std::cout << __FILE__ << "::"  << __LINE__ << " NumericTime = " << NumericTime - SymbolicTime<< std::endl ; 
      std::cout << __FILE__ << "::"  << __LINE__ << " RefactorTime = " << RefactorTime - NumericTime << std::endl ; 
      std::cout << __FILE__ << "::"  << __LINE__ << " SolveTime = " << SolveTime - RefactorTime << std::endl ; 

#endif
    } else { 
      SparseDirectTimingVars::log_file << "Solver not implemented yet" << std::endl ;
      std::cerr << "\n\n####################  Requested solver not available on this platform ##################### ATS\n" << std::endl ;
      std::cout << " SparseSolver = " << SparseSolver << std::endl ; 
      std::cerr << " SparseSolver = " << SparseSolver << std::endl ; 
    }
    
    SparseDirectTimingVars::SS_Result.Set_Total_Time( TotalTime.ElapsedTime() ); 
  }  // end for (int i=0; i<special; i++ ) 

  //
  //  Compute the error = norm(xcomp - xexact )
  //
  double error;
  passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);

  passresid->Norm2(&error);
  SparseDirectTimingVars::SS_Result.Set_Error(error) ;

  //  passxexact->Norm2(&error ) ; 
  //  passx->Norm2(&error ) ; 

  //
  //  Compute the residual = norm(Ax - b)
  //
  double residual ; 

  passA->Multiply( transpose, *passx, *passtmp);
  passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0); 
  //  passresid->Update(1.0, *passtmp, -1.0, CopyB, 0.0); 
  passresid->Norm2(&residual);

  SparseDirectTimingVars::SS_Result.Set_Residual(residual) ;
    
  double bnorm; 
  passb->Norm2( &bnorm ) ; 
  SparseDirectTimingVars::SS_Result.Set_Bnorm(bnorm) ;

  double xnorm; 
  passx->Norm2( &xnorm ) ; 
  SparseDirectTimingVars::SS_Result.Set_Xnorm(xnorm) ;

  delete readA;
  delete readx;
  delete readb;
  delete readxexact;
  delete readMap;
  delete map_;
  
  Comm.Barrier();

  return 0;
}