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

void TrilinosCrsMatrix::solve(LSS::Vector& solution, LSS::Vector& rhs)
{
  cf3_assert(m_is_created);
  cf3_assert(solution.is_created());
  cf3_assert(rhs.is_created());

  LSS::TrilinosVector& tsol = dynamic_cast<LSS::TrilinosVector&>(solution);
  LSS::TrilinosVector& trhs = dynamic_cast<LSS::TrilinosVector&>(rhs);

  Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile(options().option("settings_file").value_str());

  // Build Thyra linear algebra objects
  Teuchos::RCP<const Thyra::LinearOpBase<double> > th_mat = Thyra::epetraLinearOp(m_mat);
  Teuchos::RCP<const Thyra::VectorBase<double> > th_rhs = Thyra::create_Vector(trhs.epetra_vector(),th_mat->range());
  Teuchos::RCP<Thyra::VectorBase<double> > th_sol = Thyra::create_Vector(tsol.epetra_vector(),th_mat->domain());

  // Build stratimikos solver
  /////////////////////////////////////////////////////////

  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;

  //Teko::addTekoToStratimikosBuilder(linearSolverBuilder);
  linearSolverBuilder.setParameterList(paramList);

  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > th_invA = Thyra::linearOpWithSolve(*lowsFactory, th_mat);

  Thyra::assign(th_sol.ptr(), 0.0);
  Thyra::SolveStatus<double> status = Thyra::solve<double>(*th_invA, Thyra::NOTRANS, *th_rhs, th_sol.ptr());
  CFinfo << "Thyra::solve finished with status " << status.message << CFendl;
}

  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(tStratimikosFactory, test_RelatedFunctions) 
{
  using Teuchos::RCP;
  using Teuchos::ParameterList;

  // build global (or serial communicator)
  #ifdef HAVE_MPI
     Epetra_MpiComm comm(MPI_COMM_WORLD);
  #else
     Epetra_SerialComm comm;
  #endif

  // build epetra operator
  RCP<Epetra_Operator> eA = buildStridedSystem(comm,5);
  RCP<Thyra::LinearOpBase<double> > tA = Thyra::nonconstEpetraLinearOp(eA);

  {
     // build stratimikos factory, adding Teko's version
     Stratimikos::DefaultLinearSolverBuilder stratFactory;
     
     Teko::addTekoToStratimikosBuilder(stratFactory);
     TEST_THROW(Teko::addTekoToStratimikosBuilder(stratFactory),std::logic_error);
     Teko::addTekoToStratimikosBuilder(stratFactory,"Teko-2");

     TEST_NOTHROW(stratFactory.getValidParameters()->sublist("Preconditioner Types").sublist("Teko"));
     TEST_NOTHROW(stratFactory.getValidParameters()->sublist("Preconditioner Types").sublist("Teko-2"));
  }

  {
     Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler);

     // build stratimikos factory, adding Teko's version
     Stratimikos::DefaultLinearSolverBuilder stratFactory;
     
     Teko::addTekoToStratimikosBuilder(stratFactory,rh);
     TEST_THROW(Teko::addTekoToStratimikosBuilder(stratFactory,rh),std::logic_error);
     Teko::addTekoToStratimikosBuilder(stratFactory,rh,"Teko-2");

     TEST_NOTHROW(stratFactory.getValidParameters()->sublist("Preconditioner Types").sublist("Teko"));
     TEST_NOTHROW(stratFactory.getValidParameters()->sublist("Preconditioner Types").sublist("Teko-2"));

     RCP<ParameterList> params = Teuchos::rcp(new ParameterList(*stratFactory.getValidParameters()));
     ParameterList & tekoList = params->sublist("Preconditioner Types").sublist("Teko");
     tekoList.set("Write Block Operator", false);
     tekoList.set("Test Block Operator", false);
     tekoList.set("Strided Blocking","1 1");
     tekoList.set("Inverse Type","BGS");
     ParameterList & ifl = tekoList.sublist("Inverse Factory Library");
     ifl.sublist("BGS").set("Type","Block Gauss-Seidel");
     ifl.sublist("BGS").set("Inverse Type","Amesos");
     stratFactory.setParameterList(params);

     RCP<Teko::StratimikosFactory> precFactory
           = Teuchos::rcp_dynamic_cast<Teko::StratimikosFactory>(stratFactory.createPreconditioningStrategy("Teko-2"));
     TEST_EQUALITY(precFactory->getRequestHandler(),rh);
  }
}

