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

TEUCHOS_UNIT_TEST( EpetraLinearOp, rectangular )
{
  using Teuchos::null;
  using Teuchos::inOutArg;
  using Teuchos::updateSuccess;

  const RCP<const Epetra_Comm> comm = getEpetraComm();
  const int numProcs = comm->NumProc();

  const int numLocalRows = g_localDim;
  const int numRows = numLocalRows * comm->NumProc();
  const int numCols = numLocalRows / 2;

  const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);

  const RCP<const LinearOpBase<double> > epetraOp = epetraLinearOp(epetraCrsM);

  LinearOpTester<double> linearOpTester;
  linearOpTester.check_adjoint(numProcs == 1);
  linearOpTester.show_all_tests(g_show_all_tests);
  linearOpTester.dump_all(g_dumpAll);
  updateSuccess(linearOpTester.check(*epetraOp, inOutArg(out)), success);

  // NOTE: Above, it would seem the Epetra_CrsMatrix::Apply(...) does not work
  // when doing and adjoint where the RowMap has empty processes.

}

SolveStatus<Scalar>
DefaultMultiVectorLinearOpWithSolve<Scalar>::solveImpl(
  const EOpTransp transp,
  const MultiVectorBase<Scalar> &BB,
  const Ptr<MultiVectorBase<Scalar> > &XX,
  const Ptr<const SolveCriteria<Scalar> > solveCriteria
  ) const
{

  using Teuchos::dyn_cast;
  using Teuchos::outArg;
  using Teuchos::inOutArg;
  typedef DefaultMultiVectorProductVector<Scalar> MVPV;

  const Ordinal numCols = BB.domain()->dim();

  SolveStatus<Scalar> overallSolveStatus;
  accumulateSolveStatusInit(outArg(overallSolveStatus));
  
  for (Ordinal col_j = 0; col_j < numCols; ++col_j) {

    const RCP<const VectorBase<Scalar> > b = BB.col(col_j);
    const RCP<VectorBase<Scalar> > x = XX->col(col_j);

    RCP<const MultiVectorBase<Scalar> >
      B = dyn_cast<const MVPV>(*b).getMultiVector().assert_not_null();
    RCP<MultiVectorBase<Scalar> >
      X = dyn_cast<MVPV>(*x).getNonconstMultiVector().assert_not_null();

    const SolveStatus<Scalar> solveStatus =
      Thyra::solve(*lows_.getConstObj(), transp, *B, X.ptr(), solveCriteria);

    accumulateSolveStatus(
      SolveCriteria<Scalar>(), // Never used
      solveStatus, inOutArg(overallSolveStatus) );

  }
  
  return overallSolveStatus;

}

  int
