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

Epetra_CrsMatrix *
NOX::Epetra::DebugTools::compute_matrix_using_operator( const Epetra_Operator * op)
{

  const Epetra_Map & rowMap  =     op->OperatorDomainMap();
  Epetra_CrsMatrix * p_mat   = new Epetra_CrsMatrix(Copy, rowMap, 0);
  Epetra_Vector *    tempVec = new Epetra_Vector(rowMap);
  Epetra_Vector *    tempRes = new Epetra_Vector(rowMap);

  // Show progress on what could be a long operation
  if( 0 == op->Comm().MyPID() )
  {
    std::cout << "****************  CREATING MATRIX FROM OPERATOR ************************ "
          << std::endl;
    std::cout << NOX::Utils::fill(72) << std::endl;
  }
  int totalPerturbations = tempVec->GlobalLength();
  int outFreq = totalPerturbations / 71;
  if( 1 > outFreq )
    outFreq = 1;

  for( int col = 0; col < tempVec->GlobalLength(); ++col )
  {
    tempVec->PutScalar(0.0);
    if( rowMap.MyGID(col) )
      (*tempVec)[col] = 1.0;

    op->Apply(*tempVec, *tempRes);

    // Fill in only the nonzero entries from the apply
    for( int row = 0; row < p_mat->NumMyRows(); ++row)
    {
      if( fabs( (*tempRes)[row] ) > 1.e-12 )
      {
    int ierr = p_mat->InsertGlobalValues( rowMap.GID(row), 1, &(*tempRes)[row], &col );
    if( ierr < 0 )
    {
          std::string msg = //"ERROR (" + ierr + ") : "
           "NOX::Epetra::DebugTools::compute_matrix_using_operator crsMatrix.ExtractGlobalRowView(...) failed for row : ";// + row;
          throw msg;
    }
      }
    }
    if( (0 == op->Comm().MyPID()) && (0 == (col % outFreq)) )
      std::cout << "-" << std::flush;
  }
  if( 0 == op->Comm().MyPID() )
    std::cout << "*" << std::endl;

  p_mat->FillComplete();

  delete tempRes;
  delete tempVec;

  return p_mat;
}

void EpetraLinearOp::computeAbsRowSum(Epetra_Vector & x) const
{
  TEUCHOS_ASSERT(!is_null(rowMatrix_));

  RCP<Epetra_CrsMatrix> crsMatrix 
    = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(rowMatrix_);

  TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
    Exceptions::OpNotSupported,
    "EpetraLinearOp::computeAbsRowSum(...): wrapped matrix must be of type "
    "Epetra_CrsMatrix for this method. Other operator types are not supported."
    );

  //
  // Put inverse of the sum of absolute values of the ith row of A in x[i].
  // (this is a modified copy of Epetra_CrsMatrix::InvRowSums)
  //

  if (crsMatrix->Filled()) {
    TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
      std::invalid_argument,
      "EpetraLinearOp::computeAbsRowSum(...): Epetra_CrsMatrix must be filled"
    );
  } 
  int i, j;
  x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
  double * xp = (double*)x.Values();
  if (crsMatrix->Graph().RangeMap().SameAs(x.Map()) && crsMatrix->Exporter() != 0) {
    Epetra_Vector x_tmp(crsMatrix->RowMap());
    x_tmp.PutScalar(0.0);
    double * x_tmp_p = (double*)x_tmp.Values();
    for (i=0; i < crsMatrix->NumMyRows(); i++) {
      int      NumEntries = 0;
      double * RowValues  = 0;
      crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
      for (j=0; j < NumEntries; j++)  x_tmp_p[i] += std::abs(RowValues[j]);
    }
    TEUCHOS_TEST_FOR_EXCEPT(0!=x.Export(x_tmp, *crsMatrix->Exporter(), Add)); //Export partial row sums to x.
  }
  else if (crsMatrix->Graph().RowMap().SameAs(x.Map())) {
    for (i=0; i < crsMatrix->NumMyRows(); i++) {
      int      NumEntries = 0;
      double * RowValues  = 0;
      crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
      double scale = 0.0;
      for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
      xp[i] = scale;
    }
  }
  else { // x.Map different than both crsMatrix->Graph().RowMap() and crsMatrix->Graph().RangeMap()
    TEUCHOS_TEST_FOR_EXCEPT(true); // The map of x must be the RowMap or RangeMap of A.
  }
}

