C++ DefaultLinearSolverBuilder::createLinearSolveStrategy方法代码示例

  virtual
  Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<double>>
  get_W_factory() const
  {
    Stratimikos::DefaultLinearSolverBuilder builder;

    const std::map<std::string, boost::any> linear_solver_params = {
      {"package", std::string("Belos")},
      {"method", std::string("Pseudo Block CG")},
      {"parameters", dict{
        {"Output Frequency", 1},
        {"Verbosity", 0}
      }}
      };
    auto p = Teuchos::rcp(new Teuchos::ParameterList());
    mikado::std_map_to_teuchos_list(
        mikado::convert_to_belos_parameters(linear_solver_params),
        *p
        );
    builder.setParameterList(p);

    auto lowsFactory = builder.createLinearSolveStrategy("");

    lowsFactory->setVerbLevel(Teuchos::VERB_LOW);

    return lowsFactory;
  }

TEUCHOS_UNIT_TEST(belos_gcrodr, multiple_solves)
{
   // build global (or serial communicator)
   #ifdef HAVE_MPI
      Epetra_MpiComm Comm(MPI_COMM_WORLD);
   #else
      Epetra_SerialComm Comm;
   #endif

   // build and allocate linear system
   Teuchos::RCP<Epetra_CrsMatrix> mat = buildMatrix(100,Comm);
   Teuchos::RCP<Epetra_Vector> x0 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
   Teuchos::RCP<Epetra_Vector> x1 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
   Teuchos::RCP<Epetra_Vector> b = rcp(new Epetra_Vector(mat->OperatorRangeMap()));

   x0->Random();
   x1->Random();
   b->PutScalar(0.0);

   // sanity check
   // EpetraExt::RowMatrixToMatrixMarketFile("mat_output.mm",*mat);

   // build Thyra wrappers
   RCP<const Thyra::LinearOpBase<double> >
      tA = Thyra::epetraLinearOp( mat );
   RCP<Thyra::VectorBase<double> >
      tx0 = Thyra::create_Vector( x0, tA->domain() );
   RCP<Thyra::VectorBase<double> >
      tx1 = Thyra::create_Vector( x1, tA->domain() );
   RCP<const Thyra::VectorBase<double> >
      tb = Thyra::create_Vector( b, tA->range() );

   // now comes Stratimikos
   RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("BelosGCRODRTest.xml");
   Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
   linearSolverBuilder.setParameterList(paramList);
 
   // Create a linear solver factory given information read from the
   // parameter list.
   RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
         linearSolverBuilder.createLinearSolveStrategy("");

   // Create a linear solver based on the forward operator A
   RCP<Thyra::LinearOpWithSolveBase<double> > lows =
         Thyra::linearOpWithSolve(*lowsFactory, tA);

   // Solve the linear system 
   Thyra::SolveStatus<double> status; 
   status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx0.ptr());
   status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx1.ptr());
}

  virtual
  Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<double>>
  get_W_factory() const
  {
    Stratimikos::DefaultLinearSolverBuilder builder;
    auto p = Teuchos::rcp(new Teuchos::ParameterList());
    mikado::std_map_to_teuchos_list(
        mikado::convert_to_belos_parameters(this->linear_solver_params_),
        *p
        );
    builder.setParameterList(p);

    auto lowsFactory = builder.createLinearSolveStrategy("");

    lowsFactory->setVerbLevel(Teuchos::VERB_LOW);

    return lowsFactory;
  }