main (int argc, char *argv[])
{
  using Teuchos::inOutArg;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using std::endl;
  typedef double                          ST;
  typedef Epetra_Operator                 OP;
  typedef Epetra_MultiVector              MV;
  typedef Belos::OperatorTraits<ST,MV,OP> OPT;
  typedef Belos::MultiVecTraits<ST,MV>    MVT;

  // This calls MPI_Init and MPI_Finalize as necessary.
  Belos::Test::MPISession session (inOutArg (argc), inOutArg (argv));
  RCP<const Epetra_Comm> comm = session.getComm ();

  bool success = false;
  bool verbose = false;
  try {
    int MyPID = comm->MyPID ();

    //
    // Parameters to read from command-line processor
    //
    int frequency = -1;  // how often residuals are printed by solver
    int numRHS = 1;  // total number of right-hand sides to solve for
    int maxIters = 13000;  // maximum number of iterations for solver to use
    std::string filename ("bcsstk14.hb");
    double tol = 1.0e-5; // relative residual tolerance

    //
    // Read in command-line arguments
    //
    Teuchos::CommandLineProcessor cmdp (false, true);
    cmdp.setOption ("verbose", "quiet", &verbose, "Print messages and results.");
    cmdp.setOption ("frequency", &frequency, "Solvers frequency for printing "
        "residuals (#iters).");
    cmdp.setOption ("tol", &tol, "Relative residual tolerance used by MINRES "
        "solver.");
    cmdp.setOption ("filename", &filename, "Filename for Harwell-Boeing test "
        "matrix.");
    cmdp.setOption ("num-rhs", &numRHS, "Number of right-hand sides to solve.");
    cmdp.setOption ("max-iters", &maxIters, "Maximum number of iterations per "
        "linear system (-1 means \"adapt to problem/block size\").");
    if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      return EXIT_FAILURE;
    }
    Teuchos::oblackholestream blackHole;
    std::ostream& verbOut = (verbose && MyPID == 0) ? std::cout : blackHole;

    //
    // Generate the linear system(s) to solve.
    //
    verbOut << "Generating the linear system(s) to solve" << endl << endl;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    RCP<Epetra_Map> rowMap;
    try {
      // This might change the number of right-hand sides, if we read in
      // a right-hand side from the Harwell-Boeing file.
      Belos::Util::createEpetraProblem (filename, &rowMap, &A, &B, &X, &MyPID, numRHS);
    } catch (std::exception& e) {
      TEUCHOS_TEST_FOR_EXCEPTION (true, std::runtime_error,
          "Failed to create Epetra problem for matrix "
          "filename \"" << filename << "\".  "
          "createEpetraProblem() reports the following "
          "error: " << e.what());
    }
    //
    // Compute the initial residual norm of the problem, so we can see
    // by how much it improved after the solve.
    //
    std::vector<double> initialResidualNorms (numRHS);
    std::vector<double> initialResidualInfNorms (numRHS);
    Epetra_MultiVector R (*rowMap, numRHS);
    OPT::Apply (*A, *X, R);
    MVT::MvAddMv (-1.0, R, 1.0, *B, R); // R := -(A*X) + B.
    MVT::MvNorm (R, initialResidualNorms);
    MVT::MvNorm (R, initialResidualInfNorms, Belos::InfNorm);
    if (verbose) {
      verbOut << "Initial residual 2-norms:            \t";
      for (int i = 0; i < numRHS; ++i) {
        verbOut << initialResidualNorms[i];
        if (i < numRHS-1) {
          verbOut << ", ";
        }
      }
      verbOut << endl << "Initial residual Inf-norms:          \t";
      for (int i = 0; i < numRHS; ++i) {
        verbOut << initialResidualInfNorms[i];
        if (i < numRHS-1) {
          verbOut << ", ";
        }
      }
      verbOut << endl;
    }
//.........这里部分代码省略.........

bool BrentsLineSearch<Scalar>::doLineSearch(
  const MeritFunc1DBase<Scalar> &phi,
  const PointEval1D<Scalar> &point_k,
  const Ptr<PointEval1D<Scalar> > &point_kp1,
  const Ptr<int> &numIters
  ) const
{

  using Teuchos::as;
  using Teuchos::OSTab;
  using Teuchos::outArg;
  using Teuchos::inOutArg;
  typedef ScalarTraits<Scalar> ST;
  typedef PointEval1D<Scalar> PE1D;

#ifdef TEUCHOS_DEBUG
  TEUCHOS_ASSERT_EQUALITY(point_k.alpha, ST::zero());
  TEUCHOS_ASSERT_INEQUALITY(point_k.phi, !=, PE1D::valNotGiven());
  TEUCHOS_ASSERT_EQUALITY(point_k.Dphi, PE1D::valNotGiven());
  TEUCHOS_ASSERT(!is_null(point_kp1));
  TEUCHOS_ASSERT_INEQUALITY(point_kp1->alpha, >, ST::zero());
  TEUCHOS_ASSERT_INEQUALITY(point_kp1->phi, !=, PE1D::valNotGiven());
  TEUCHOS_ASSERT_EQUALITY(point_kp1->Dphi, PE1D::valNotGiven());
#endif

  const RCP<Teuchos::FancyOStream> out = this->getOStream();
  bracket_.setOStream(out);
  brentsMin_.setOStream(out);

  *out << "\nStarting bracketing and brents 1D minimization linesearch ...\n";

  OSTab tab(out);

  int totalNumIters = 0;

  PointEval1D<Scalar> p_l = point_k;
  PointEval1D<Scalar> &p_m = *point_kp1; // Update in place!
  PointEval1D<Scalar> p_u;

  bool success = true;

  // A) Bracket the minimum

  int numBracketIters = -1;

  const bool bracketSuccess = bracket_.bracketMinimum(
    phi, inOutArg(p_l), inOutArg(p_m), outArg(p_u), outArg(numBracketIters) );

  if (!bracketSuccess) success = false;

  totalNumIters += numBracketIters;

  // B) Do approximate mimimization in the bracket

  if (bracketSuccess) {

    int numBrentsIters = -1;
    
    const bool brentsSuccess = brentsMin_.approxMinimize(
      phi, p_l, inOutArg(p_m), p_u, outArg(numBrentsIters) );
    
    if (!brentsSuccess) success = false;
    
    totalNumIters += numBrentsIters;

  }

  // C) Overall success?

  if (!is_null(numIters))
    *numIters = totalNumIters;

  return success;
  
}

NonlinearCGUtils::ESolveReturn
NonlinearCG<Scalar>::doSolve(
  const Ptr<Thyra::VectorBase<Scalar> > &p_inout,
  const Ptr<ScalarMag> &g_opt_out,
  const Ptr<const ScalarMag> &g_reduct_tol_in,
  const Ptr<const ScalarMag> &g_grad_tol_in,
  const Ptr<const ScalarMag> &alpha_init_in,
  const Ptr<int> &numIters_out
  )
{

  typedef ScalarTraits<Scalar> ST;
  typedef ScalarTraits<ScalarMag> SMT;
  
  using Teuchos::null;
  using Teuchos::as;
  using Teuchos::tuple;
  using Teuchos::rcpFromPtr;
  using Teuchos::optInArg;
  using Teuchos::inOutArg;
  using GlobiPack::computeValue;
  using GlobiPack::PointEval1D;
  using Thyra::VectorSpaceBase;
  using Thyra::VectorBase;
  using Thyra::MultiVectorBase;
  using Thyra::scalarProd;
  using Thyra::createMember;
  using Thyra::createMembers;
  using Thyra::get_ele;
  using Thyra::norm;
  using Thyra::V_StV;
  using Thyra::Vt_S;
  using Thyra::eval_g_DgDp;
  typedef Thyra::Ordinal Ordinal;
  typedef Thyra::ModelEvaluatorBase MEB;
  namespace NCGU = NonlinearCGUtils;
  using std::max;

  // Validate input

  g_opt_out.assert_not_null();

  // Set streams

  const RCP<Teuchos::FancyOStream> out = this->getOStream();
  linesearch_->setOStream(out);

  // Determine what step constants will be computed

  const bool compute_beta_PR =
    (
      solverType_ == NCGU::NONLINEAR_CG_PR_PLUS
      ||
      solverType_ == NCGU::NONLINEAR_CG_FR_PR
      );

  const bool compute_beta_HS = (solverType_ == NCGU::NONLINEAR_CG_HS);

  //
  // A) Set up the storage for the algorithm
  //
  
  const RCP<DefaultPolyLineSearchPointEvaluator<Scalar> >
    pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();

  const RCP<UnconstrainedOptMeritFunc1D<Scalar> >
    meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
      model_, paramIndex_, responseIndex_ );

  const RCP<const VectorSpaceBase<Scalar> >
    p_space = model_->get_p_space(paramIndex_),
    g_space = model_->get_g_space(responseIndex_);

  // Stoarge for current iteration
  RCP<VectorBase<Scalar> >
    p_k = rcpFromPtr(p_inout),        // Current solution for p
    p_kp1 = createMember(p_space),    // Trial point for p (in line search)
    g_vec = createMember(g_space),    // Vector (size 1) form of objective g(p) 
    g_grad_k = createMember(p_space), // Gradient of g DgDp^T
    d_k = createMember(p_space),      // Search direction
    g_grad_k_diff_km1 = null;         // g_grad_k - g_grad_km1 (if needed)

  // Storage for previous iteration
  RCP<VectorBase<Scalar> >
    g_grad_km1 = null, // Will allocate if we need it!
    d_km1 = null; // Will allocate if we need it!
  ScalarMag
    alpha_km1 = SMT::zero(),
    g_km1 = SMT::zero(),
    g_grad_km1_inner_g_grad_km1 = SMT::zero(),
    g_grad_km1_inner_d_km1 = SMT::zero();
  
  if (compute_beta_PR || compute_beta_HS) {
    g_grad_km1 = createMember(p_space);
    g_grad_k_diff_km1 = createMember(p_space);
  }
  
  if (compute_beta_HS) {
    d_km1 = createMember(p_space);
  }
//.........这里部分代码省略.........

void VanderPolModel::evalModelImpl(
  const ModelEvaluatorBase::InArgs<double> &inArgs,
  const ModelEvaluatorBase::OutArgs<double> &outArgs
  ) const
{

  using Teuchos::as;
  using Teuchos::outArg;
  using Teuchos::optInArg;
  using Teuchos::inOutArg;
  using Sacado::Fad::DFad;

  TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
      "Error, setParameterList must be called first!\n"
      );

  const RCP<const VectorBase<double> > x_in = inArgs.get_x().assert_not_null();
  Thyra::ConstDetachedVectorView<double> x_in_view( *x_in ); 

  double t = inArgs.get_t();
  double eps = epsilon_;
  if (acceptModelParams_) {
    const RCP<const VectorBase<double> > p_in = inArgs.get_p(0).assert_not_null();
    Thyra::ConstDetachedVectorView<double> p_in_view( *p_in ); 
    eps = p_in_view[0];
  }

  RCP<const VectorBase<double> > x_dot_in;
  double alpha = -1.0;
  double beta = -1.0;
  if (isImplicit_) {
    x_dot_in = inArgs.get_x_dot().assert_not_null();
    alpha = inArgs.get_alpha();
    beta = inArgs.get_beta();
  }

  const RCP<VectorBase<double> > f_out = outArgs.get_f();
  const RCP<Thyra::LinearOpBase<double> > W_out = outArgs.get_W_op();
  RCP<Thyra::MultiVectorBase<double> > DfDp_out; 
  if (acceptModelParams_) {
    Derivative<double> DfDp = outArgs.get_DfDp(0); 
    DfDp_out = DfDp.getMultiVector();
  }

  // Determine how many derivatives we will compute
  
  int num_derivs = 0;
  if (nonnull(W_out)) {
    num_derivs += 2;
    if (isImplicit_) {
      num_derivs += 2;
    }
  }
  if (nonnull(DfDp_out))
    num_derivs += 1;

  // Set up the FAD derivative objects

  int deriv_i = 0;

  Array<DFad<double> > x_dot_fad;
  int x_dot_idx_offset = 0;
  if (isImplicit_) {
    Thyra::ConstDetachedVectorView<double> x_dot_in_view( *x_dot_in );
    if (nonnull(W_out)) {
      x_dot_idx_offset = deriv_i;
      x_dot_fad = convertToIndepVarFadArray<double>(x_dot_in_view.sv().values()(),
        num_derivs, inOutArg(deriv_i));
    }
    else {
      x_dot_fad = convertToPassiveFadArray<double>(x_dot_in_view.sv().values()());
    }
  }

  Array<DFad<double> > x_fad;
  int x_idx_offset = 0;
  if (nonnull(W_out)) {
    x_idx_offset = deriv_i;
    x_fad = convertToIndepVarFadArray<double>(x_in_view.sv().values()(),
      num_derivs, inOutArg(deriv_i));
  }
  else {
    x_fad = convertToPassiveFadArray<double>(x_in_view.sv().values()());
  }

  DFad<double> eps_fad(eps); // Default passive
  int eps_idx_offset = 0;
  if (nonnull(DfDp_out)) {
    eps_idx_offset = deriv_i;
    eps_fad = DFad<double>(num_derivs, deriv_i++, eps);
  }
  
  // Compute the function

  Array<DFad<double> > f_fad(2);
  this->eval_f<DFad<double> >( x_dot_fad, x_fad, eps_fad, t, f_fad );

  // Extract the output

  if (nonnull(f_out)) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   const RCP<const LinearOpBase<double> > A =
     Thyra::block2x2<double>(
       Thyra::epetraLinearOp(F),
       Thyra::epetraLinearOp(Bt),
       Thyra::epetraLinearOp(B),
       Thyra::epetraLinearOp(C),
       "A"
       );

   const RCP<Thyra::EpetraOperatorWrapper> epetra_A =
     rcp(new Thyra::EpetraOperatorWrapper(A));

   // begin the tests!
   const Epetra_Map & rangeMap  = epetra_A->OperatorRangeMap();
   const Epetra_Map & domainMap = epetra_A->OperatorDomainMap();

   // check to see that the number of global elements is correct
   TEST_EQUALITY(rangeMap.NumGlobalElements(), 2*nx*ny);
   TEST_EQUALITY(domainMap.NumGlobalElements(), 2*nx*ny);

   // largest global ID should be one less then the # of elements
   TEST_EQUALITY(rangeMap.NumGlobalElements()-1, rangeMap.MaxAllGID());
   TEST_EQUALITY(domainMap.NumGlobalElements()-1, domainMap.MaxAllGID());

   // create a vector to test: copyThyraIntoEpetra
   {
      const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
      Thyra::randomize(-100.0, 100.0, tv.ptr());
      const RCP<const VectorBase<double> > tv_0 =
        Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
      const RCP<const VectorBase<double> > tv_1 =
        Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
      const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
      const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);

      int off_0 = vv_0.globalOffset();
      int off_1 = vv_1.globalOffset();
      
      // create its Epetra counter part
      Epetra_Vector ev(epetra_A->OperatorDomainMap());
      epetra_A->copyThyraIntoEpetra(*tv, ev);

      // compare handle_tv to ev!
      TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
      const int numMyElements = domainMap.NumMyElements();
      double tval = 0.0;
      for(int i=0; i < numMyElements; i++) {
         int gid = domainMap.GID(i);
         if(gid<nx*ny)
            tval = vv_0[gid-off_0];
         else
            tval = vv_1[gid-off_1-nx*ny];
         TEST_EQUALITY(ev[i], tval);
      }
   }

   // create a vector to test: copyEpetraIntoThyra
   {
      // create an Epetra vector
     Epetra_Vector ev(epetra_A->OperatorDomainMap());
     ev.Random();

      // create its thyra counterpart
      const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
      const RCP<const VectorBase<double> > tv_0 =
        Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
      const RCP<const VectorBase<double> > tv_1 =
        Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
      const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
      const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);

      int off_0 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
        tv_0->space())->localOffset();
      int off_1 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
        tv_1->space())->localOffset();

      epetra_A->copyEpetraIntoThyra(ev, tv.ptr());
   
      // compare tv to ev!
      TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
      int numMyElements = domainMap.NumMyElements();
      double tval = 0.0;
      for(int i=0;i<numMyElements;i++) {
         int gid = domainMap.GID(i);
         if(gid<nx*ny)
            tval = vv_0[gid-off_0];
         else
            tval = vv_1[gid-off_1-nx*ny];
         TEST_EQUALITY(ev[i], tval);
      }
   }

   // Test using Thyra::LinearOpTester
   const RCP<const LinearOpBase<double> > thyraEpetraOp = epetraLinearOp(epetra_A);
   LinearOpTester<double> linearOpTester;
   linearOpTester.show_all_tests(true);
   const bool checkResult = linearOpTester.check(*thyraEpetraOp, inOutArg(out));
   TEST_ASSERT(checkResult);

}