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

//EpetraVector_To_TpetraVectorNonConst: copies Epetra_Vector to non-const Tpetra_Vector
Teuchos::RCP<Tpetra_Vector> Petra::EpetraVector_To_TpetraVectorNonConst(const Epetra_Vector& epetraVector_,
                                                               const Teuchos::RCP<const Teuchos::Comm<int> >& commT_)
{
  //get map from epetraVector_ and convert to Tpetra::Map
  auto mapT = EpetraMap_To_TpetraMap(epetraVector_.Map(), commT_);
  ST *values;
  epetraVector_.ExtractView(&values);
  Teuchos::ArrayView<ST> valuesAV = Teuchos::arrayView(values, mapT->getNodeNumElements());
  Teuchos::RCP<Tpetra_Vector> tpetraVector_ = Teuchos::rcp(new Tpetra_Vector(mapT, valuesAV));
  return tpetraVector_;
}

// =============================================================================
Epetra_NumPyVector::Epetra_NumPyVector(const Epetra_Vector & source):
  Epetra_Vector(source)
{
  map = new Epetra_BlockMap(source.Map());
  npy_intp dims[ ] = { map->NumMyPoints() };
  double *v = NULL;
  Epetra_Vector::ExtractView(&v);
  array = (PyArrayObject *)
    PyArray_SimpleNewFromData(1,dims,PyArray_DOUBLE,(void *)v);
  if (!array)
  {
    cleanup();
    throw PythonException();
  }
}

//=============================================================================
//=============================================================================
int Epetra_MsrMatrix::InvColSums(Epetra_Vector& x) const {
//
// Put inverse of the sum of absolute values of the jth column of A in x[j].
//

  if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
  if (!OperatorDomainMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
  

  Epetra_Vector * xp = 0;
  Epetra_Vector * x_tmp = 0;
  

  // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
  if (RowMatrixImporter()!=0) {
    x_tmp = new Epetra_Vector(RowMatrixColMap()); // Create import vector if needed
    xp = x_tmp;
  }
  int ierr = 0;
  int i, j;

  for (i=0; i < NumMyCols_; i++) (*xp)[i] = 0.0;

  for (i=0; i < NumMyRows_; i++) {
    int NumEntries = GetRow(i);// Copies ith row of matrix into Values_ and Indices_
    for (j=0; j < NumEntries; j++) (*xp)[Indices_[j]] += fabs(Values_[j]);
  }

  if (RowMatrixImporter()!=0){
    x.Export(*x_tmp, *RowMatrixImporter(), Add); // Fill x with Values from import vector
    delete x_tmp;
    xp = &x;
  }
  // Invert values, don't allow them to get too large
  for (i=0; i < NumMyRows_; i++) {
    double scale = (*xp)[i];
    if (scale<Epetra_MinDouble) {
      if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
      else if (ierr!=1) ierr = 2;
      (*xp)[i] = Epetra_MaxDouble;
    }
    else
      (*xp)[i] = 1.0/scale;
  }
  UpdateFlops(NumGlobalNonzeros());
  EPETRA_CHK_ERR(ierr);
  return(0);
}

void
NOX::Epetra::DebugTools::writeVector( std::string baseName, const Epetra_Vector & vec, FORMAT_TYPE outFormat, bool writeMap )
{
  if( ASCII == outFormat )
  {
    vec.Print( std::cout );
  }
  else
  {
    std::string fileName = baseName + "_vec";
    EpetraExt::VectorToMatrixMarketFile(fileName.c_str(), vec);
    if( writeMap )
    {
      fileName = baseName + "_map";
      EpetraExt::BlockMapToMatrixMarketFile(fileName.c_str(), vec.Map());
    }
  }
}

NOX::Epetra::Vector::Vector(const Epetra_Vector& source, NOX::CopyType type,
                Teuchos::RCP<NOX::Epetra::VectorSpace> vs)
{
  if (Teuchos::is_null(vs))
    vectorSpace = Teuchos::rcp(new NOX::Epetra::VectorSpaceL2);
  else
    vectorSpace = vs;

  switch (type) {

  case DeepCopy:        // default behavior

    epetraVec = Teuchos::rcp(new Epetra_Vector(source));
    break;

  case ShapeCopy:

    epetraVec = Teuchos::rcp(new Epetra_Vector(source.Map()));
    break;

  }
}

void
Albany::SolutionMaxValueResponseFunction::
evaluateGradient(const double current_time,
		 const Epetra_Vector* xdot,
		 const Epetra_Vector* xdotdot,
		 const Epetra_Vector& x,
		 const Teuchos::Array<ParamVec>& p,
		 ParamVec* deriv_p,
		 Epetra_Vector* g,
		 Epetra_MultiVector* dg_dx,
		 Epetra_MultiVector* dg_dxdot,
		 Epetra_MultiVector* dg_dxdotdot,
		 Epetra_MultiVector* dg_dp)
{
  int global_index;
  double mxv;
  computeMaxValue(x, mxv, global_index);

  // Evaluate response g
  if (g != NULL)
    (*g)[0] = mxv;

  // Evaluate dg/dx
  if (dg_dx != NULL) {
    dg_dx->PutScalar(0.0);
    int lid = x.Map().LID(global_index);
    if(lid >= 0) (*dg_dx)[0][lid] = 1.0;
  }

  // Evaluate dg/dxdot
  if (dg_dxdot != NULL)
    dg_dxdot->PutScalar(0.0);
  if (dg_dxdotdot != NULL)
    dg_dxdotdot->PutScalar(0.0);

  // Evaluate dg/dp
  if (dg_dp != NULL)
    dg_dp->PutScalar(0.0);
}

void  EpetraGhostView::import(const Epetra_Import& importer,
  const Epetra_Vector& srcObject)
{
  /* If my vector does not yet exist, create it using the target map of the
   * importer */
  if (ghostView_.get()==0)
  {
    ghostView_ = rcp(new Epetra_Vector(importer.TargetMap()));
  }

  /* do the import */
  int ierr = ghostView_->Import(srcObject, importer, Insert);

  if (ierr < 0)
  {
    Out::os() << "target map=" << endl;
    importer.TargetMap().Print(Out::os());
    Out::os() << "source map=" << endl;
    srcObject.Map().Print(Out::os());
  }
  TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error, "ierr=" << ierr << " in EpetraGhostView::import()");
}