int main(int argc,char * argv[])
{
   // calls MPI_Init and MPI_Finalize
   Teuchos::GlobalMPISession mpiSession(&argc,&argv);

   // read in parameter list
   Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("strat_example.xml"); 

   // build global communicator
#ifdef HAVE_MPI
   Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
   Epetra_SerialComm Comm;
#endif

   // Read in the matrix, store pointer as an RCP
   Epetra_CrsMatrix * ptrA = 0;
   EpetraExt::MatrixMarketFileToCrsMatrix("../data/nsjac.mm",Comm,ptrA);
   RCP<Epetra_CrsMatrix> A = rcp(ptrA);

   // read in the RHS vector
   Epetra_Vector * ptrb = 0;
   EpetraExt::MatrixMarketFileToVector("../data/nsrhs_test.mm",A->OperatorRangeMap(),ptrb);
   RCP<const Epetra_Vector> b = rcp(ptrb);

   // allocate vectors
   RCP<Epetra_Vector> x = rcp(new Epetra_Vector(A->OperatorDomainMap()));
   x->PutScalar(0.0);

   // Build Thyra linear algebra objects
   RCP<const Thyra::LinearOpBase<double> > th_A = Thyra::epetraLinearOp(A);
   RCP<const Thyra::VectorBase<double> > th_b = Thyra::create_Vector(b,th_A->range());
   RCP<Thyra::VectorBase<double> > th_x = Thyra::create_Vector(x,th_A->domain());

   // Build stratimikos solver
   /////////////////////////////////////////////////////////

   Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;

   Teko::addTekoToStratimikosBuilder(linearSolverBuilder);
   linearSolverBuilder.setParameterList(paramList);

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

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

   Thyra::assign(th_x.ptr(), 0.0);
   Thyra::SolveStatus<double> status 
       = Thyra::solve<double>(*th_invA, Thyra::NOTRANS, *th_b, th_x.ptr());

   return 0;
}

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

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

        // Note: Use --help on the command line to see the above documentation

        CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
        if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;


        //
        *out << "\nA) Reading in the matrix ...\n";
        //

#ifdef HAVE_MPI
        Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
        Epetra_SerialComm comm;
#endif

        const LinearOpPtr A = readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
                                  matrixFile, comm, "A");
        *out << "\nA = " << describe(*A,verbLevel) << "\n";

        const RCP<ParameterList> paramList = getParametersFromXmlFile(paramListFile);
        if (showParams) {
            *out << "\nRead in parameter list:\n\n";
            paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
        }

        //
        *out << "\nB) Get the preconditioner as a forward solver\n";
        //

        const RCP<ParameterList> precParamList = sublist(paramList, "Preconditioner Solver");
        Stratimikos::DefaultLinearSolverBuilder precSolverBuilder;
        precSolverBuilder.setParameterList(precParamList);
        const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > precSolverStrategy
            = createLinearSolveStrategy(precSolverBuilder);
        //precSolverStrategy->setVerbLevel(verbLevel);

        const LinearOpPtr A_inv_prec = inverse<double>(*precSolverStrategy, A,
                                       Thyra::SUPPORT_SOLVE_FORWARD_ONLY,
                                       Teuchos::null, // Use internal solve criteria
                                       Thyra::IGNORE_SOLVE_FAILURE // Ignore solve failures since this is just a prec
                                                      );
        *out << "\nA_inv_prec = " << describe(*A_inv_prec, verbLevel) << "\n";

        if (testPrecIsLinearOp) {
            *out << "\nTest that the preconditioner A_inv_prec is indeed a linear operator.\n";
            Thyra::LinearOpTester<double> linearOpTester;
            linearOpTester.check_adjoint(false);
            const bool linearOpCheck = linearOpTester.check(*A_inv_prec, out.ptr());
            if (!linearOpCheck) {
                success = false;
            }
        }

        //
        *out << "\nC) Create the forward solver using the created preconditioner ...\n";
        //

        const RCP<ParameterList> fwdSolverParamList = sublist(paramList, "Forward Solver");
        Stratimikos::DefaultLinearSolverBuilder fwdSolverSolverBuilder;
        fwdSolverSolverBuilder.setParameterList(fwdSolverParamList);
        const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > fwdSolverSolverStrategy
            = createLinearSolveStrategy(fwdSolverSolverBuilder);

        const RCP<Thyra::LinearOpWithSolveBase<double> >

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