//.........这里部分代码省略.........
  double abstol = 1.e-4;
  double reltol = 1.e-4 ;

  // Test NOX::Epetra::BroydenOperator
  int numGlobalElems = 3 * NumProc;
  Epetra_Map      broydenRowMap   ( numGlobalElems, 0, Comm );
  Epetra_Vector   broydenWorkVec  ( broydenRowMap );
  Epetra_CrsGraph broydenWorkGraph( Copy, broydenRowMap, 0 );
  std::vector<int> globalIndices(3);
  for( int lcol = 0; lcol < 3; ++lcol )
    globalIndices[lcol] = 3 * MyPID + lcol;

  std::vector<int> myGlobalIndices(2);

  // Row 1 structure
  myGlobalIndices[0] = globalIndices[0];
  myGlobalIndices[1] = globalIndices[2];
  broydenWorkGraph.InsertGlobalIndices( globalIndices[0], 2, &myGlobalIndices[0] );
  // Row 2 structure
  myGlobalIndices[0] = globalIndices[0];
  myGlobalIndices[1] = globalIndices[1];
  broydenWorkGraph.InsertGlobalIndices( globalIndices[1], 2, &myGlobalIndices[0] );
  // Row 3 structure
  myGlobalIndices[0] = globalIndices[1];
  myGlobalIndices[1] = globalIndices[2];
  broydenWorkGraph.InsertGlobalIndices( globalIndices[2], 2, &myGlobalIndices[0] );

  broydenWorkGraph.FillComplete();

  Teuchos::RCP<Epetra_CrsMatrix> broydenWorkMatrix =
    Teuchos::rcp( new Epetra_CrsMatrix( Copy, broydenWorkGraph ) );

  // Create an identity matrix
  broydenWorkVec.PutScalar(1.0);
  broydenWorkMatrix->ReplaceDiagonalValues(broydenWorkVec);

  NOX::Epetra::BroydenOperator broydenOp( noxParams, printing, broydenWorkVec, broydenWorkMatrix, true );

  broydenWorkVec[0] =  1.0;
  broydenWorkVec[1] = -1.0;
  broydenWorkVec[2] =  2.0;
  broydenOp.setStepVector( broydenWorkVec );

  broydenWorkVec[0] =  2.0;
  broydenWorkVec[1] =  1.0;
  broydenWorkVec[2] =  3.0;
  broydenOp.setYieldVector( broydenWorkVec );

  broydenOp.computeSparseBroydenUpdate();

  // Create the gold matrix for comparison
  Teuchos::RCP<Epetra_CrsMatrix> goldMatrix = Teuchos::rcp( new Epetra_CrsMatrix( Copy, broydenWorkGraph ) );

  int      numCols ;
  double * values  ;

  // Row 1 answers
  goldMatrix->ExtractMyRowView( 0, numCols, values );
  values[0] =  6.0 ;
  values[1] =  2.0 ;
  // Row 2 answers
  goldMatrix->ExtractMyRowView( 1, numCols, values );
  values[0] =  5.0 ;
  values[1] =  0.0 ;
  // Row 3 structure
  goldMatrix->ExtractMyRowView( 2, numCols, values );

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

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

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