//=============================================================================
int Epetra_MsrMatrix::RightScale(const Epetra_Vector& x) {
//
// This function scales the jth row of A by x[j].
//

  // For this method, we have no choice but to work with the UGLY MSR data structures.

  if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
  if (!OperatorDomainMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A

  int * bindx = Amat_->bindx;
  double * val = Amat_->val;
  Epetra_Vector * xp = 0;
  Epetra_Vector * x_tmp = 0;

  // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
  if (RowMatrixImporter()!=0) {
    x_tmp = new Epetra_Vector(RowMatrixColMap()); // Create import vector if needed
    x_tmp->Import(x,*RowMatrixImporter(), Insert); // x_tmp will have all the values we need
    xp = x_tmp;
  }

  int i, j;

  for (i=0; i < NumMyRows_; i++) {
    int NumEntries = bindx[i+1] - bindx[i];
    double scale = (*xp)[i];
    val[i] *= scale;
    double * Values = val + bindx[i];
    int * Indices = bindx + bindx[i];
    for (j=0; j < NumEntries; j++)  Values[j] *=  (*xp)[Indices[j]];
  }
  if (x_tmp!=0) delete x_tmp;
  NormOne_ = -1.0; // Reset Norm so it will be recomputed.
  NormInf_ = -1.0; // Reset Norm so it will be recomputed.
  UpdateFlops(NumGlobalNonzeros());
  return(0);
}

void gather_in_block(const std::string & blockId, const panzer::UniqueGlobalIndexer<int,GlobalOrdinal> & dofMngr,
                     const Epetra_Vector & x,const std::vector<std::size_t> & localCellIds,
                     std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > & fc)
{
   const std::vector<int> & fieldNums = dofMngr.getBlockFieldNumbers(blockId);

   for(std::size_t fieldIndex=0;fieldIndex<fieldNums.size();fieldIndex++) {
      int fieldNum = fieldNums[fieldIndex];
      std::string fieldStr = dofMngr.getFieldString(fieldNum);

      // grab the field
      const std::vector<int> & elmtOffset = dofMngr.getGIDFieldOffsets(blockId,fieldNum);
      fc[fieldStr] = Kokkos::DynRankView<double,PHX::Device>("fc",localCellIds.size(),elmtOffset.size());

      // gather operation for each cell in workset
      for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
         std::vector<GlobalOrdinal> GIDs;
         std::vector<int> LIDs;
         std::size_t cellLocalId = localCellIds[worksetCellIndex];
      
         dofMngr.getElementGIDs(cellLocalId,GIDs);
      
         // caculate the local IDs for this element
         LIDs.resize(GIDs.size());
         for(std::size_t i=0;i<GIDs.size();i++)
            LIDs[i] = x.Map().LID(GIDs[i]);
   
         // loop over basis functions and fill the fields
         for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
            int offset = elmtOffset[basis];
            int lid = LIDs[offset];
            fc[fieldStr](worksetCellIndex,basis) = x[lid];
         }
      }
   }
}