void Piro::RythmosSolver<Scalar>::initialize(
#endif
    const Teuchos::RCP<Teuchos::ParameterList> &appParams,
    const Teuchos::RCP< Thyra::ModelEvaluator<Scalar> > &in_model,
    const Teuchos::RCP<Rythmos::IntegrationObserverBase<Scalar> > &observer)
{

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

    // set some internals
    model = in_model;
    num_p = in_model->Np();
    num_g = in_model->Ng();

    //
    *out << "\nA) Get the base parameter list ...\n";
    //


    if (appParams->isSublist("Rythmos")) {
        RCP<Teuchos::ParameterList> rythmosPL = sublist(appParams, "Rythmos", true);
        rythmosPL->validateParameters(*getValidRythmosParameters(),0);

        {
            const std::string verbosity = rythmosPL->get("Verbosity Level", "VERB_DEFAULT");
            if      (verbosity == "VERB_NONE")    solnVerbLevel = Teuchos::VERB_NONE;
            else if (verbosity == "VERB_DEFAULT") solnVerbLevel = Teuchos::VERB_DEFAULT;
            else if (verbosity == "VERB_LOW")     solnVerbLevel = Teuchos::VERB_LOW;
            else if (verbosity == "VERB_MEDIUM")  solnVerbLevel = Teuchos::VERB_MEDIUM;
            else if (verbosity == "VERB_HIGH")    solnVerbLevel = Teuchos::VERB_HIGH;
            else if (verbosity == "VERB_EXTREME") solnVerbLevel = Teuchos::VERB_EXTREME;
            else TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,"Unknown verbosity option specified in Piro_RythmosSolver.");
        }

        t_initial = rythmosPL->get("Initial Time", 0.0);
        t_final = rythmosPL->get("Final Time", 0.1);

        const std::string stepperType = rythmosPL->get("Stepper Type", "Backward Euler");

        //
        *out << "\nC) Create and initalize the forward model ...\n";
        //

        *out << "\nD) Create the stepper and integrator for the forward problem ...\n";
        //

        if (rythmosPL->get<std::string>("Nonlinear Solver Type") == "Rythmos") {
            Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> > rythmosTimeStepSolver =
                Rythmos::timeStepNonlinearSolver<Scalar>();
            if (rythmosPL->getEntryPtr("NonLinear Solver")) {
                RCP<Teuchos::ParameterList> nonlinePL =
                    sublist(rythmosPL, "NonLinear Solver", true);
                rythmosTimeStepSolver->setParameterList(nonlinePL);
            }
            fwdTimeStepSolver = rythmosTimeStepSolver;
        }
        else if (rythmosPL->get<std::string>("Nonlinear Solver Type") == "NOX") {
#ifdef HAVE_PIRO_NOX
            Teuchos::RCP<Thyra::NOXNonlinearSolver> nox_solver =  Teuchos::rcp(new Thyra::NOXNonlinearSolver);
            Teuchos::RCP<Teuchos::ParameterList> nox_params = Teuchos::rcp(new Teuchos::ParameterList);
            *nox_params = appParams->sublist("NOX");
            nox_solver->setParameterList(nox_params);
            fwdTimeStepSolver = nox_solver;
#else
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,"Requested NOX solver for a Rythmos Transient solve, Trilinos was not built with NOX enabled.  Please rebuild Trilinos or use the native Rythmos nonlinear solver.");
#endif

        }

        if (stepperType == "Backward Euler") {
            fwdStateStepper = Rythmos::backwardEulerStepper<Scalar> (model, fwdTimeStepSolver);
            fwdStateStepper->setParameterList(sublist(rythmosPL, "Rythmos Stepper", true));
        }
        else if (stepperType == "Forward Euler") {
            fwdStateStepper = Rythmos::forwardEulerStepper<Scalar> (model);
            fwdStateStepper->setParameterList(sublist(rythmosPL, "Rythmos Stepper", true));
        }
        else if (stepperType == "Explicit RK") {
            fwdStateStepper = Rythmos::explicitRKStepper<Scalar>(model);
            fwdStateStepper->setParameterList(sublist(rythmosPL, "Rythmos Stepper", true));
        }
        else if (stepperType == "BDF") {
            Teuchos::RCP<Teuchos::ParameterList> BDFparams =
                Teuchos::sublist(rythmosPL, "Rythmos Stepper", true);
            Teuchos::RCP<Teuchos::ParameterList> BDFStepControlPL =
                Teuchos::sublist(BDFparams,"Step Control Settings");

            fwdStateStepper = Teuchos::rcp( new Rythmos::ImplicitBDFStepper<Scalar>(model,fwdTimeStepSolver,BDFparams) );
            fwdStateStepper->setInitialCondition(model->getNominalValues());

        }
        else {
            // first (before failing) check to see if the user has added stepper factory
            typename std::map<std::string,Teuchos::RCP<Piro::RythmosStepperFactory<Scalar> > >::const_iterator
            stepFactItr = stepperFactories.find(stepperType);
            if(stepFactItr!=stepperFactories.end()) {
                // the user has added it, hot dog lets build a new stepper!
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      verbLevel = Teuchos::VERB_LOW;

    const Teuchos::EVerbosityLevel
      solnVerbLevel = ( dumpFinalSolutions ? Teuchos::VERB_EXTREME : verbLevel );

    //
    // Get the base parameter list that all other parameter lists will be read
    // from.
    //

    RCP<ParameterList>
      paramList = Teuchos::parameterList();
    if (paramsFileName.length())
      updateParametersFromXmlFile( paramsFileName, paramList.ptr() );
    if(extraParamsFile.length())
      Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile, paramList.ptr() );
    if (extraParamsString.length())
      updateParametersFromXmlString( extraParamsString, paramList.ptr() );

    if (testExactSensitivity) {
      paramList->sublist(DiagonalTransientModel_name).set("Exact Solution as Response",true);
    }

    paramList->validateParameters(*getValidParameters(),0); // Only validate top level lists!

    //
    // Create the Stratimikos linear solver factory.
    //
    // This is the linear solve strategy that will be used to solve for the
    // linear system with the W.
    //

    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
    linearSolverBuilder.setParameterList(sublist(paramList,Stratimikos_name));
    RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
      W_factory = createLinearSolveStrategy(linearSolverBuilder);

    //
    // Create the underlying EpetraExt::ModelEvaluator
    //

    RCP<EpetraExt::DiagonalTransientModel>
      epetraStateModel = EpetraExt::diagonalTransientModel(
        epetra_comm,
        sublist(paramList,DiagonalTransientModel_name)
        );

    *out <<"\nepetraStateModel valid options:\n";
    epetraStateModel->getValidParameters()->print(
      *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
      );

    //
    // Create the Thyra-wrapped ModelEvaluator
    //

    RCP<Thyra::ModelEvaluator<double> >
      stateModel = epetraModelEvaluator(epetraStateModel,W_factory);

    *out << "\nParameter names = " << *stateModel->get_p_names(0) << "\n";

    //
    // Create the Rythmos stateStepper
    //

    RCP<Rythmos::TimeStepNonlinearSolver<double> >

//.........这里部分代码省略.........
                       "The maximum allowed error in the integrated state in relation to the exact state solution" );

        // ToDo: Read in more parameters

        CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
        if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

        if ( Teuchos::VERB_DEFAULT == verbLevel )
            verbLevel = Teuchos::VERB_LOW;

        const Teuchos::EVerbosityLevel
        solnVerbLevel = ( dumpFinalSolutions ? Teuchos::VERB_EXTREME : verbLevel );

        //
        // B) Get the base parameter list that all other parameter lists will be
        // read from.
        //

        RCP<ParameterList> paramList = Teuchos::parameterList();
        if (paramsFileName.length())
            updateParametersFromXmlFile( paramsFileName, &*paramList );
        if (extraParamsString.length())
            updateParametersFromXmlString( extraParamsString, &*paramList );

        paramList->validateParameters(*getValidParameters());

        //
        // C) Create the Stratimikos linear solver factories.
        //

        // Get the linear solve strategy that will be used to solve for the linear
        // system with the dae's W matrix.
        Stratimikos::DefaultLinearSolverBuilder daeLinearSolverBuilder;
        daeLinearSolverBuilder.setParameterList(sublist(paramList,DAELinearSolver_name));
        RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
        daeLOWSF = createLinearSolveStrategy(daeLinearSolverBuilder);

        // Get the linear solve strategy that can be used to override the overall
        // linear system solve
        Stratimikos::DefaultLinearSolverBuilder overallLinearSolverBuilder;
        overallLinearSolverBuilder.setParameterList(sublist(paramList,OverallLinearSolver_name));
        RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
        overallLOWSF = createLinearSolveStrategy(overallLinearSolverBuilder);

        //
        // D) Create the underlying EpetraExt::ModelEvaluator
        //

        RCP<EpetraExt::DiagonalTransientModel> epetraDaeModel =
            EpetraExt::diagonalTransientModel(
                epetra_comm,
                sublist(paramList,DiagonalTransientModel_name)
            );

        *out <<"\nepetraDaeModel valid options:\n";
        epetraDaeModel->getValidParameters()->print(
            *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
        );

        //
        // E) Create the Thyra-wrapped ModelEvaluator
        //

        RCP<Thyra::ModelEvaluator<double> > daeModel =
            epetraModelEvaluator(epetraDaeModel,daeLOWSF);

