C++ Teuchos::rcp方法代码示例

TEUCHOS_UNIT_TEST(tUniqueGlobalIndexer_Utilities,ArrayToFieldVector)
{
   typedef Intrepid::FieldContainer<int> IntFieldContainer;

   // build global (or serial communicator)
   #ifdef HAVE_MPI
      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
   #else
      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
   #endif

   int myRank = eComm->MyPID();
   int numProcs = eComm->NumProc();

   TEUCHOS_ASSERT(numProcs==2);

   RCP<panzer::UniqueGlobalIndexer<short,int> > globalIndexer 
         = rcp(new panzer::unit_test::UniqueGlobalIndexer(myRank,numProcs));

   panzer::ArrayToFieldVector<short,int,KokkosClassic::DefaultNode::DefaultNodeType> atfv(globalIndexer);

   std::map<std::string,IntFieldContainer> dataU, dataT;
   fillFieldContainer(globalIndexer->getFieldNum("U"),"block_0",*globalIndexer,dataU["block_0"]);
   fillFieldContainer(globalIndexer->getFieldNum("U"),"block_1",*globalIndexer,dataU["block_1"]);
   fillFieldContainer(globalIndexer->getFieldNum("T"),"block_0",*globalIndexer,dataT["block_0"]);

   Teuchos::RCP<Tpetra::MultiVector<int,int,int> > reducedUDataVector = atfv.getDataVector<int>("U",dataU);
   Teuchos::RCP<Tpetra::MultiVector<int,int,int> > reducedTDataVector = atfv.getDataVector<int>("T",dataT);

   std::vector<int> fields_u(reducedUDataVector->getLocalLength());
   std::vector<int> fields_t(reducedTDataVector->getLocalLength());

   reducedUDataVector->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(fields_u));
   reducedTDataVector->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(fields_t));
   
   if(myRank==0) {
      TEST_EQUALITY(reducedUDataVector->getLocalLength(),6);
      TEST_EQUALITY(reducedTDataVector->getLocalLength(),5);

      TEST_EQUALITY(fields_u[0], 6);
      TEST_EQUALITY(fields_u[1], 0);
      TEST_EQUALITY(fields_u[2], 2);
      TEST_EQUALITY(fields_u[3], 8);
      TEST_EQUALITY(fields_u[4],10);
      TEST_EQUALITY(fields_u[5],13);

      TEST_EQUALITY(fields_t[0], 7);
      TEST_EQUALITY(fields_t[1], 1);
      TEST_EQUALITY(fields_t[2], 3);
      TEST_EQUALITY(fields_t[3], 9);
      TEST_EQUALITY(fields_t[4],11);
   }
   else if(myRank==1) {
      TEST_EQUALITY(reducedUDataVector->getLocalLength(),4);
      TEST_EQUALITY(reducedTDataVector->getLocalLength(),1);

      TEST_EQUALITY(fields_u[0], 4);
      TEST_EQUALITY(fields_u[1],12);
      TEST_EQUALITY(fields_u[2],15);
      TEST_EQUALITY(fields_u[3],14);

      TEST_EQUALITY(fields_t[0], 5);
   }
   else 
      TEUCHOS_ASSERT(false);
}

//----------------------------------------------------------------------------
Teuchos::RCP<AAdapt::AbstractAdapter>
AAdapt::AdaptationFactory::createAdapter() {
  using Teuchos::rcp;

  Teuchos::RCP<AAdapt::AbstractAdapter> strategy;
  std::string& method = adapt_params_->get("Method", "");

#if defined(ALBANY_STK)
  if(method == "Copy Remesh") {
    strategy = rcp(new AAdapt::CopyRemesh(adapt_params_,
                                          param_lib_,
                                          state_mgr_,
                                          commT_));
  }

#if defined(ALBANY_LCM) && defined(ALBANY_BGL)

  else if(method == "Topmod") {
    strategy = rcp(new AAdapt::TopologyMod(adapt_params_,
                                           param_lib_,
                                           state_mgr_,
                                           commT_));
  }

#endif

#if defined(ALBANY_LCM) && defined(LCM_SPECULATIVE)

  else if(method == "Random") {
    strategy = rcp(new AAdapt::RandomFracture(adapt_params_,
                   param_lib_,
                   state_mgr_,
                   commT_));
  }

#endif
#endif
#if 0
#if defined(ALBANY_LCM) && defined(ALBANY_STK_PERCEPT)

  else if(method == "Unif Size") {
    strategy = rcp(new AAdapt::STKAdapt<AAdapt::STKUnifRefineField>(adapt_params_,
                   param_lib_,
                   state_mgr_,
                   epetra_comm_));
  }

#endif
#endif

#if defined(ALBANY_STK)
  else 
#endif
    TEUCHOS_TEST_FOR_EXCEPTION(true,
                               Teuchos::Exceptions::InvalidParameter,
                               std::endl <<
                               "Error! Unknown adaptivity method requested:"
                               << method <<
                               " !" << std::endl
                               << "Supplied parameter list is " <<
                               std::endl << *adapt_params_);

  return strategy;
}

int main(int argc, char *argv[]) {
  //
  int MyPID = 0;
#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
  MyPID = Comm.MyPID();
#else
  Epetra_SerialComm Comm;
#endif
  //
  typedef double                            ST;
  typedef Teuchos::ScalarTraits<ST>        SCT;
  typedef SCT::magnitudeType                MT;
  typedef Epetra_MultiVector                MV;
  typedef Epetra_Operator                   OP;
  typedef Belos::MultiVecTraits<ST,MV>     MVT;
  typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  bool verbose = false, debug = false, proc_verbose = false;
  int frequency = -1;        // frequency of status test output.
  int blocksize = 1;         // blocksize
  int numrhs = 1;            // number of right-hand sides to solve for
  int maxiters = -1;         // maximum number of iterations allowed per linear system
  int maxsubspace = 50;      // maximum number of blocks the solver can use for the subspace
  int maxrestarts = 15;      // number of restarts allowed 
  std::string filename("orsirr1.hb");
  MT tol = 1.0e-5;           // relative residual tolerance

  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("debug","nondebug",&debug,"Print debugging information from solver.");
  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
  cmdp.setOption("filename",&filename,"Filename for test matrix.  Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
  cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
  cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
  cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }
  if (!verbose)
    frequency = -1;  // reset frequency if test is not verbose
  //
  // Get the problem
  //
  RCP<Epetra_Map> Map;
  RCP<Epetra_CrsMatrix> A;
  RCP<Epetra_MultiVector> B, X;
  RCP<Epetra_Vector> vecB, vecX;
  EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB);
  A->OptimizeStorage();
  proc_verbose = verbose && (MyPID==0);  /* Only print on the zero processor */

  // Check to see if the number of right-hand sides is the same as requested.
  if (numrhs>1) {
    X = rcp( new Epetra_MultiVector( *Map, numrhs ) );
    B = rcp( new Epetra_MultiVector( *Map, numrhs ) );
    X->Seed();
    X->Random();
    OPT::Apply( *A, *X, *B );
    X->PutScalar( 0.0 );
  }
  else {
    X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
    B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
  }
  //
  // ********Other information used by block solver***********
  // *****************(can be user specified)******************
  //
  const int NumGlobalElements = B->GlobalLength();
  if (maxiters == -1)
    maxiters = NumGlobalElements/blocksize - 1; // maximum number of iterations to run
  //
  ParameterList belosList;
  belosList.set( "Num Blocks", maxsubspace);             // Maximum number of blocks in Krylov factorization
  belosList.set( "Block Size", blocksize );              // Blocksize to be used by iterative solver
  belosList.set( "Maximum Iterations", maxiters );       // Maximum number of iterations allowed
  belosList.set( "Maximum Restarts", maxrestarts );      // Maximum number of restarts allowed
  belosList.set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested
  int verbosity = Belos::Errors + Belos::Warnings;
  if (verbose) {
    verbosity += Belos::TimingDetails + Belos::StatusTestDetails;
    if (frequency > 0)
      belosList.set( "Output Frequency", frequency );
  }
  if (debug) {
    verbosity += Belos::Debug;
  }
  belosList.set( "Verbosity", verbosity );
  //
  // Construct an unpreconditioned linear problem instance.