//.........这里部分代码省略.........
    CommandLineProcessor  clp(false); // Don't throw exceptions

    lowsfCreator.setupCLP(&clp);
    solver.setupCLP(&clp);

    clp.setOption( "xt0", &xt0 );
    clp.setOption( "xt1", &xt1 );
    clp.setOption( "pt0", &pt0 );
    clp.setOption( "pt1", &pt1 );
    clp.setOption( "d", &d );
    clp.setOption( "x00", &x00 );
    clp.setOption( "x01", &x01 );
    clp.setOption( "p00", &p00 );
    clp.setOption( "p01", &p01 );
    clp.setOption( "pL0", &pL0 );
    clp.setOption( "pL1", &pL1 );
    clp.setOption( "pU0", &pU0 );
    clp.setOption( "pU1", &pU1 );
    clp.setOption( "xL0", &xL0 );
    clp.setOption( "xL1", &xL1 );
    clp.setOption( "xU0", &xU0 );
    clp.setOption( "xU1", &xU1 );
    clp.setOption( "support-derivs", "no-support-derivs", &supportDerivs );
    clp.setOption("extra-xml-file",&extraXmlFile,"File with extra XML text that will modify the initial XML read in");
 
    std::string line("");
    if(extraXmlFile.length()) {
      std::ifstream myfile(extraXmlFile.c_str());
      if (myfile.is_open())
      {
        getline (myfile,line);
        solver.extraParamsXmlStringOption(line);
        std::cout << line << "\n";
        myfile.close();
      }
    }

    CommandLineProcessor::EParseCommandLineReturn
      parse_return = clp.parse(argc,argv,&std::cerr);

    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
      return parse_return;

    lowsfCreator.readParameters(out.get());
    solver.readParameters(out.get());

    //
    // Create the NLP
    //
    
    // Create the EpetraExt::ModelEvaluator object

    Teuchos::RCP<EpetraModelEval4DOpt>
      epetraModel = Teuchos::rcp(new EpetraModelEval4DOpt(xt0,xt1,pt0,pt1,d,x00,x01,p00,p01));
    epetraModel->setSupportDerivs(supportDerivs);
    epetraModel->set_p_bounds(pL0,pL1,pU0,pU1);
    epetraModel->set_x_bounds(xL0,xL1,xU0,xU1);

    // Create the Thyra::EpetraModelEvaluator object

    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
      lowsFactory = lowsfCreator.createLinearSolveStrategy("");

    Teuchos::RCP<Thyra::EpetraModelEvaluator>
      epetraThyraModel = Teuchos::rcp(new Thyra::EpetraModelEvaluator());
    
    epetraThyraModel->initialize(epetraModel,lowsFactory);
    
    //
    // Solve the NLP
    //
    
    // Set the model
    solver.setModel(epetraThyraModel);

    // Read the initial guess if one exists
    solver.readInitialGuess(out.get());

    // Solve the NLP
    const MoochoSolver::ESolutionStatus	solution_status = solver.solve();

    // Write the final solution if requested
    solver.writeFinalSolution(out.get());

    // Write the parameters that where read
    lowsfCreator.writeParamsFile(*lowsFactory);
    solver.writeParamsFile();
    
    //
    // Return the solution status (0 if sucessfull)
    //

    return solution_status;

  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,dummySuccess)

  return MoochoSolver::SOLVE_RETURN_EXCEPTION;

}

  void solvetriadmatrixwithtrilinos(int& nnz, int& order, int* row, 
              int* col, double* val, double* rhs, double* solution) {
#else
  void solvetriadmatrixwithtrilinos_(int& nnz, int& order, int* row, 
              int* col, double* val, double* rhs, double* solution) {
#endif

    try{
    
#ifdef _MPI
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm Comm;
#endif
    
    int i, j, ierr;
    int MyPID = Comm.MyPID();
    bool verbose = (MyPID == 0);
    Epetra_Map RowMap(order, 0, Comm);
    int NumMyElements = RowMap.NumMyElements();
    int *MyGlobalElements = new int[NumMyElements];
    RowMap.MyGlobalElements(&MyGlobalElements[0]);

#ifdef _MPI
    int nPEs;
    MPI_Comm_size(MPI_COMM_WORLD, &nPEs);
#endif

    int anEst = nnz / order + 1;
    Epetra_CrsMatrix A(Copy, RowMap, anEst);
    
    for (j=0; j<nnz; ++j) {
      if (RowMap.MyGID(row[j]) ) {
	ierr = A.InsertGlobalValues(row[j], 1, &(val[j]), &(col[j]) );
	assert(ierr >= 0);
      }
    }
    
    ierr = A.FillComplete();
    assert(ierr == 0);
    
    //-------------------------------------------------------------------------
    // RN_20091221: Taking care of the rhs
    //-------------------------------------------------------------------------
    Epetra_Vector b(RowMap);

    // Inserting values into the rhs
    double *MyGlobalValues = new double[NumMyElements];
    for (j=0; j<NumMyElements; ++j) {
      MyGlobalValues[j] = rhs[MyGlobalElements[j] ];
    }
    ierr = b.ReplaceGlobalValues(NumMyElements, &MyGlobalValues[0],
				 &MyGlobalElements[0]);

    //-------------------------------------------------------------------------
    // RN_20091221: Taking care of the solution
    //-------------------------------------------------------------------------
    Epetra_Vector x(RowMap);

    Teuchos::ParameterList paramList;

    Teuchos::RCP<Teuchos::ParameterList>
      paramList1 = Teuchos::rcp(&paramList, false);
    Teuchos::updateParametersFromXmlFile("./strat1.xml", paramList1.get() );

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

    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;

    Teuchos::RCP<Epetra_CrsMatrix> epetraOper = Teuchos::rcp(&A, false);
    Teuchos::RCP<Epetra_Vector> epetraRhs = Teuchos::rcp(&b, false);
    Teuchos::RCP<Epetra_Vector> epetraSol = Teuchos::rcp(&x, false);

    Teuchos::RCP<const Thyra::LinearOpBase<double> >
      thyraOper = Thyra::epetraLinearOp(epetraOper);
    Teuchos::RCP<Thyra::VectorBase<double> >
      thyraRhs = Thyra::create_Vector(epetraRhs, thyraOper->range() );
    Teuchos::RCP<Thyra::VectorBase<double> >
      thyraSol = Thyra::create_Vector(epetraSol, thyraOper->domain() );

    linearSolverBuilder.setParameterList(Teuchos::rcp(&paramList, false) );

    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
      lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");

    lowsFactory->setOStream(out);
    lowsFactory->setVerbLevel(Teuchos::VERB_LOW);

    Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> >
      lows = Thyra::linearOpWithSolve(*lowsFactory, thyraOper);

    Thyra::SolveStatus<double>
      status = Thyra::solve(*lows, Thyra::NOTRANS, *thyraRhs, &*thyraSol);

    thyraSol = Teuchos::null;

    // For debugging =)
    //    cout << "A: " << A << endl;
    //    cout << "b: " << b << endl;
//.........这里部分代码省略.........

  //================================================================
  //================================================================
  // RN_20091215: This needs to be called only once per time step
  // in the beginning to set up the problem.
  //================================================================
  void FC_FUNC(inittrilinos,INITTRILINOS) (int& bandwidth, int& mySize,
	       int* myIndicies, double* myX, double* myY, double* myZ,
	       int* mpi_comm_f) {
// mpi_comm_f: CISM's fortran mpi communicator

#ifdef GLIMMER_MPI
    // Make sure the MPI_Init in Fortran is recognized by C++.
    // We used to call an extra MPI_Init if (!flag), but the behavior of doing so is uncertain,
    // especially if CISM's MPI communicator is a subset of MPI_COMM_WORLD (as can be the case in CESM).
    // Thus, for now, we die with an error message if C++ perceives MPI to be uninitialized.
    // If this causes problems (e.g., if certain MPI implementations seem not to recognize 
    // that MPI has already been initialized), then we will revisit how to handle this.
       int flag;
       MPI_Initialized(&flag);
       if (!flag) {
	 cout << "ERROR in inittrilinos: MPI not initialized according to C++ code" << endl;
	 exit(1);
       }
    MPI_Comm mpi_comm_c = MPI_Comm_f2c(*mpi_comm_f);
    Epetra_MpiComm comm(mpi_comm_c);
    Teuchos::MpiComm<int> tcomm(Teuchos::opaqueWrapper(mpi_comm_c));
#else
    Epetra_SerialComm comm;
    Teuchos::SerialComm<int> tcomm;
#endif

    Teuchos::RCP<const Epetra_Map> rowMap = 
      Teuchos::rcp(new Epetra_Map(-1,mySize,myIndicies,1,comm) );

    TEUCHOS_TEST_FOR_EXCEPTION(!rowMap->UniqueGIDs(), std::logic_error,
       "Error: inittrilinos, myIndices array needs to have Unique entries" 
        << " across all processor.");

    // Diagnostic output for partitioning
    int minSize, maxSize;
    comm.MinAll(&mySize, &minSize, 1);
    comm.MaxAll(&mySize, &maxSize, 1);
    if (comm.MyPID()==0) 
      cout << "\nPartition Info in init_trilinos: Total nodes = " << rowMap->NumGlobalElements()
           << "  Max = " << maxSize << "  Min = " << minSize 
           << "  Ave = " << rowMap->NumGlobalElements() / comm.NumProc() << endl;

    soln = Teuchos::rcp(new Epetra_Vector(*rowMap));

    // Read parameter list once
    try { 
       pl = Teuchos::rcp(new Teuchos::ParameterList("Trilinos Options"));
       Teuchos::updateParametersFromXmlFileAndBroadcast("trilinosOptions.xml", pl.ptr(), tcomm);

       Teuchos::ParameterList validPL("Valid List");;
       validPL.sublist("Stratimikos"); validPL.sublist("Piro");
       pl->validateParameters(validPL, 0);
    }
    catch (std::exception& e) {
      cout << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n" 
           << e.what() << "\nExiting: Invalid trilinosOptions.xml file."
           << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << endl;
      exit(1);
    }
    catch (...) {
      cout << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n" 
           << "\nExiting: Invalid trilinosOptions.xml file."
           << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << endl;
      exit(1);
    }

    try { 
      // Set the coordinate position of the nodes for ML for repartitioning (important for #procs > 100s)
      if (pl->sublist("Stratimikos").isParameter("Preconditioner Type")) {
         if ("ML" == pl->sublist("Stratimikos").get<string>("Preconditioner Type")) {
           Teuchos::ParameterList& mlList =
              pl->sublist("Stratimikos").sublist("Preconditioner Types").sublist("ML").sublist("ML Settings");
           mlList.set("x-coordinates",myX);
           mlList.set("y-coordinates",myY);
           mlList.set("z-coordinates",myZ);
           mlList.set("PDE equations", 1);
         }
      }

      out = Teuchos::VerboseObjectBase::getDefaultOStream();

      // Reset counters every time step: can remove these lines to have averages over entire run
      linearSolveIters_total = 0;
      linearSolveCount=0;
      linearSolveSuccessCount = 0;

      // Create an interface that holds a CrsMatrix instance and some useful methods.
      interface = Teuchos::rcp(new TrilinosMatrix_Interface(rowMap, bandwidth, comm));

      Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
      linearSolverBuilder.setParameterList(Teuchos::sublist(pl, "Stratimikos"));
      lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
      lowsFactory->setOStream(out);
      lowsFactory->setVerbLevel(Teuchos::VERB_LOW);

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

  TEUCHOS_UNIT_TEST(explicit_model_evaluator, basic)
  {
    using Teuchos::RCP;

    PHX::KokkosDeviceSession session;

    bool parameter_on = true;
    Teuchos::RCP<panzer::FieldManagerBuilder> fmb;  
    Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > rLibrary; 
    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof;
    Teuchos::RCP<panzer::GlobalData> gd;
  
    buildAssemblyPieces(parameter_on,fmb,rLibrary,gd,lof);

    // Test a transient me
    {
      typedef Thyra::ModelEvaluatorBase MEB;
      typedef Thyra::ModelEvaluatorBase::InArgs<double> InArgs;
      typedef Thyra::ModelEvaluatorBase::OutArgs<double> OutArgs;
      typedef Thyra::VectorBase<double> VectorType;
      typedef Thyra::LinearOpBase<double> OperatorType;
      typedef panzer::ModelEvaluator<double> PME;
      typedef panzer::ExplicitModelEvaluator<double> ExpPME;

      std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names;
      std::vector<Teuchos::RCP<Teuchos::Array<double> > > p_values;
      bool build_transient_support = true;

      Stratimikos::DefaultLinearSolverBuilder builder;
      Teuchos::RCP<Teuchos::ParameterList> validList = Teuchos::rcp(new Teuchos::ParameterList(*builder.getValidParameters()));
      builder.setParameterList(validList);
      RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = builder.createLinearSolveStrategy("Amesos");
    
      RCP<PME> me = Teuchos::rcp(new PME(fmb,rLibrary,lof,p_names,p_values,lowsFactory,gd,build_transient_support,0.0));
      RCP<ExpPME> exp_me = Teuchos::rcp(new ExpPME(me,true,false)); // constant mass, use lumped

      RCP<VectorType> exp_f, f;

      // explicit evaluation
      {
        // set the nominal values
        InArgs nom_vals = exp_me->getNominalValues();
        TEST_ASSERT(nom_vals.supports(MEB::IN_ARG_x));
        TEST_ASSERT(nom_vals.supports(MEB::IN_ARG_x_dot)); // this is supported for stabilization purposes
        TEST_ASSERT(nom_vals.supports(MEB::IN_ARG_alpha)); // alpha and beta support needed for outputting responses
        TEST_ASSERT(nom_vals.supports(MEB::IN_ARG_beta));

        // create in args
        InArgs in_args = exp_me->createInArgs();
        TEST_ASSERT(in_args.supports(MEB::IN_ARG_x));
        TEST_ASSERT(in_args.supports(MEB::IN_ARG_x_dot)); // this is supported for stabilization purposes
        TEST_ASSERT(in_args.supports(MEB::IN_ARG_alpha)); // alpha and beta support needed for outputting responses
        TEST_ASSERT(in_args.supports(MEB::IN_ARG_beta));
        InArgs nomValues = exp_me->getNominalValues();
        RCP<VectorType> x = Thyra::createMember(*exp_me->get_x_space());
        RCP<VectorType> x_dot = Thyra::createMember(*exp_me->get_x_space());
        Thyra::assign(x_dot.ptr(),0.0);
        Thyra::assign(x.ptr(),5.0);
        in_args.set_x(x);
        in_args.set_x_dot(x_dot);
        
        // create out args
        OutArgs out_args = exp_me->createOutArgs();
        TEST_ASSERT(out_args.supports(MEB::OUT_ARG_f));
        TEST_ASSERT(!out_args.supports(MEB::OUT_ARG_W_op));
        TEST_ASSERT(!out_args.supports(MEB::OUT_ARG_W));
        exp_f = Thyra::createMember(*exp_me->get_f_space());
        out_args.set_f(exp_f);

        exp_me->evalModel(in_args, out_args);
      }

      // implicit evaluation
      RCP<OperatorType> mass = me->create_W_op();
      {
        // create in args
        InArgs in_args = me->createInArgs();
        InArgs nomValues = me->getNominalValues();
        RCP<VectorType> x = Thyra::createMember(*me->get_x_space());
        RCP<VectorType> x_dot = Thyra::createMember(*me->get_x_space());
        Thyra::assign(x_dot.ptr(),0.0);
        Thyra::assign(x.ptr(),5.0);
        in_args.set_x(x);
        in_args.set_x_dot(x_dot);
        
        // create out args
        OutArgs out_args = me->createOutArgs();
        f = Thyra::createMember(*me->get_f_space());
        out_args.set_f(f);

        me->evalModel(in_args, out_args);

        in_args.set_x(x);
        in_args.set_x_dot(x_dot);
        in_args.set_alpha(1.0);
        in_args.set_beta(0.0);

        out_args.set_f(Teuchos::null);
        out_args.set_W_op(mass);

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

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

   // setup field manager build
   /////////////////////////////////////////////////////////////
   out << "SETUP FMB" << std::endl;
 
    // Add in the application specific closure model factory
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory; 
    user_app::MyModelFactory_TemplateBuilder cm_builder;
    cm_factory.buildObjects(cm_builder);

    Teuchos::ParameterList closure_models("Closure Models");
    closure_models.sublist("solid").sublist("SOURCE_TEMPERATURE").set<double>("Value",1.0);
    closure_models.sublist("solid").sublist("DENSITY").set<double>("Value",1.0);
    closure_models.sublist("solid").sublist("HEAT_CAPACITY").set<double>("Value",1.0);
    closure_models.sublist("ion solid").sublist("SOURCE_ION_TEMPERATURE").set<double>("Value",1.0);
    closure_models.sublist("ion solid").sublist("ION_DENSITY").set<double>("Value",1.0);
    closure_models.sublist("ion solid").sublist("ION_HEAT_CAPACITY").set<double>("Value",1.0);

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

    fmb->setWorksetContainer(wkstContainer);
    fmb->setupVolumeFieldManagers(physicsBlocks,cm_factory,closure_models,*linObjFactory,user_data);
    fmb->setupBCFieldManagers(bcs,physicsBlocks,*eqset_factory,cm_factory,bc_factory,closure_models,*linObjFactory,user_data);

   // setup assembly engine
   /////////////////////////////////////////////////////////////

   // build assembly engine
   panzer::AssemblyEngine_TemplateManager<panzer::Traits> ae_tm;
   panzer::AssemblyEngine_TemplateBuilder builder(fmb,linObjFactory);
   ae_tm.buildObjects(builder);

   // setup linear algebra and solve 
   /////////////////////////////////////////////////////////////

   // build ghosted variables
   out << "BUILD LA" << std::endl;
   RCP<panzer::EpetraLinearObjContainer> ghostCont 
         = Teuchos::rcp_dynamic_cast<panzer::EpetraLinearObjContainer>(linObjFactory->buildGhostedLinearObjContainer());
   RCP<panzer::EpetraLinearObjContainer> container 
         = Teuchos::rcp_dynamic_cast<panzer::EpetraLinearObjContainer>(linObjFactory->buildLinearObjContainer());
   eLinObjFactory->initializeContainer(panzer::EpetraLinearObjContainer::X |
                                       panzer::EpetraLinearObjContainer::DxDt |
                                       panzer::EpetraLinearObjContainer::F |
                                       panzer::EpetraLinearObjContainer::Mat,*container);
   eLinObjFactory->initializeGhostedContainer(panzer::EpetraLinearObjContainer::X |
                                              panzer::EpetraLinearObjContainer::DxDt |
                                              panzer::EpetraLinearObjContainer::F |
                                              panzer::EpetraLinearObjContainer::Mat,*ghostCont);

   panzer::AssemblyEngineInArgs input(ghostCont,container);
   input.alpha = 0;
   input.beta = 1;

   // evaluate physics
   out << "EVALUTE" << std::endl;
   ae_tm.getAsObject<panzer::Traits::Residual>()->evaluate(input);
   ae_tm.getAsObject<panzer::Traits::Jacobian>()->evaluate(input);

   out << "RAN SUCCESSFULLY!" << std::endl;

   out << "SOLVE" << std::endl;

   // notice that this should be called by the assembly driver!
   // linObjFactory->ghostToGlobalContainer(*ghostCont,*container);

   Teuchos::RCP<const Thyra::LinearOpBase<double> > th_A = Thyra::epetraLinearOp(container->get_A());
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > range  = th_A->range();
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domain = th_A->domain();

   Teuchos::RCP<Thyra::VectorBase<double> > th_x = Thyra::create_Vector(container->get_x(),domain);
   Teuchos::RCP<Thyra::VectorBase<double> > th_f = Thyra::create_Vector(container->get_f(),range);

   // solve with amesos
   Stratimikos::DefaultLinearSolverBuilder solverBuilder;
   Teuchos::RCP<Teuchos::ParameterList> validList = Teuchos::rcp(new Teuchos::ParameterList(*solverBuilder.getValidParameters()));
   solverBuilder.setParameterList(validList);
   
   RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = solverBuilder.createLinearSolveStrategy("Amesos");
   RCP<Thyra::LinearOpWithSolveBase<double> > lows = Thyra::linearOpWithSolve(*lowsFactory, th_A.getConst());
   Thyra::solve<double>(*lows, Thyra::NOTRANS, *th_f, th_x.ptr());

   if(false) {
      EpetraExt::RowMatrixToMatrixMarketFile("a_op.mm",*container->get_A());
      EpetraExt::VectorToMatrixMarketFile("x_vec.mm",*container->get_x());
      EpetraExt::VectorToMatrixMarketFile("b_vec.mm",*container->get_f());
   }

   out << "WRITE" << std::endl;

   // redistribute solution vector
   linObjFactory->globalToGhostContainer(*container,*ghostCont,panzer::EpetraLinearObjContainer::X | panzer::EpetraLinearObjContainer::DxDt); 

   panzer_stk::write_solution_data(*dofManager,*mesh,*ghostCont->get_x());
   mesh->writeToExodus("output.exo");

   PHX::FinalizeKokkosDevice();

   return 0;
}

//.........这里部分代码省略.........
          ,"  "
          )
        );
    journalOut->copyAllOutputOptions(*out);

    *out
      << "\n***"
      << "\n*** NLPThyraEpetraAdvDiffReactOptMain, Global numProcs = "<<numProcs
      << "\n***\n";

#ifdef HAVE_MPI
    MPI_Comm mpiComm = MPI_COMM_WORLD;
#endif

    Teuchos::RCP<Epetra_Comm> comm = Teuchos::null;
#ifdef HAVE_MPI
    comm = Teuchos::rcp(new Epetra_MpiComm(mpiComm));
#else
    comm = Teuchos::rcp(new Epetra_SerialComm());
#endif
    
    //
    // Create the Thyra::ModelEvaluator object
    //
    
    *out << "\nCreate the GLpApp::AdvDiffReactOptModel wrapper object ...\n";
    
    Teuchos::RCP<GLpApp::AdvDiffReactOptModel>
      epetraModel = advDiffReacModelCreator.createModel(comm);
    epetraModel->setOStream(journalOut);

    *out << "\nCreate the Thyra::LinearOpWithSolveFactory object ...\n";

    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
      lowsFactory = lowsfCreator.createLinearSolveStrategy("");
    // ToDo: Set the output stream before calling above!
    ///lowsFactory = lowsfCreator.createLOWSF(OSTab(journalOut).get());
    
    *out << "\nCreate the Thyra::EpetraModelEvaluator wrapper object ...\n";
    
    Teuchos::RCP<Thyra::EpetraModelEvaluator>
      epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator()); // Sets default options!
    epetraThyraModel->setOStream(journalOut);
    epetraThyraModel->initialize(epetraModel,lowsFactory);

    *out
      << "\nnx = " << epetraThyraModel->get_x_space()->dim()
      << "\nnp = " << epetraThyraModel->get_p_space(0)->dim() << "\n";

    if(matchingVecFile != "") {
      *out << "\nReading the matching vector \'q\' from the file(s) with base name \""<<matchingVecFile<<"\" ...\n";
      Thyra::DefaultSpmdMultiVectorFileIO<Scalar> fileIO;
      epetraModel->set_q(
        Thyra::get_Epetra_Vector(
          *epetraModel->get_x_map()
          ,readVectorFromFile(fileIO,matchingVecFile,*epetraThyraModel->get_x_space())
          )
        );
    }

    //
    // Solve the NLP
    //

    // Set the journal file
    solver.getSolver().set_journal_out(journalOut);
    
    // Set the model
    solver.setModel(epetraThyraModel);

    // Read the initial guess if one exists
    solver.readInitialGuess(out.get());

    // Solve the NLP
    const MoochoSolver::ESolutionStatus	solution_status = solver.solve();

    // Write the solution to file
    solver.writeFinalSolution(out.get());
    
    // Write the final parameters to file
    lowsfCreator.writeParamsFile(*lowsFactory);
    solver.writeParamsFile();
    
    //
    // Return the solution status (0 if successful)
    //

    if(solution_status == MoochoSolver::SOLVE_RETURN_SOLVED)
      *out << "\nEnd Result: TEST PASSED\n";
    else
      *out << "\nEnd Result: TEST FAILED\n";
    
    return solution_status;

  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true, *out, dummySuccess)

  return MoochoSolver::SOLVE_RETURN_EXCEPTION;

}