// calls MPI_Init and MPI_Finalize
int main(int argc,char * argv[])
{
   using Teuchos::RCP;
   using panzer::StrPureBasisPair;
   using panzer::StrPureBasisComp;

   PHX::InitializeKokkosDevice();

   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
   Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
   out.setOutputToRootOnly(0);
   out.setShowProcRank(true);

   // variable declarations
   ////////////////////////////////////////////////////

   // factory definitions
   Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
   panzer_stk::SquareQuadMeshFactory mesh_factory;
   user_app::BCFactory bc_factory;

   // other declarations
   const std::size_t workset_size = 20;
   Teuchos::RCP<panzer::FieldManagerBuilder> fmb = 
         Teuchos::rcp(new panzer::FieldManagerBuilder);

   RCP<panzer_stk::STK_Interface> mesh;

   // construction of uncommitted (no elements) mesh 
   ////////////////////////////////////////////////////////

   // set mesh factory parameters
   RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
   pl->set("X Blocks",2);
   pl->set("Y Blocks",1);
   pl->set("X Elements",10);
   pl->set("Y Elements",10);
   mesh_factory.setParameterList(pl);
   mesh = mesh_factory.buildUncommitedMesh(MPI_COMM_WORLD);

   // construct input physics and physics block
   ////////////////////////////////////////////////////////

   out << "BUILD PHYSICS" << std::endl;
   Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList("Physics Blocks");
   std::vector<panzer::BC> bcs;
   std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;
   {
      std::map<std::string,std::string> block_ids_to_physics_ids;
      std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;

      testInitialzation(ipb, bcs);

      block_ids_to_physics_ids["eblock-0_0"] = "test physics";
      block_ids_to_physics_ids["eblock-1_0"] = "test physics";

      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");
   
      Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();

      int default_integration_order = 1;
      
      // build physicsBlocks map
      panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
                                 block_ids_to_cell_topo,
				 ipb,
				 default_integration_order,
				 workset_size,
			         eqset_factory,
				 gd,
			         true,
			         physicsBlocks);
   }

   // finish building mesh, set required field variables and mesh bulk data
   ////////////////////////////////////////////////////////////////////////

   {
      std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator physIter;
      for(physIter=physicsBlocks.begin();physIter!=physicsBlocks.end();++physIter) {
         Teuchos::RCP<const panzer::PhysicsBlock> pb = *physIter;
         const std::vector<StrPureBasisPair> & blockFields = pb->getProvidedDOFs();

         // insert all fields into a set
         std::set<StrPureBasisPair,StrPureBasisComp> fieldNames;
         fieldNames.insert(blockFields.begin(),blockFields.end());

         // add basis to DOF manager: block specific
         std::set<StrPureBasisPair,StrPureBasisComp>::const_iterator fieldItr;
         for (fieldItr=fieldNames.begin();fieldItr!=fieldNames.end();++fieldItr) {

            mesh->addSolutionField(fieldItr->first,pb->elementBlockID());
         }
      }

      mesh_factory.completeMeshConstruction(*mesh,MPI_COMM_WORLD);
   }

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

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

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

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

  using std::endl;
  typedef double Scalar;
  // typedef double ScalarMag; // unused
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Teuchos::describe;
  using Teuchos::Array;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::outArg;
  using Teuchos::rcp_implicit_cast;
  using Teuchos::rcp_dynamic_cast;
  using Teuchos::as;
  using Teuchos::ParameterList;
  using Teuchos::CommandLineProcessor;
  typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
  typedef Thyra::ModelEvaluatorBase MEB;

  bool result, success = true;

  Teuchos::GlobalMPISession mpiSession(&argc,&argv);

  RCP<Epetra_Comm> epetra_comm;