//.........这里部分代码省略.........

/*! \brief Create a mesh of approximately the desired size.
 *
 *  We want 3 dimensions close to equal in length.
 */
const RCP<tMVector_t> getMeshCoordinates(
    const RCP<const Teuchos::Comm<int> > & comm,
    zgno_t numGlobalCoords)
{
  int rank = comm->getRank();
  int nprocs = comm->getSize();

  double k = log(numGlobalCoords) / 3;
  double xdimf = exp(k) + 0.5;
  ssize_t xdim = static_cast<ssize_t>(floor(xdimf));
  ssize_t ydim = xdim;
  ssize_t zdim = numGlobalCoords / (xdim*ydim);
  ssize_t num=xdim*ydim*zdim;
  ssize_t diff = numGlobalCoords - num;
  ssize_t newdiff = 0;

  while (diff > 0){
    if (zdim > xdim && zdim > ydim){
      zdim++;
      newdiff = diff - (xdim*ydim);
      if (newdiff < 0)
        if (diff < -newdiff)
          zdim--;
    }
    else if (ydim > xdim && ydim > zdim){
      ydim++;
      newdiff = diff - (xdim*zdim);
      if (newdiff < 0)
        if (diff < -newdiff)
          ydim--;
    }
    else{
      xdim++;
      newdiff = diff - (ydim*zdim);
      if (newdiff < 0)
        if (diff < -newdiff)
          xdim--;
    }

    diff = newdiff;
  }

  num=xdim*ydim*zdim;
  diff = numGlobalCoords - num;
  if (diff < 0)
    diff /= -numGlobalCoords;
  else
    diff /= numGlobalCoords;

  if (rank == 0){
    if (diff > .01)
      cout << "Warning: Difference " << diff*100 << " percent" << endl;
    cout << "Mesh size: " << xdim << "x" << ydim << "x" <<
      zdim << ", " << num << " vertices." << endl;
  }

  // Divide coordinates.

  ssize_t numLocalCoords = num / nprocs;
  ssize_t leftOver = num % nprocs;
  ssize_t gid0 = 0;

  if (rank <= leftOver)
    gid0 = zgno_t(rank) * (numLocalCoords+1);
  else
    gid0 = (leftOver * (numLocalCoords+1)) + 
           ((zgno_t(rank) - leftOver) * numLocalCoords);

  if (rank < leftOver)
    numLocalCoords++;

  ssize_t gid1 = gid0 + numLocalCoords;

  zgno_t *ids = new zgno_t [numLocalCoords];
  if (!ids)
    throw bad_alloc();
  ArrayRCP<zgno_t> idArray(ids, 0, numLocalCoords, true);

  for (ssize_t i=gid0; i < gid1; i++)
    *ids++ = zgno_t(i);   

  RCP<const tMap_t> idMap = rcp(
    new tMap_t(num, idArray.view(0, numLocalCoords), 0, comm));

  // Create a Tpetra::MultiVector of coordinates.

  zscalar_t *x = new zscalar_t [numLocalCoords*3]; 
  if (!x)
    throw bad_alloc();
  ArrayRCP<zscalar_t> coordArray(x, 0, numLocalCoords*3, true);

  zscalar_t *y = x + numLocalCoords;
  zscalar_t *z = y + numLocalCoords;

  zgno_t xStart = 0;
  zgno_t yStart = 0;
//.........这里部分代码省略.........

bool run_composite_linear_ops_tests(
  const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > comm,
  const int n,
  const bool useSpmd,
  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType &tol,
  const bool dumpAll,
  Teuchos::FancyOStream *out_arg
  )
{

  using Teuchos::as;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType    ScalarMag;
  typedef Teuchos::ScalarTraits<ScalarMag> STM;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::null;
  using Teuchos::rcp_const_cast;
  using Teuchos::rcp_dynamic_cast;
  using Teuchos::dyn_cast;
  using Teuchos::OSTab;
  using Thyra::relErr;
  using Thyra::passfail;

  RCP<Teuchos::FancyOStream>
    out = rcp(new Teuchos::FancyOStream(rcp(out_arg,false)));

  const Teuchos::EVerbosityLevel
    verbLevel = dumpAll?Teuchos::VERB_EXTREME:Teuchos::VERB_HIGH;

  if (nonnull(out)) *out
    << "\n*** Entering run_composite_linear_ops_tests<"<<ST::name()<<">(...) ...\n";

  bool success = true, result;

  const ScalarMag warning_tol = ScalarMag(1e-2)*tol, error_tol = tol;
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.linear_properties_warning_tol(warning_tol);
  linearOpTester.linear_properties_error_tol(error_tol);
  linearOpTester.adjoint_warning_tol(warning_tol);
  linearOpTester.adjoint_error_tol(error_tol);
  linearOpTester.dump_all(dumpAll);
  Thyra::LinearOpTester<Scalar> symLinearOpTester(linearOpTester);
  symLinearOpTester.check_for_symmetry(true);
  symLinearOpTester.symmetry_warning_tol(STM::squareroot(warning_tol));
  symLinearOpTester.symmetry_error_tol(STM::squareroot(error_tol));

  RCP<const Thyra::VectorSpaceBase<Scalar> > space;
  if(useSpmd) space = Thyra::defaultSpmdVectorSpace<Scalar>(comm,n,-1);
  else space = Thyra::defaultSpmdVectorSpace<Scalar>(n);
  if (nonnull(out)) *out
    << "\nUsing a basic vector space described as " << describe(*space,verbLevel) << " ...\n";

  if (nonnull(out)) *out << "\nCreating random n x (n/2) multi-vector origA ...\n";
  RCP<Thyra::MultiVectorBase<Scalar> >
    mvOrigA = createMembers(space,n/2,"origA");
  Thyra::seed_randomize<Scalar>(0);
  //RTOpPack::show_spmd_apply_op_dump = true;
  Thyra::randomize( as<Scalar>(as<Scalar>(-1)*ST::one()), as<Scalar>(as<Scalar>(+1)*ST::one()),
    mvOrigA.ptr() );
  RCP<const Thyra::LinearOpBase<Scalar> >
    origA = mvOrigA;
  if (nonnull(out)) *out << "\norigA =\n" << describe(*origA,verbLevel);
  //RTOpPack::show_spmd_apply_op_dump = false;

  if (nonnull(out)) *out << "\nTesting origA ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.check(*origA, out.ptr());
  if(!result) success = false;

  if (nonnull(out)) *out
    << "\nCreating implicit scaled linear operator A1 = scale(0.5,origA) ...\n";
  RCP<const Thyra::LinearOpBase<Scalar> >
    A1 = scale(as<Scalar>(0.5),origA);
  if (nonnull(out)) *out << "\nA1 =\n" << describe(*A1,verbLevel);

  if (nonnull(out)) *out << "\nTesting A1 ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.check(*A1,out.ptr());
  if(!result) success = false;

  if (nonnull(out)) *out << "\nTesting that A1.getOp() == origA ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.compare(
    *dyn_cast<const Thyra::DefaultScaledAdjointLinearOp<Scalar> >(*A1).getOp(),
    *origA,out.ptr());
  if(!result) success = false;

  {

    if (nonnull(out)) *out
      << "\nUnwrapping origA to get non-persisting pointer to origA_1, scalar and transp ...\n";
    Scalar  scalar;
    Thyra::EOpTransp transp;
    const Thyra::LinearOpBase<Scalar> *origA_1 = NULL;
    unwrap( *origA, &scalar, &transp, &origA_1 );
    TEUCHOS_TEST_FOR_EXCEPT( origA_1 == NULL );

    if (nonnull(out)) *out << "\nscalar = " << scalar << " == 1 ? ";
    result = (scalar == ST::one());
//.........这里部分代码省略.........

void
Albany::ResponseFactory::
createResponseFunction(
  const std::string& name,
  Teuchos::ParameterList& responseParams,
  Teuchos::Array< Teuchos::RCP<AbstractResponseFunction> >& responses) const
{
  using std::string;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;
  using Teuchos::Array;

  RCP<const Epetra_Comm> comm = app->getComm();

  if (name == "Solution Average") {
    responses.push_back(rcp(new Albany::SolutionAverageResponseFunction(comm)));
  }

  else if (name == "Solution Two Norm") {
    responses.push_back(rcp(new Albany::SolutionTwoNormResponseFunction(comm)));
  }

  else if (name == "Solution Values") {
    responses.push_back(rcp(new Albany::SolutionValuesResponseFunction(app, responseParams)));
  }

  else if (name == "Solution Max Value") {
    int eq = responseParams.get("Equation", 0);
    int neq = responseParams.get("Num Equations", 1);
    bool inor =  responseParams.get("Interleaved Ordering", true);

    responses.push_back(
      rcp(new Albany::SolutionMaxValueResponseFunction(comm, neq, eq, inor)));
  }

  else if (name == "Solution Two Norm File") {
    responses.push_back(
      rcp(new Albany::SolutionFileResponseFunction<Albany::NormTwo>(comm)));
  }

  else if (name == "Solution Inf Norm File") {
    responses.push_back(
      rcp(new Albany::SolutionFileResponseFunction<Albany::NormInf>(comm)));
  }

  else if (name == "Aggregated") {
    int num_aggregate_responses = responseParams.get<int>("Number");
    Array< RCP<AbstractResponseFunction> > aggregated_responses;
    Array< RCP<ScalarResponseFunction> > scalar_responses;
    for (int i=0; i<num_aggregate_responses; i++) {
      std::string id = Albany::strint("Response",i);
      std::string name = responseParams.get<std::string>(id);
      std::string sublist_name = Albany::strint("ResponseParams",i);
      createResponseFunction(name, responseParams.sublist(sublist_name),
			     aggregated_responses);

    }
    scalar_responses.resize(aggregated_responses.size());
    for (int i=0; i<aggregated_responses.size(); i++) {
      TEUCHOS_TEST_FOR_EXCEPTION(
	aggregated_responses[i]->isScalarResponse() != true, std::logic_error,
	"Response function " << i << " is not a scalar response function." <<
	std::endl <<
	"The aggregated response can only aggregate scalar response " <<
	"functions!");
      scalar_responses[i] =
	Teuchos::rcp_dynamic_cast<ScalarResponseFunction>(
	  aggregated_responses[i]);
    }
    responses.push_back(
      rcp(new Albany::AggregateScalarResponseFunction(comm, scalar_responses)));
  }

  else if (name == "Field Integral" ||
	   name == "Field Value" ||
	   name == "Field Average" ||
	   name == "Surface Velocity Mismatch" ||
	   name == "Center Of Mass" ||
	   name == "Save Field" ||
	   name == "Region Boundary" ||
	   name == "Element Size Field" ||
	   name == "IP to Nodal Field" ||
	   name == "Save Nodal Fields" ||
	   name == "PHAL Field Integral") {
    responseParams.set("Name", name);
    for (int i=0; i<meshSpecs.size(); i++) {
      responses.push_back(
	rcp(new Albany::FieldManagerScalarResponseFunction(
	      app, prob, meshSpecs[i], stateMgr, responseParams)));
    }
  }

  else if (name == "Solution") {
    responses.push_back(
      rcp(new Albany::SolutionResponseFunction(app, responseParams)));
  }

  else if (name == "KL") {
    Array< RCP<AbstractResponseFunction> > base_responses;
//.........这里部分代码省略.........

int Ifpack_SORa::Compute(){
  if(!IsInitialized_) Initialize();
  Epetra_Map *RowMap=const_cast<Epetra_Map*>(&A_->RowMatrixRowMap());
  Epetra_Vector Adiag(*RowMap);
  Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
  int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
  double *vals_s,*vals_h;
  bool RowMatrixMode=(Acrs_==Teuchos::null);

  // Label
  sprintf(Label_, "IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_); 

  if(RowMatrixMode){
    if(!A_->Comm().MyPID()) cout<<"SORa: RowMatrix mode enabled"<<endl;
    // RowMatrix mode, build Acrs_ the hard way.
    Epetra_RowMatrix *Arow=&*A_;
    Epetra_Map *ColMap=const_cast<Epetra_Map*>(&A_->RowMatrixColMap());

    int Nmax=A_->MaxNumEntries();
    int length;
    vector<int> indices(Nmax);
    vector<double> values(Nmax); 
    Epetra_CrsMatrix *Acrs=new Epetra_CrsMatrix(Copy,*RowMap,Nmax);

    for(int i=0;i<Arow->NumMyRows();i++){
      Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices[0]);
      for(int j=0;j<length;j++)
	indices[j]=ColMap->GID(indices[j]);
      Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices[0]);
    }
    Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
    Acrs_=rcp(Acrs,true);
  }

  // Create Askew2
  // Note: Here I let EpetraExt do the thinking for me.  Since it gets all the maps correct for the E + F^T stencil.
  // There are probably more efficient ways to do this but this method has the bonus of allowing code reuse.
  IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,-1,Askew2));
  Askew2->FillComplete();
  
  // Create Aherm2
  IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,1,Aherm2));
  Aherm2->FillComplete();

  int nnz2=Askew2->NumMyNonzeros();
  int N=Askew2->NumMyRows();


  // Grab pointers
  IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
  IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));

  // Sanity checking: Make sure the sparsity patterns are the same