void 
HMX_PDE::computeSourceTerm(map<string, Epetra_Vector*> fields, 
                                Epetra_Vector& result)
{
  // All dependent variables should be in place before now

  // We assume all other registered dependent problems are species which
  // affect the reaction rate of this specie

  Epetra_Vector *TvecPtr = 0;
  map<string, Epetra_Vector*>::iterator iter = fields.find(tempFieldName);
  if( iter == fields.end() ) 
  {
    std::cout << "ERROR: Cannot find Temperature field \"" << tempFieldName 
         << "\" for use in computeSourceTerm for problem \""
         << getName() << std::endl;
    throw "HMX_PDE ERROR";
  }
  else
    TvecPtr = (*iter).second;

  Epetra_Vector &T = *TvecPtr;

  // If this problem is the temperature equation, don't compute a source
  // term.  This would be where a volumetric heating term would go.
  if( getName() == tempFieldName ) 
  {
    result.PutScalar(0.0);
    return;
  }
  else 
  {
    double rateK;
    
    map<string, double>::iterator requiredFieldIter;
    map<string, double>::iterator requiredFieldEnd = SrcTermExponent.end();

    for( int i = 0; i<result.MyLength(); i++ ) {

      rateK = pow(T[i],StericCoef) * PreExp_A * 
              exp( -ActEnergy / (Const_R * T[i]) );

      result[i] = 1.0;  // Start point for product

      // Loop over all required fields and contribute to product
      for( requiredFieldIter = SrcTermExponent.begin();
           requiredFieldIter != requiredFieldEnd; requiredFieldIter++) 
      {
        iter = fields.find( (*requiredFieldIter).first );
        if( iter == fields.end() ) 
        {
          std::cout << "ERROR: Cannot find required field \"" 
               << (*requiredFieldIter).first 
               << "\" for use in computeSourceTerm for problem \""
               << getName() << std::endl;
          throw "HMX_PDE ERROR";
        }
	Epetra_Vector &reqFieldVec = *((*iter).second);

        result[i] *= pow( reqFieldVec[i], (*requiredFieldIter).second );
      }
      result[i] *= rateK;
    }
  }
}