/*----------------------------------------------------------------------*
 |  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 Amesos_Lapack::GEEV(Epetra_Vector& Er, Epetra_Vector& Ei)
{
  if (IsSymbolicFactorizationOK_ == false)
    AMESOS_CHK_ERR(SymbolicFactorization());

  if (MyPID_ == 0)
    AMESOS_CHK_ERR(DenseMatrix_.Shape(static_cast<int>(NumGlobalRows64()),static_cast<int>(NumGlobalRows64())));

  AMESOS_CHK_ERR(DistributedToSerial());
  AMESOS_CHK_ERR(SerialToDense());

  Teuchos::RCP<Epetra_Vector> LocalEr;
  Teuchos::RCP<Epetra_Vector> LocalEi;

  if (NumProcs_ == 1)
  {
    LocalEr = Teuchos::rcp(&Er, false);
    LocalEi = Teuchos::rcp(&Ei, false);
  }
  else
  {
    LocalEr = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
    LocalEi = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
  }

  if (MyPID_ == 0) 
  {
    int n = static_cast<int>(NumGlobalRows64());
    char jobvl = 'N'; /* V/N to calculate/not calculate left eigenvectors
                         of matrix H.*/
    char jobvr = 'N'; /* As above, but for right eigenvectors. */
    int info = 0;
    int ldvl = n;
    int ldvr = n;

    Er.PutScalar(0.0);
    Ei.PutScalar(0.0);

    Epetra_LAPACK LAPACK;

    std::vector<double> work(1);
    int lwork = -1;

    LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
                LocalEr->Values(), LocalEi->Values(), NULL,
                ldvl, NULL, 
                ldvr, &work[0], 
                lwork, &info);

    lwork = (int)work[0];
    work.resize(lwork);
    LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
                LocalEr->Values(), LocalEi->Values(), NULL,
                ldvl, NULL, 
                ldvr, &work[0], 
                lwork, &info);

    if (info)
      AMESOS_CHK_ERR(info);
  }

  if (NumProcs_ != 1)
  {
    // I am not really sure that exporting the results make sense... 
    // It is just to be coherent with the other parts of the code.
    Er.Import(*LocalEr, Epetra_Import(Er.Map(), SerialMap()), Insert);
    Ei.Import(*LocalEi, Epetra_Import(Ei.Map(), SerialMap()), Insert);
  }

  return(0);
}

AztecOO_StatusType
AztecOO_StatusTestResNorm::CheckStatus(int CurrentIter, 
                                       Epetra_MultiVector * CurrentResVector, 
                                       double CurrentResNormEst,  
                                       bool SolutionUpdated)
{
  (void)CurrentIter;
  Epetra_Vector * crv = dynamic_cast<Epetra_Vector *>(CurrentResVector);

  // This section computes the norm of the residual vector
  if (restype_==Implicit && resnormtype_==TwoNorm && CurrentResNormEst!=-1.0) 
    resvalue_ = CurrentResNormEst;
  else if (crv==0) { // Cannot proceed because there is no norm est or res vector
    status_ = Failed;
    return(status_);
  }
  else if (restype_==Explicit && SolutionUpdated) {
    curresvecexplicit_ = true;
    if (localresvector_==0) localresvector_ = new Epetra_Vector(crv->Map());
    // Compute explicit residual
    operator_.Apply(lhs_, *localresvector_);
    localresvector_->Update(1.0, rhs_, -1.0); // localresvector_ = rhs_ - operator_* lhs_
    if (resweights_!=0) { // Check if we should scale the vector
      // localresvector_ = resweights_ * localresvector_
      localresvector_->Multiply(1.0, *resweights_, *localresvector_, 0.0);
    }
    resvalue_ = ComputeNorm(*localresvector_, resnormtype_);
  }
  else {
    curresvecexplicit_ = false;
    if (resweights_!=0) { // Check if we should scale the vector
      if (localresvector_==0) localresvector_ = new Epetra_Vector(crv->Map());
      // localresvector_ = resweights_ * localresvector_
      localresvector_->Multiply(1.0, *resweights_, *crv, 0.0);
      resvalue_ = ComputeNorm(*localresvector_, resnormtype_);
    }
    else
      resvalue_ = ComputeNorm(*crv, resnormtype_);
  }


  // Compute scaling term (done once)
  if (firstcallCheckStatus_) {
    if (scaletype_==NormOfRHS) {
      if (scaleweights_!=0) { // Check if we should scale the vector
	if (localresvector_==0) localresvector_ = new Epetra_Vector(rhs_.Map());
	// localresvector = scaleweights_ * rhs_
	localresvector_->Multiply(1.0, *scaleweights_, rhs_, 0.0);
	scalevalue_ = ComputeNorm(*localresvector_, resnormtype_);
      }
      else {
	scalevalue_ = ComputeNorm(rhs_, scalenormtype_);
      }
    }
    else if (scaletype_==NormOfInitRes) {
      if (restype_==Implicit && scalenormtype_==TwoNorm && CurrentResNormEst!=-1.0) 
	scalevalue_ = CurrentResNormEst;
      else {
	if (scaleweights_!=0) { // Check if we should scale the vector
	  if (localresvector_==0) localresvector_ = new Epetra_Vector(crv->Map());
	  // weightedrhs = scaleweights_ * initial residual
	  localresvector_->Multiply(1.0, *scaleweights_, *crv, 0.0);
	  scalevalue_ = ComputeNorm(*localresvector_, resnormtype_);
	}
	else {
	  scalevalue_ = ComputeNorm(rhs_, scalenormtype_);
	}
      }
    }
    if (scalevalue_==0.0) {
      status_ = Failed;
      return(status_);
    } 
  }

  testvalue_ = resvalue_/scalevalue_;
  if (testvalue_>tolerance_) 
    if (convergedOnce_) {
      if (numExtraIterations_<maxNumExtraIterations_) {
	numExtraIterations_++;
	status_ = Unconverged;
      }
      else {
	status_ = PartialFailed;
      }
    }
    else {
      status_ = Unconverged;
    }
  else if (testvalue_<=tolerance_) {
    convergedOnce_ = true;	  
    status_ = Converged;
  }
  else
    status_ = NaN;

  firstcallCheckStatus_ = false;
  return status_;
}