#define SANITY_CHECK
#ifdef SANITY_CHECK
  for(int i=0;i<N;i++)
    if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
  for(int i=0;i<nnz2;i++)
    if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
#endif

  // Dirichlet Detection & Nuking of Aherm2 and Askew2
  // Note: Relies on Aherm2/Askew2 having identical sparsity patterns (see sanity check above)
  if(HaveOAZBoundaries_){ 
    int numBCRows;
    int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
    Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
    Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
    Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
    delete [] dirRows;
    delete dirCols;
  }

  // Grab diagonal of A
  A_->ExtractDiagonalCopy(Adiag);

  // Allocate the diagonal for W
  Epetra_Vector *Wdiag = new Epetra_Vector(*RowMap);

  // Build the W matrix (lower triangle only)
  // Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
  int maxentries=Askew2->MaxNumEntries();
  int* gids=new int [maxentries];
  double* newvals=new double[maxentries];
  W=new Epetra_CrsMatrix(Copy,*RowMap,0);
  for(int i=0;i<N;i++){
    // Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
    int rowgid=Acrs_->GRID(i);
    double c_data=0.0;
    double ipdamp=0.0;
    int idx=0;

    for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){      
      int colgid=Askew2->GCID(colind_s[j]);
      c_data+=fabs(vals_s[j]);
      if(rowgid>colgid){
	// Rely on the fact that off-diagonal entries are always numbered last, dropping the entry entirely.
	if(colind_s[j] < N) {       
	  gids[idx]=colgid;
	  newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

  int status=0; // 0 = pass, failures are incremented

  // Initialize MPI and timer
  int Proc=0;
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  double total_time = -MPI_Wtime();
  (void) MPI_Comm_rank(MPI_COMM_WORLD, &Proc);
  MPI_Comm appComm = MPI_COMM_WORLD;
#else
  int appComm=0;
#endif

  using Teuchos::RCP;
  using Teuchos::rcp;

  // Command-line argument for input file
  //char* defaultfile="input_1.xml";
  
  try {

    RCP<EpetraExt::ModelEvaluator> Model = rcp(new MockModelEval_A(appComm));

    // Set input arguments to evalModel call
    EpetraExt::ModelEvaluator::InArgs inArgs = Model->createInArgs();

    RCP<Epetra_Vector> x = rcp(new Epetra_Vector(*(Model->get_x_init())));
    inArgs.set_x(x);

    int num_p = inArgs.Np();     // Number of *vectors* of parameters
    RCP<Epetra_Vector> p1;
    if (num_p > 0) {
      p1 = rcp(new Epetra_Vector(*(Model->get_p_init(0))));
      inArgs.set_p(0,p1);
    }
    int numParams = p1->MyLength(); // Number of parameters in p1 vector

    // Set output arguments to evalModel call
    EpetraExt::ModelEvaluator::OutArgs outArgs = Model->createOutArgs();

    RCP<Epetra_Vector> f = rcp(new Epetra_Vector(x->Map()));
    outArgs.set_f(f);

    int num_g = outArgs.Ng(); // Number of *vectors* of responses
    RCP<Epetra_Vector> g1;
    if (num_g > 0) {
      g1 = rcp(new Epetra_Vector(*(Model->get_g_map(0))));
      outArgs.set_g(0,g1);
    }

    RCP<Epetra_Operator> W_op = Model->create_W();
    outArgs.set_W(W_op);

    RCP<Epetra_MultiVector> dfdp = rcp(new Epetra_MultiVector(
                             *(Model->get_x_map()), numParams));
    outArgs.set_DfDp(0, dfdp);

    RCP<Epetra_MultiVector> dgdp = rcp(new Epetra_MultiVector(g1->Map(), numParams));
    outArgs.set_DgDp(0, 0, dgdp);

    RCP<Epetra_MultiVector> dgdx = rcp(new Epetra_MultiVector(x->Map(), g1->MyLength()));
    outArgs.set_DgDx(0, dgdx);

    // Now, evaluate the model!
    Model->evalModel(inArgs, outArgs);

    // Print out everything
   if (Proc == 0)
    cout << "Finished Model Evaluation: Printing everything {Exact in brackets}" 
         << "\n-----------------------------------------------------------------"
         << std::setprecision(9) << endl;
    x->Print(cout << "\nSolution vector! {3,3,3,3}\n");
    if (num_p>0) p1->Print(cout << "\nParameters! {1,1}\n");
    f->Print(cout << "\nResidual! {8,5,0,-7}\n");
    if (num_g>0) g1->Print(cout << "\nResponses! {2}\n");
    RCP<Epetra_CrsMatrix> W = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_op, true);
    W->Print(cout << "\nJacobian! {6 on diags}\n");
    dfdp->Print(cout << "\nDfDp sensitivity MultiVector! {-1,0,0,0}{0,-4,-6,-8}\n");
    dgdp->Print(cout << "\nDgDp response sensitivity MultiVector!{2,2}\n");
    dgdx->Print(cout << "\nDgDx^T response gradient MultiVector! {-2,-2,-2,-2}\n");

    if (Proc == 0)
     cout <<
      "\n-----------------------------------------------------------------\n";
  }

  catch (std::exception& e) {
    cout << e.what() << endl;
    status = 10;
  }
  catch (string& s) {
    cout << s << endl;
    status = 20;
  }
  catch (char *s) {
    cout << s << endl;
    status = 30;
  }
//.........这里部分代码省略.........

  TEUCHOS_UNIT_TEST(initial_condition_control, control)
  {
    using Teuchos::RCP;
    using Teuchos::rcp;

    Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));

    // setup mesh
    /////////////////////////////////////////////
    RCP<panzer_stk_classic::STK_Interface> mesh;
    {
      RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
      pl->set<int>("X Elements",2);
      pl->set<int>("Y Elements",2);
      pl->set<int>("X Blocks",2);
      pl->set<int>("Y Blocks",1);

      panzer_stk_classic::SquareQuadMeshFactory mesh_factory;
      mesh_factory.setParameterList(pl);
      mesh = mesh_factory.buildMesh(MPI_COMM_WORLD);
      mesh->writeToExodus("test.exo");
    }

    RCP<const shards::CellTopology> ct = mesh->getCellTopology("eblock-0_0");

    panzer::CellData cellData(4,ct);

    RCP<panzer::IntegrationRule> int_rule = rcp(new panzer::IntegrationRule(2,cellData));

    ICFieldDescriptor densDesc;
    densDesc.fieldName  = "DENSITY";
    densDesc.basisName  = "Const";
    densDesc.basisOrder = 0;

    ICFieldDescriptor condDesc;
    condDesc.fieldName  = "CONDUCTIVITY";
    condDesc.basisName  = "HGrad";
    condDesc.basisOrder = 1;

    RCP<panzer::PureBasis> const_basis = rcp(new panzer::PureBasis(densDesc.basisName,densDesc.basisOrder,cellData));
    RCP<panzer::PureBasis> hgrad_basis = rcp(new panzer::PureBasis(condDesc.basisName,condDesc.basisOrder,cellData));

    RCP<const panzer::FieldPattern> constFP = rcp(new panzer::Intrepid2FieldPattern(const_basis->getIntrepid2Basis()));
    RCP<const panzer::FieldPattern> hgradFP = rcp(new panzer::Intrepid2FieldPattern(hgrad_basis->getIntrepid2Basis()));
 
    // setup DOF manager
    /////////////////////////////////////////////
    RCP<panzer::ConnManager<int,int> > conn_manager 
           = Teuchos::rcp(new panzer_stk_classic::STKConnManager<int>(mesh));

    RCP<panzer::DOFManager<int,int> > dofManager
        = rcp(new panzer::DOFManager<int,int>(conn_manager,MPI_COMM_WORLD));
    dofManager->addField(densDesc.fieldName, constFP);
    dofManager->addField(condDesc.fieldName, hgradFP);
    dofManager->buildGlobalUnknowns();

    Teuchos::RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > elof 
          = Teuchos::rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(comm.getConst(),dofManager));

    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof = elof;

    // setup worksets
    /////////////////////////////////////////////
    
    std::map<std::string,panzer::WorksetNeeds> needs;
    needs["eblock-0_0"].cellData = cellData;
    needs["eblock-0_0"].int_rules.push_back(int_rule);
    needs["eblock-0_0"].bases          = { const_basis,    hgrad_basis};
    needs["eblock-0_0"].rep_field_name = {   densDesc.fieldName, condDesc.fieldName};

    needs["eblock-1_0"].cellData = cellData;
    needs["eblock-1_0"].int_rules.push_back(int_rule);
    needs["eblock-1_0"].bases          = { const_basis,    hgrad_basis};
    needs["eblock-1_0"].rep_field_name = {   densDesc.fieldName, condDesc.fieldName};

    Teuchos::RCP<panzer_stk_classic::WorksetFactory> wkstFactory 
        = Teuchos::rcp(new panzer_stk_classic::WorksetFactory(mesh)); // build STK workset factory
    Teuchos::RCP<panzer::WorksetContainer> wkstContainer              // attach it to a workset container (uses lazy evaluation)
        = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,needs));

    // setup field manager builder
    /////////////////////////////////////////////
      
    // Add in the application specific closure model factory
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory; 
    user_app::STKModelFactory_TemplateBuilder cm_builder;
    cm_factory.buildObjects(cm_builder);

    Teuchos::ParameterList user_data("User Data");

    Teuchos::ParameterList ic_closure_models("Initial Conditions");
    ic_closure_models.sublist("eblock-0_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-0_0").sublist(condDesc.fieldName).set<double>("Value",9.0);
    ic_closure_models.sublist("eblock-1_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-1_0").sublist(condDesc.fieldName).set<double>("Value",9.0);    

    std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
    block_ids_to_cell_topo["eblock-0_0"] = mesh->getCellTopology("eblock-0_0");
    block_ids_to_cell_topo["eblock-1_0"] = mesh->getCellTopology("eblock-1_0");

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