int test_azoo_scaling(Epetra_CrsMatrix& A,
                      Epetra_Vector& x,
                      Epetra_Vector& b,
                      bool verbose)
{
  Epetra_Vector vec1(x);
  Epetra_Vector vec2(x);
  Epetra_Vector diag(x);
  Epetra_Vector vec3(x);
  Epetra_Vector vec4(x);
  Epetra_Vector rhs(x);
  Epetra_Vector soln_none(x);
  Epetra_Vector soln_jacobi(x);
  Epetra_Vector soln_rowsum(x);
  Epetra_Vector soln_symdiag(x);

  vec1.PutScalar(1.0);

  A.Multiply(false, vec1, vec2);

  A.ExtractDiagonalCopy(diag);

  double* diag_vals = NULL;
  diag.ExtractView(&diag_vals);

  int* options = new int[AZ_OPTIONS_SIZE];
  double* params = new double[AZ_PARAMS_SIZE];
  AZ_defaults(options, params);

  options[AZ_output] = verbose ? 1 : AZ_none;

  options[AZ_scaling] = AZ_Jacobi;
  AztecOO::MatrixData mdata(&A);
  AZ_MATRIX* Amat = AZ_matrix_create(vec1.Map().NumMyElements());
  AZ_set_MATFREE(Amat, (void*)(&mdata), Epetra_Aztec_matvec);

  AZ_SCALING* scaling = AZ_scaling_create();

  double* xvals = NULL, *bvals = NULL;
  x.ExtractView(&xvals);
  b.ExtractView(&bvals);

  int err = AztecOO_scale_epetra(AZ_SCALE_MAT_RHS_SOL, Amat,
                                 options, bvals, xvals, NULL, scaling);
  if (err != 0) {
    if (verbose) {
      cout << "AztecOO_scale_epetra returned err="<<err<<endl;
    }
    return(err);
  }

  A.Multiply(false, vec1, vec3);

  vec4.Multiply(1.0, diag, vec3, 0.0);

  double vec2nrm, vec4nrm;

  vec2.Norm2(&vec2nrm);
  vec4.Norm2(&vec4nrm);

  if (fabs(vec2nrm - vec4nrm) > 1.e-6) {
    return(-1);
  }

  //now call the scaling function again, just to allow for
  //testing memory-leak issues.
  err = AztecOO_scale_epetra(AZ_SCALE_MAT_RHS_SOL, Amat,
                             options, bvals, xvals, NULL, scaling);
  if (err != 0) {
    if (verbose) {
      cout << "AztecOO_scale_epetra returned err="<<err<<endl;
    }
    return(err);
  }

  AztecOO_scale_epetra(AZ_DESTROY_SCALING_DATA, Amat, options,
                       bvals, xvals, NULL, scaling);

  x.PutScalar(1.0);

  Epetra_CrsMatrix* Atmp = create_and_fill_crs_matrix(A.RowMap());
  Atmp->Multiply(false, x, rhs);

  x.PutScalar(0.0);

  AztecOO azoo(&A, &x, &b);

  azoo.SetAztecOption(AZ_scaling, AZ_Jacobi);
  if (verbose) {
    azoo.SetAztecOption(AZ_output, 1);
  }
  else {
    azoo.SetAztecOption(AZ_output, AZ_none);
  }

  azoo.Iterate(100, 1.e-6);

  delete Atmp;

  Epetra_CrsMatrix* Atmp1 = create_and_fill_crs_matrix(A.RowMap());
//.........这里部分代码省略.........

int test_azoo_with_ilut(Epetra_CrsMatrix& A,
                        Epetra_Vector& x,
                        Epetra_Vector& b,
                        bool verbose)
{
  x.PutScalar(0.0);

  AztecOO* azoo0 = new AztecOO(&A, &x, &b);

  azoo0->SetAztecOption(AZ_solver, AZ_gmres);
  azoo0->SetAztecOption(AZ_output, AZ_none);
  azoo0->SetAztecOption(AZ_conv, AZ_Anorm);
  azoo0->SetAztecOption(AZ_precond, AZ_dom_decomp);
  azoo0->SetAztecOption(AZ_subdomain_solve, AZ_ilut);
  azoo0->SetAztecOption(AZ_keep_info, 1);

  if (verbose) {
    cout << "testing AztecOO with GMRES and ILUT, AZ_keep_info==1" << endl;
  }

  int maxiters = 100;
  double tolerance = 1.e-12;
  int err = azoo0->Iterate(maxiters, tolerance);
  if (err != 0) {
    if (verbose) cout << "AztecOO::Iterate return err="<<err<<endl;
    return(err);
  }

  double resid = resid2norm(A, x, b);
  if (verbose) {
    cout << "residual 2-norm after GMRES/ILUT solve: "
        << resid << endl;
  }

  if (resid > 1.e-6) {
    return(-1);
  }

  if (verbose) {
    cout << "solving with GMRES/ILUT again, AZ_pre_calc==AZ_reuse"
        << endl << "(will error out if factors weren't kept from"
         << " previous solve)"<<endl;
  }

  azoo0->SetAztecOption(AZ_pre_calc, AZ_reuse);
  x.PutScalar(0.0);
  err = azoo0->Iterate(maxiters, tolerance);
  if (err != 0) {
    if (verbose) cout << "AztecOO::Iterate return err="<<err<<endl;
    return(err);
  }

  double resid2 = resid2norm(A, x, b);
  if (verbose) {
    cout << "after second GMRES/ILUT solve, residual 2-norm: "
      << resid2 <<endl;
  }
  if (resid2 > 1.e-6) {
    return(-1);
  }

  if (verbose) {
    cout << "azoo0->SolveTime(): " << azoo0->SolveTime() << endl;
  }

  delete azoo0;

  return(0);
}

int test_azoo_as_precond_op(Epetra_CrsMatrix& A,
                            Epetra_Vector& x,
                            Epetra_Vector& b,
                            bool verbose)
{
  AztecOO* azoo0 = new AztecOO(&A, &x, &b);

  azoo0->SetAztecOption(AZ_solver, AZ_cg);
  azoo0->SetAztecOption(AZ_output, AZ_none);
  azoo0->SetAztecOption(AZ_conv, AZ_none);

  AztecOO_Operator azooOp(azoo0, 10);

  Epetra_LinearProblem elp(&A, &x, &b);

  x.PutScalar(0.0);
  AztecOO* azoo1 = new AztecOO(elp);
  azoo1->SetPrecOperator(&azooOp);
  azoo1->SetAztecOption(AZ_solver, AZ_gmres);
  azoo1->SetAztecOption(AZ_output, AZ_none);
  azoo1->SetAztecOption(AZ_conv, AZ_none);

  if (verbose) {
    cout << "testing recursive solve (AztecOO as"
    << " preconditioner for another AztecOO)."<<endl;
  }

  int maxiters = 100;
  double tolerance = 1.e-12;
  azoo1->Iterate(maxiters, tolerance);

  double resid1 = resid2norm(A, x, b);
  if (verbose) {
    cout << "residual 2-norm after recursive solve: "
        << resid1 << endl;
  }

  if (resid1 > 1.e-6) {
    return(-1);
  }

  if (verbose) {
    cout << "now make sure the precond AztecOO instance"
      << " hasn't been corrupted."<<endl;
  }

//  AZ_manage_memory(0, -43, 0, 0, 0);

  x.PutScalar(0.0);
  azoo0->Iterate(&A, &x, &b, maxiters, tolerance);

  double resid0 = resid2norm(A, x, b);
  if (verbose) {
    cout << "residual 2-norm: " << resid0 << endl;
  }
  if (resid0 > 1.e-6) {
    return(-1);
  }

  delete azoo0;
  delete azoo1;

  return(0);
}

//.........这里部分代码省略.........
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }
   
   
   // set the smoother
   if (level_==0)
      Set_Smoother(ml,ag,level_,nlevel,thislevel_ml_,thislevel_ag_,fsmoothertype,nsmooth_fine);

   else if (level_ != nlevel_-1) // set the smoother from the input
      Set_Smoother(ml,ag,level_,nlevel,thislevel_ml_,thislevel_ag_,smoothertype,nsmooth);

   else // set the coarse solver from the input
      Set_Smoother(ml,ag,level_,nlevel,thislevel_ml_,thislevel_ag_,coarsesolvetype,nsmooth_coarse);
  
   // create this level's preconditioner class
   ML_Epetra::MultiLevelOperator* ml_tmp = new ML_Epetra::MultiLevelOperator(
                                                       thislevel_ml_,comm_,
                                                       Mat->OperatorDomainMap(),
                                                       Mat->OperatorRangeMap());
   thislevel_prec_ = new ML_NOX::ML_Nox_ConstrainedMultiLevelOperator(ml_tmp,*coarseinterface_);

   if (!thislevel_prec_)
   {
      cout << "**ERR**: ML_NOX::ML_Nox_NonlinearLevel::ML_Nox_NonlinearLevel:\n"
           << "**ERR**: thislevel_prec_==NULL on level " << level_ << "\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }
                                                                 
   // intensive test of this level's ML-smoother
#if 0
   {
   cout << "Test of smoother on level " << level_ << endl;
   Epetra_Vector *out = new Epetra_Vector(Copy,*xthis_,0);
   out->PutScalar(0.0);
   cout << "Input\n";
   xthis_->PutScalar(1.0);
   Mat->Multiply(false,*xthis_,*out);
   xthis_->PutScalar(3.0);
   cout << "rhs\n";
   cout << *out;
   double norm = 0.0;
   out->Norm1(&norm);
   cout << "Norm of rhs = " << norm << endl;
   thislevel_prec_->ApplyInverse(*out,*xthis_);
   cout << "result after smoother\n";
   cout << *xthis_;
   delete out; out = 0;
   }
   if (level_==2) exit(0);
#endif   

   // ------------------------------------------------------------------------
   // generate this level's coarse prepostoperator
   if (level_==0)
      coarseprepost_ = new ML_NOX::Ml_Nox_CoarsePrePostOperator(*coarseinterface_,
                                                                fineinterface_);    

   // ------------------------------------------------------------------------
   // set up NOX on this level   
   // ------------------------------------------------------------------------
   nlParams_ = new Teuchos::ParameterList();
   Teuchos::ParameterList& printParams = nlParams_->sublist("Printing");        
   printParams.setParameter("MyPID", comm_.MyPID()); 
   printParams.setParameter("Output Precision", 9);
   printParams.setParameter("Output Processor", 0);
   if (ml_printlevel_>9)

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