C++ LinearOpTester::check方法代码示例

TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( SimpleDenseLinearOp, basic,
  Scalar )
{

  using Teuchos::rcp_dynamic_cast;
  typedef ScalarTraits<Scalar> ST;

  const RCP<MultiVectorBase<Scalar> > mv =
    createSerialMultiVector<Scalar>(g_dim, g_dim/2, ST::one());
  const RCP<LinearOpBase<Scalar> > op =
    createNonconstSimpleDenseLinearOp<Scalar>(mv);

  TEST_EQUALITY(
    mv,
    rcp_dynamic_cast<SimpleDenseLinearOp<Scalar> >(op)->getNonconstMultiVector()
    );

  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.dump_all(g_dumpAll);
  TEST_ASSERT(linearOpTester.check(*op, ptrFromRef(out)));

}

bool runCgSolveExample(
  const int dim,
  const Scalar diagScale,
  const bool symOp,
  const bool showAllTests,
  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
  const int maxNumIters
  )
{

  using Teuchos::as;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::multiply;
  using Thyra::scale;
  typedef typename ST::magnitudeType  ScalarMag;
  bool success = true;
  bool result;
  Teuchos::RCP<Teuchos::FancyOStream> out =
    Teuchos::VerboseObjectBase::getDefaultOStream();
  *out << "\n***\n*** Running silly CG solver using scalar type = \'"
       << ST::name() << "\' ...\n***\n";
  Teuchos::Time timer("");
  timer.start(true);

  //
  // (A) Setup a simple linear system with tridiagonal operator:
  //
  //       [   a*2   -1                         ]
  //       [ -r(1)  a*2       -1                ]
  //  A =  [          .        .        .       ]
  //       [             -r(n-2)      a*2    -1 ]
  //       [                      -r(n-1)   a*2 ]
  //

  // (A.1) Create the tridiagonal matrix operator
  *out << "\nConstructing tridiagonal matrix A of dimension = " << dim
       << " and diagonal multiplier = " << diagScale << " ...\n";
  Teuchos::Array<Scalar> lower(dim-1), diag(dim), upper(dim-1);
  const Scalar
    up = -ST::one(),
    diagTerm = as<Scalar>(2.0) * diagScale * ST::one(),
    low = -(symOp ? ST::one() : ST::random());
  int k = 0;
  // First row
  diag[k] = diagTerm; upper[k] = up;
  // Middle rows
  for( k = 1; k < dim - 1; ++k ) {
    lower[k-1] = low; diag[k] = diagTerm; upper[k] = up;
  }
  // Last row
  lower[k-1] = low; diag[k] = diagTerm;
  RCP<const Thyra::LinearOpBase<Scalar> > A =
    rcp(new ExampleTridiagSerialLinearOp<Scalar>(dim, lower, diag, upper));

  // (A.2) Testing the linear operator constructed linear operator
  *out << "\nTesting the constructed linear operator A ...\n";
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.enable_all_tests(false);
  linearOpTester.check_linear_properties(true);
  linearOpTester.set_all_error_tol(tolerance);
  linearOpTester.set_all_warning_tol(1e-2*tolerance);
  linearOpTester.show_all_tests(showAllTests);
  result = linearOpTester.check(*A, out.ptr());
  if(!result) success = false;

  // (A.3) Create RHS vector b and set to a random value
  RCP<Thyra::VectorBase<Scalar> > b = createMember(A->range());
  Thyra::seed_randomize<Scalar>(0);
  Thyra::randomize( -ST::one(), +ST::one(), b.ptr() );

  // (A.4) Create LHS vector x and set to zero
  RCP<Thyra::VectorBase<Scalar> > x = createMember(A->domain());
  Thyra::V_S( x.ptr(), ST::zero() );

  // (A.5) Create the final linear system
  if(!symOp) {
    *out << "\nSetting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...\n";
    // A^H*A
    RCP<const Thyra::LinearOpBase<Scalar> > AtA = multiply(adjoint(A), A);
    // A^H*b
    RCP<Thyra::VectorBase<Scalar> > nb = createMember(AtA->range());
    Thyra::apply<Scalar>(*A, Thyra::CONJTRANS, *b, nb.ptr());
    A = AtA;
    b = nb;
  }

  // (A.6) Testing the linear operator used with the solve
  *out << "\nTesting the linear operator used with the solve ...\n";
  linearOpTester.check_for_symmetry(true);
  result = linearOpTester.check(*A, out.ptr());
  if(!result) success = false;

  //
  // (B) Solve the linear system with the silly CG solver
  //
  *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    E = mvE,
    F = mvF,
    J = mvJ,
    K = mvK,
    L = mvL,
    N = mvN,
    P = mvP,
    Q = mvQ;

  out << describe(*A, verbLevel);
  out << describe(*B, verbLevel);
  out << describe(*C, verbLevel);
  out << describe(*E, verbLevel);
  out << describe(*F, verbLevel);
  out << describe(*J, verbLevel);
  out << describe(*K, verbLevel);
  out << describe(*L, verbLevel);
  out << describe(*N, verbLevel);
  out << describe(*P, verbLevel);
  out << describe(*Q, verbLevel);

  //
  // B) Create the composed linear operators
  //

  // I
  const RCP<const LinearOpBase<Scalar> > I = identity(space1, "I");

  // D = diag(d)
  const RCP<const LinearOpBase<Scalar> > D = diagonal(d, "D");

  // M00 = [ gama*B*A + C,  E + F ] ^H
  //       [ J^H * A,       I     ]
  const RCP<const LinearOpBase<Scalar> > M00 =
    adjoint(
      block2x2(
        add( scale(gamma,multiply(B,A)), C ),  add( E, F ),
        multiply(adjoint(J),A),                I
        ),
      "M00"
      );

  out << "\nM00 = " << describe(*M00, verbLevel);

  // M01 = beta * [ Q ]
  //              [ K ]
  const RCP<const LinearOpBase<Scalar> > M01 =
    scale(
      beta,
      block2x1( Q, K ),
      "M01"
      );

  out << "\nM01 = "  << describe(*M01, verbLevel);
            
  // M10 = [ L * N^H,  eta*P ]
  const RCP<const LinearOpBase<Scalar> > M10 =
    block1x2(
      multiply(L,adjoint(N)),  scale(eta,P),
      "M10"
      );

  out << "\nM10 = " << describe(*M10, verbLevel);

  // M11 = D - Q^H*Q
  const RCP<const LinearOpBase<Scalar> > M11 =
    subtract( D, multiply(adjoint(Q),Q), "M11" );

  out << "\nM11 = "  << describe(*M11, verbLevel);
  

  // M = [ M00, M01 ]
  //     [ M10, M11 ]
  const RCP<const LinearOpBase<Scalar> > M =
    block2x2(
      M00, M01,
      M10, M11,
      "M"
      );

  out << "\nM = " << describe(*M, verbLevel);

  
  //
  // C) Test the final composed operator
  //

  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.set_all_error_tol(errorTol);
  linearOpTester.check_adjoint(testAdjoint);
  if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
    linearOpTester.show_all_tests(true);
  if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME))
    linearOpTester.dump_all(true);

  const bool result = linearOpTester.check(*M,&out);

  return result;

}

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

//.........这里部分代码省略.........
#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> >
        A_lows = fwdSolverSolverStrategy->createOp();

        initializePreconditionedOp<double>( *fwdSolverSolverStrategy, A,
                                            unspecifiedPrec(A_inv_prec), A_lows.ptr());
        //A_lows->setVerbLevel(verbLevel);
        *out << "\nA_lows = " << describe(*A_lows, verbLevel) << "\n";

        //
        *out << "\nD) Solve the linear system for a random RHS ...\n";
        //

        VectorPtr x = createMember(A->domain());
        VectorPtr b = createMember(A->range());
        Thyra::randomize(-1.0, +1.0, b.ptr());
        Thyra::assign(x.ptr(), 0.0); // Must give an initial guess!