// =============================================================================
int main (int argc, char *argv[])
{
  Teuchos::GlobalMPISession(&argc, &argv, NULL);

  // Create a communicator for Epetra objects
#ifdef HAVE_MPI
  Epetra_MpiComm eComm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm eComm();
#endif

  const RCP<Teuchos::FancyOStream> out =
    Teuchos::VerboseObjectBase::getDefaultOStream();

  bool success = true;
  try {
    // ===========================================================================
    // handle command line arguments
    Teuchos::CommandLineProcessor My_CLP;

    My_CLP.setDocString("Linear solver testbed for the 1D Poisson matrix.\n");

    std::string action("matvec");
    My_CLP.setOption("action",
                     &action,
                     "Which action to perform with the operator (matvec, solve_cg, solve_minres, solve_gmres)"
                    );

    std::string solver("cg");
//       My_CLP.setOption("solver", &solver, "Krylov subspace method (cg, minres, gmres)");

//       Operator op = JAC;
//       Operator allOpts[] = {JAC, KEO, KEOREG, POISSON1D};
//       std::string allOptNames[] = {"jac", "keo", "keoreg", "poisson1d"};
//       My_CLP.setOption("operator", &op, 4, allOpts, allOptNames);

    bool verbose = true;
    My_CLP.setOption("verbose", "quiet",
                     &verbose,
                     "Print messages and results.");

    int frequency = 10;
    My_CLP.setOption("frequency",
                     &frequency,
                     "Solvers frequency for printing residuals (#iters).");

    // Make sure this value is large enough to keep the cores busy for a while.
    int n = 1000;
    My_CLP.setOption("size",
                     &n,
                     "Size of the equation system (default: 1000).");

    // print warning for unrecognized arguments
    My_CLP.recogniseAllOptions(true);
    My_CLP.throwExceptions(false);

    // finally, parse the command line
    TEUCHOS_ASSERT_EQUALITY(My_CLP.parse (argc, argv),
                            Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL
                           );
    // =========================================================================
    // Construct Epetra matrix.
    RCP<Teuchos::Time> matrixConstructTime =
      Teuchos::TimeMonitor::getNewTimer("Epetra matrix construction");
    RCP<Epetra_CrsMatrix> A;
    {
      Teuchos::TimeMonitor tm(*matrixConstructTime);
      A = contructEpetraMatrix(n, eComm);
    }
//       A->Print(std::cout);

    // create initial guess and right-hand side
    RCP<Epetra_Vector> x =
      rcp(new Epetra_Vector(A->OperatorDomainMap()));
    RCP<Epetra_MultiVector> b =
      rcp(new Epetra_Vector(A->OperatorRangeMap()));
    // b->Random();
    TEUCHOS_ASSERT_EQUALITY(0, b->PutScalar(1.0));

    if (action.compare("matvec") == 0) {
      TEUCHOS_ASSERT_EQUALITY(0, x->PutScalar(1.0));
      RCP<Teuchos::Time> mvTime = Teuchos::TimeMonitor::getNewTimer("Epetra operator apply");
      {
        Teuchos::TimeMonitor tm(*mvTime);
        // Don't TEUCHOS_ASSERT_EQUALITY() here for speed.
        A->Apply(*x, *b);
      }

      // print timing data
      Teuchos::TimeMonitor::summarize();
    } else {
      // -----------------------------------------------------------------------
      // Belos part
      Teuchos::ParameterList belosList;
      // Relative convergence tolerance requested
      belosList.set("Convergence Tolerance", 1.0e-12);
      if (verbose) {
        belosList.set("Verbosity",
                      Belos::Errors +
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifndef HAVE_TPETRA_COMPLEX_DOUBLE
#  error "Anasazi: This test requires Scalar = std::complex<double> to be enabled in Tpetra."
#else
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::tuple;
  using std::cout;
  using std::endl;

  typedef std::complex<double>                ST;
  typedef Teuchos::ScalarTraits<ST>          SCT;
  typedef SCT::magnitudeType                  MT;
  typedef Tpetra::MultiVector<ST>             MV;
  typedef MV::global_ordinal_type             GO;
  typedef Tpetra::Operator<ST>                OP;
  typedef Anasazi::MultiVecTraits<ST,MV>     MVT;
  typedef Anasazi::OperatorTraits<ST,MV,OP>  OPT;

  Teuchos::GlobalMPISession mpisess (&argc, &argv, &std::cout);

  bool success = false;

  const ST ONE = SCT::one ();

  int info = 0;

  RCP<const Teuchos::Comm<int> > comm =
    Tpetra::DefaultPlatform::getDefaultPlatform ().getComm ();

  const int MyPID = comm->getRank ();

  bool verbose = false;
  bool debug = false;
  bool insitu = false;
  bool herm = false;
  std::string which("LM");
  std::string filename;
  int nev = 4;
  int blockSize = 4;
  MT tol = 1.0e-6;

  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
  cmdp.setOption("insitu","exsitu",&insitu,"Perform in situ restarting.");
  cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
  cmdp.setOption("herm","nonherm",&herm,"Solve Hermitian or non-Hermitian problem.");
  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix (assumes non-Hermitian unless specified otherwise).");
  cmdp.setOption("nev",&nev,"Number of eigenvalues to compute.");
  cmdp.setOption("blockSize",&blockSize,"Block size for the algorithm.");
  cmdp.setOption("tol",&tol,"Tolerance for convergence.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }
  if (debug) verbose = true;
  if (filename == "") {
    // get default based on herm
    if (herm) {
      filename = "mhd1280b.cua";
    }
    else {
      filename = "mhd1280a.cua";
    }
  }

  if (MyPID == 0) {
    cout << Anasazi::Anasazi_Version() << endl << endl;
  }

  // Get the data from the HB file
  int dim,dim2,nnz;
  int rnnzmax;
  double *dvals;
  int *colptr,*rowind;
  nnz = -1;
  if (MyPID == 0) {
    info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,&colptr,&rowind,&dvals);
    // find maximum NNZ over all rows
    vector<int> rnnz(dim,0);
    for (int *ri=rowind; ri<rowind+nnz; ++ri) {
      ++rnnz[*ri-1];
    }
    rnnzmax = *std::max_element(rnnz.begin(),rnnz.end());
  }
  else {
    // address uninitialized data warnings
    dvals = NULL;
    colptr = NULL;
    rowind = NULL;
  }
  Teuchos::broadcast(*comm,0,&info);
  Teuchos::broadcast(*comm,0,&nnz);
  Teuchos::broadcast(*comm,0,&dim);
  Teuchos::broadcast(*comm,0,&rnnzmax);
  if (info == 0 || nnz < 0) {
    if (MyPID == 0) {
      cout << "Error reading '" << filename << "'" << endl
           << "End Result: TEST FAILED" << endl;
//.........这里部分代码省略.........

  void IncSVDPOD::Initialize( const Teuchos::RCP< Teuchos::ParameterList >& params,
                              const Teuchos::RCP< const Epetra_MultiVector >& ss,
                              const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > >& fileio ) {

    using Teuchos::rcp;

    // Get the "Reduced Basis Method" sublist.
    Teuchos::ParameterList rbmethod_params = params->sublist( "Reduced Basis Method" );

    // Get the maximum basis size
    maxBasisSize_ = rbmethod_params.get<int>("Max Basis Size");
    TEUCHOS_TEST_FOR_EXCEPTION(maxBasisSize_ < 2,std::invalid_argument,"\"Max Basis Size\" must be at least 2.");

    // Get a filter
    filter_ = rbmethod_params.get<Teuchos::RCP<Filter<double> > >("Filter",Teuchos::null);
    if (filter_ == Teuchos::null) {
      int k = rbmethod_params.get("Rank",(int)5);
      filter_ = rcp( new RangeFilter<double>(LARGEST,k,k) );
    }

    // Get convergence tolerance
    tol_ = rbmethod_params.get<double>("Convergence Tolerance",tol_);

    // Get debugging flag
    debug_ = rbmethod_params.get<bool>("IncSVD Debug",debug_);

    // Get verbosity level
    verbLevel_ = rbmethod_params.get<int>("IncSVD Verbosity Level",verbLevel_);

    // Get an Anasazi orthomanager
    if (rbmethod_params.isType<
          Teuchos::RCP< Anasazi::OrthoManager<double,Epetra_MultiVector> > 
        >("Ortho Manager")
       ) 
    {
      ortho_ = rbmethod_params.get< 
                Teuchos::RCP<Anasazi::OrthoManager<double,Epetra_MultiVector> >
               >("Ortho Manager");
      TEUCHOS_TEST_FOR_EXCEPTION(ortho_ == Teuchos::null,std::invalid_argument,"User specified null ortho manager.");
    }
    else {
      std::string omstr = rbmethod_params.get("Ortho Manager","DGKS");
      if (omstr == "SVQB") {
        ortho_ = rcp( new Anasazi::SVQBOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
      }
      else { // if omstr == "DGKS"
        ortho_ = rcp( new Anasazi::BasicOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
      }
    }

    // Lmin,Lmax,Kstart
    lmin_ = rbmethod_params.get("Min Update Size",1);
    TEUCHOS_TEST_FOR_EXCEPTION(lmin_ < 1 || lmin_ >= maxBasisSize_,std::invalid_argument,
                       "Method requires 1 <= min update size < max basis size.");
    lmax_ = rbmethod_params.get("Max Update Size",maxBasisSize_);
    TEUCHOS_TEST_FOR_EXCEPTION(lmin_ > lmax_,std::invalid_argument,"Max update size must be >= min update size.");

    startRank_ = rbmethod_params.get("Start Rank",lmin_);
    TEUCHOS_TEST_FOR_EXCEPTION(startRank_ < 1 || startRank_ > maxBasisSize_,std::invalid_argument,
                       "Starting rank must be in [1,maxBasisSize_)");
    // MaxNumPasses
    maxNumPasses_ = rbmethod_params.get("Maximum Number Passes",maxNumPasses_);
    TEUCHOS_TEST_FOR_EXCEPTION(maxNumPasses_ != -1 && maxNumPasses_ <= 0, std::invalid_argument,
                       "Maximum number passes must be -1 or > 0.");
    // Save the pointer to the snapshot matrix
    TEUCHOS_TEST_FOR_EXCEPTION(ss == Teuchos::null,std::invalid_argument,"Input snapshot matrix cannot be null.");
    A_ = ss;

    // MaxNumCols
    maxNumCols_ = A_->NumVectors();
    maxNumCols_ = rbmethod_params.get("Maximum Number Columns",maxNumCols_);
    TEUCHOS_TEST_FOR_EXCEPTION(maxNumCols_ < A_->NumVectors(), std::invalid_argument,
                       "Maximum number of columns must be at least as many as in the initializing data set.");

    // V locally replicated or globally distributed
    // this must be true for now
    // Vlocal_ = rbmethod_params.get("V Local",Vlocal_);

    // Allocate space for the factorization
    sigma_.reserve( maxBasisSize_ );
    U_ = Teuchos::null;
    V_ = Teuchos::null;
    U_ = rcp( new Epetra_MultiVector(ss->Map(),maxBasisSize_,false) );
    if (Vlocal_) {
      Epetra_LocalMap lclmap(maxNumCols_,0,ss->Comm());
      V_ = rcp( new Epetra_MultiVector(lclmap,maxBasisSize_,false) );
    }
    else {
      Epetra_Map gblmap(maxNumCols_,0,ss->Comm());
      V_ = rcp( new Epetra_MultiVector(gblmap,maxBasisSize_,false) );
    }
    B_ = rcp( new Epetra_SerialDenseMatrix(maxBasisSize_,maxBasisSize_) );
    resNorms_.reserve(maxBasisSize_);

    // clear counters
    numProc_ = 0;
    curNumPasses_ = 0;

    // we are now initialized, albeit with null rank
    isInitialized_ = true;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  info("Desired number of eigenvalues: %d.", NUMBER_OF_EIGENVALUES);

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform initial mesh refinements (optional).
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Enter boundary markers. 
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_ALL);

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_zero(BDY_ALL);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof: %d.", ndof);

  // Initialize the weak formulation for the left hand side, i.e., H.
  WeakForm wf_left, wf_right;
  wf_left.add_matrix_form(callback(bilinear_form_left));
  wf_right.add_matrix_form(callback(bilinear_form_right));

  // Initialize matrices.
  RCP<SparseMatrix> matrix_left = rcp(new CSCMatrix());
  RCP<SparseMatrix> matrix_right = rcp(new CSCMatrix());

  // Assemble the matrices.
  bool is_linear = true;
  DiscreteProblem dp_left(&wf_left, &space, is_linear);
  dp_left.assemble(matrix_left.get());
  DiscreteProblem dp_right(&wf_right, &space, is_linear);
  dp_right.assemble(matrix_right.get());

  EigenSolver es(matrix_left, matrix_right);
  info("Calling Pysparse...");
  es.solve(NUMBER_OF_EIGENVALUES, TARGET_VALUE, TOL, MAX_ITER);
  info("Pysparse finished.");
  es.print_eigenvalues();

  // Initializing solution vector, solution and ScalarView.
  double* coeff_vec;
  Solution sln;
  ScalarView view("Solution", new WinGeom(0, 0, 440, 350));

  // Reading solution vectors and visualizing.
  double* eigenval = new double[NUMBER_OF_EIGENVALUES];
  int neig = es.get_n_eigs();
  if (neig != NUMBER_OF_EIGENVALUES) error("Mismatched number of eigenvectors in the eigensolver output file.");  
  for (int ieig = 0; ieig < neig; ieig++) {
    eigenval[ieig] = es.get_eigenvalue(ieig);
    int n;
    es.get_eigenvector(ieig, &coeff_vec, &n);
    // Convert coefficient vector into a Solution.
    Solution::vector_to_solution(coeff_vec, &space, &sln);

    // Visualize the solution.
    char title[100];
    sprintf(title, "Solution %d, val = %g", ieig, eigenval[ieig]);
    view.set_title(title);
    view.show(&sln);

    // Wait for keypress.
    View::wait(HERMES_WAIT_KEYPRESS);
  }

  return 0; 
};

  //
  // Constructor
  //
  MyOp(const GlobalOrdinal n, const RCP< const Teuchos::Comm<int> > comm)
  {
    //
    // Construct a map for our block row distribution
    //
    opMap_ = rcp( new Map(n, 0, comm) );

    //
    // Get the rank of this process and the number of processes
    // We're going to have to do something special with the first and last processes
    //
    myRank_ = comm->getRank();
    numProcs_ = comm->getSize();

    //
    // Get the local number of rows
    //
    LocalOrdinal nlocal = opMap_->getNodeNumElements();

    //
    // Define the distribution that you need for the matvec.
    // When you define this for your own operator, it is helpful to draw pictures
    // on a sheet of paper to keep track of who needs to receive which elements.
    // Here, each process needs to receive one element from each of its neighbors.
    //

    // All processes but the first will receive one entry from the
    // previous process.
    if (myRank_ > 0) {
      nlocal++;
    }
    // All processes but the last will receive one entry from the next
    // process.
    if (myRank_ < numProcs_-1) {
      nlocal++;
    }
    // Construct a list of columns where this process has nonzero elements.
    // For our tridiagonal matrix, this is firstRowItOwns-1:lastRowItOwns+1.
    std::vector<GlobalOrdinal> indices;
    indices.reserve (nlocal);
    // The first process is a special case...
    if (myRank_ > 0) {
      indices.push_back (opMap_->getMinGlobalIndex () - 1);
    }
    for (GlobalOrdinal i = opMap_->getMinGlobalIndex ();
         i <= opMap_->getMaxGlobalIndex (); ++i) {
      indices.push_back (i);
    }
    // So is the last process...
    if (myRank_ < numProcs_-1) {
      indices.push_back (opMap_->getMaxGlobalIndex () + 1);
    }

    // Wrap our vector in an array view, which is like a raw pointer.
    Teuchos::ArrayView<const GlobalOrdinal> elementList (indices);

    // Make a Map for handling the redistribution.  There will be some
    // redundancies (i.e., some of the entries will be owned by
    // multiple MPI processes).
    GlobalOrdinal numGlobalElements = n + 2*(numProcs_-1);
    redistMap_ = rcp (new Map (numGlobalElements, elementList, 0, comm));

    // Make an Import object that describes how data will be
    // redistributed.  It takes a Map describing who owns what
    // originally, and a Map that describes who you WANT to own what.
    importer_= rcp (new Import (opMap_, redistMap_));
  }

bool run( const Teuchos::RCP<const Teuchos::Comm<int> > & comm ,
          const CMD & cmd)
{
  typedef typename Kokkos::Compat::KokkosDeviceWrapperNode<Device> NodeType;
  bool success = true;
  try {

  const int comm_rank = comm->getRank();

  // Create Tpetra Node -- do this first as it initializes host/device
  Teuchos::RCP<NodeType> node = createKokkosNode<NodeType>( cmd , *comm );

  // Set up stochastic discretization
  using Teuchos::Array;
  using Teuchos::RCP;
  using Teuchos::rcp;
  typedef Stokhos::OneDOrthogPolyBasis<int,double> one_d_basis;
  typedef Stokhos::LegendreBasis<int,double> legendre_basis;
  typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > order_type;
  typedef Stokhos::TotalOrderBasis<int,double,order_type> product_basis;
  typedef Stokhos::Sparse3Tensor<int,double> Cijk;
  const int dim = cmd.CMD_USE_UQ_DIM;
  const int order = cmd.CMD_USE_UQ_ORDER ;
  Array< RCP<const one_d_basis> > bases(dim);
  for (int i=0; i<dim; i++)
    bases[i] = rcp(new legendre_basis(order, true));
  RCP<const product_basis> basis = rcp(new product_basis(bases));
  RCP<Cijk> cijk = basis->computeTripleProductTensor();

  typedef Stokhos::DynamicStorage<int,double,Device> Storage;
  typedef Sacado::UQ::PCE<Storage> Scalar;
  typename Scalar::cijk_type kokkos_cijk =
    Stokhos::create_product_tensor<Device>(*basis, *cijk);
  Kokkos::setGlobalCijkTensor(kokkos_cijk);

  // typedef Stokhos::TensorProductQuadrature<int,double> quadrature;
  // RCP<const quadrature> quad     = rcp(new quadrature(basis));
  // const int num_quad_points                 = quad->size();
  // const Array<double>& quad_weights         = quad->getQuadWeights();
  // const Array< Array<double> >& quad_points = quad->getQuadPoints();
  // const Array< Array<double> >& quad_values = quad->getBasisAtQuadPoints();

  // Print output headers
  const std::vector< size_t > widths =
    print_headers( std::cout , cmd , comm_rank );

  using Kokkos::Example::FENL::TrivialManufacturedSolution;
  using Kokkos::Example::FENL::ElementComputationKLCoefficient;
  using Kokkos::Example::BoxElemPart;
  using Kokkos::Example::FENL::fenl;
  using Kokkos::Example::FENL::Perf;

  const double bc_lower_value = 1 ;
  const double bc_upper_value = 2 ;
  const TrivialManufacturedSolution manufactured_solution;

  int nelem[3] = { cmd.CMD_USE_FIXTURE_X  ,
                   cmd.CMD_USE_FIXTURE_Y  ,
                   cmd.CMD_USE_FIXTURE_Z  };

  // Create KL diffusion coefficient
  const double kl_mean = cmd.CMD_USE_MEAN;
  const double kl_variance = cmd.CMD_USE_VAR;
  const double kl_correlation = cmd.CMD_USE_COR;
  typedef ElementComputationKLCoefficient< Scalar, double, Device > KL;
  KL diffusion_coefficient( kl_mean, kl_variance, kl_correlation, dim );
  typedef typename KL::RandomVariableView RV;
  typedef typename RV::HostMirror HRV;
  RV rv = diffusion_coefficient.getRandomVariables();
  HRV hrv = Kokkos::create_mirror_view(rv);

  // Set random variables
  // ith random variable \xi_i = \psi_I(\xi) / \psi_I(1.0)
  // where I is determined by the basis ordering (since the component basis
  // functions have unit two-norm, \psi_I(1.0) might not be 1.0).  We compute
  // this by finding the index of the multivariate term that is first order in
  // the ith slot, all other orders 0
  Teuchos::Array<double> point(dim, 1.0);
  Teuchos::Array<double> basis_vals(basis->size());
  basis->evaluateBases(point, basis_vals);
  for (int i=0; i<dim; ++i) {
    Stokhos::MultiIndex<int> term(dim, 0);
    term[i] = 1;
    int index = basis->index(term);
    hrv(i).fastAccessCoeff(index) = 1.0 / basis_vals[index];
  }
  Kokkos::deep_copy( rv, hrv );

  // Compute stochastic response using stochastic Galerkin method
  Scalar response = 0;
  Perf perf;
  if ( cmd.CMD_USE_FIXTURE_QUADRATIC  )
    perf = fenl< Scalar , Device , BoxElemPart::ElemQuadratic >
      ( comm , node , cmd.CMD_PRINT , cmd.CMD_USE_TRIALS ,
        cmd.CMD_USE_ATOMIC , cmd.CMD_USE_BELOS , cmd.CMD_USE_MUELU ,
        cmd.CMD_USE_MEANBASED ,
        nelem , diffusion_coefficient , manufactured_solution ,
        bc_lower_value , bc_upper_value ,
        false , response);
  else
//.........这里部分代码省略.........