bool FiniteDifferenceColoringWithUpdate::differenceProbe(const Epetra_Vector& x, Epetra_CrsMatrix& jac,const Epetra_MapColoring& colors){

  // Allocate space for perturbation, get column version of x for scaling
  Epetra_Vector xp(x);
  Epetra_Vector *xcol;
  int N=jac.NumMyRows();

  if(jac.ColMap().SameAs(x.Map()))
     xcol=const_cast<Epetra_Vector*>(&x);
  else{
    xcol=new Epetra_Vector(jac.ColMap(),true);//zeros out by default
    xcol->Import(x,*jac.Importer(),InsertAdd);
  }

  // Counters for probing diagnostics
  double tmp,probing_error_lower_bound=0.0,jc_norm=0.0;

  // Grab coloring info (being very careful to ignore color 0)
  int Ncolors=colors.MaxNumColors()+1;
  int num_c0_global,num_c0_local=colors.NumElementsWithColor(0);
  colors.Comm().MaxAll(&num_c0_local,&num_c0_global,1);
  if(num_c0_global>0) Ncolors--;

  if(Ncolors==0) return false;

  // Pointers for Matrix Info
  int entries, *indices;
  double *values;

  // NTS: Fix me
  if ( diffType == Centered ) exit(1);

  double scaleFactor = 1.0;
  if ( diffType == Backward )
    scaleFactor = -1.0;

  // Compute RHS at initial solution
  computeF(x,fo,NOX::Epetra::Interface::Required::FD_Res);

  /* Probing, vector by vector since computeF does not have a MultiVector interface */
  // Assume that anything with Color 0 gets ignored.
  for(int j=1;j<Ncolors;j++){
    xp=x;
    for(int i=0;i<N;i++){
      if(colors[i]==j)
    xp[i] += scaleFactor*(alpha*abs(x[i])+beta);
    }

    computeF(xp, fp, NOX::Epetra::Interface::Required::FD_Res);

    // Do the subtraction to estimate the Jacobian (w/o including step length)
    Jc.Update(1.0, fp, -1.0, fo, 0.0);

    // Relative error in probing
     if(use_probing_diags){
       Jc.Norm2(&tmp);
       jc_norm+=tmp*tmp;
     }

    for(int i=0;i<N;i++){
      // Skip for uncolored row/columns, else update entries
      if(colors[i]==0) continue;

      jac.ExtractMyRowView(i,entries,values,indices);
      for(int k=0;k<jac.NumMyEntries(i);k++){
    if(colors[indices[k]]==j){
      values[k]=Jc[i] / (scaleFactor*(alpha*abs((*xcol)[indices[k]])+beta));
      // If probing diagnostics are on, zero out the entries as they are used
      if(use_probing_diags) Jc[i]=0.0;
      break;// Only one value per row...
    }
      }
    }
    if(use_probing_diags){
      Jc.Norm2(&tmp);
      probing_error_lower_bound+=tmp*tmp;
    }
  }

  // If diagnostics are requested, output Frobenius norm lower bound
  if(use_probing_diags && !x.Comm().MyPID()) printf("Probing Error Lower Bound (Frobenius) abs = %6.4e rel = %6.4e\n",sqrt(probing_error_lower_bound),sqrt(probing_error_lower_bound)/sqrt(jc_norm));

  // Cleanup
  if(!jac.ColMap().SameAs(x.Map()))
    delete xcol;

  return true;
}