#ifdef HAVE_MPI
  epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
#else
  epetra_comm = rcp( new Epetra_SerialComm );
#endif // HAVE_MPI

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

  try {

    //
    // Read commandline options
    //

    CommandLineProcessor clp;
    clp.throwExceptions(false);
    clp.addOutputSetupOptions(true);

    std::string paramsFileName = "";
    clp.setOption( "params-file", &paramsFileName,
      "File name for XML parameters" );

    double t_final = 1e-3;
    clp.setOption( "final-time", &t_final,
      "Final integration time (initial time is 0.0)" );

    int numTimeSteps = 10;
    clp.setOption( "num-time-steps", &numTimeSteps,
      "Number of (fixed) time steps.  If <= 0.0, then variable time steps are taken" );

    double maxStateError = 1e-14;
    clp.setOption( "max-state-error", &maxStateError,
      "Maximum relative error in the integrated state allowed" );

    Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
    setVerbosityLevelOption( "verb-level", &verbLevel,
      "Top-level verbosity level.  By default, this gets deincremented as you go deeper into numerical objects.",
      &clp );

    Teuchos::EVerbosityLevel solnVerbLevel = Teuchos::VERB_DEFAULT;
    setVerbosityLevelOption( "soln-verb-level", &solnVerbLevel,
      "Solution verbosity level",
      &clp );

    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

    //
    *out << "\nA) Get the base parameter list ...\n";
    //

    RCP<ParameterList>
      paramList = Teuchos::parameterList();
    if (paramsFileName.length())
      updateParametersFromXmlFile( paramsFileName, paramList.ptr() );

    paramList->validateParameters(*getValidParameters());

    const Scalar t_init = 0.0;

    const Rythmos::TimeRange<Scalar> fwdTimeRange(t_init, t_final);
    const Scalar delta_t = t_final / numTimeSteps;
    *out << "\ndelta_t = " << delta_t;

    //
    *out << "\nB) Create the Stratimikos linear solver factory ...\n";
    //
    // This is the linear solve strategy that will be used to solve for the
    // linear system with the W.
    //

    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
    linearSolverBuilder.setParameterList(sublist(paramList,Stratimikos_name));
    RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
//.........这里部分代码省略.........

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

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