本文整理汇总了C++中Thyra类的典型用法代码示例。如果您正苦于以下问题:C++ Thyra类的具体用法?C++ Thyra怎么用?C++ Thyra使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Thyra类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: V_S
void DiagonalROME<Scalar>::setDiagonalBarVector(
const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
{
typedef Teuchos::ScalarTraits<Scalar> ST;
using Teuchos::rcp_dynamic_cast;
using Thyra::createMember;
using Thyra::ele_wise_divide;
using Thyra::V_S;
diag_bar_ = diag_bar.assert_not_null();
// Reset the scalar product for p_space!
RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
// s_bar[i] = diag[i] / diag_bar[i]
V_S( s_bar.ptr(), ST::zero() );
ele_wise_divide( ST::one(), *diag_, *diag_bar_, s_bar.ptr() );
s_bar_ = s_bar;
const RCP<Thyra::ScalarProdVectorSpaceBase<Scalar> > sp_p_space =
rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
//sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));
}
开发者ID:00liujj,项目名称:trilinos,代码行数:26,代码来源:Diagonal_ThyraROME_def.hpp
示例2: sillyCgSolve
bool sillyCgSolve(
const Thyra::LinearOpBase<Scalar> &A,
const Thyra::VectorBase<Scalar> &b,
const int maxNumIters,
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x,
std::ostream &out
)
{
// Create some typedefs and some other stuff to make the code cleaner
typedef Teuchos::ScalarTraits<Scalar> ST; typedef typename ST::magnitudeType ScalarMag;
const Scalar one = ST::one(), zero = ST::zero(); using Teuchos::as;
using Teuchos::RCP; using Thyra::VectorSpaceBase; using Thyra::VectorBase;
using Thyra::NOTRANS; using Thyra::V_V; using Thyra::apply;
// Validate input
THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()", A, Thyra::NOTRANS, *x, &b);
Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM;
out << "\nStarting CG solver ...\n" << std::scientific << "\ndescribe A:\n"<<describe(A, vl)
<< "\ndescribe b:\n"<<describe(b, vl)<<"\ndescribe x:\n"<<describe(*x, vl)<<"\n";
// Initialization
const RCP<const VectorSpaceBase<Scalar> > space = A.domain();
const RCP<VectorBase<Scalar> > r = createMember(space);
// r = -A*x + b
V_V(r.ptr(), b); apply<Scalar>(A, NOTRANS, *x, r.ptr(), -one, one);
const ScalarMag r0_nrm = norm(*r);
if (r0_nrm==zero) return true;
const RCP<VectorBase<Scalar> > p = createMember(space), q = createMember(space);
Scalar rho_old = -one;
// Perform the iterations
for( int iter = 0; iter <= maxNumIters; ++iter ) {
// Check convergence and output iteration
const ScalarMag r_nrm = norm(*r);
const bool isConverged = r_nrm/r0_nrm <= tolerance;
if( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) {
out << "Iter = " << iter << ", ||b-A*x||/||b-A*x0|| = " << (r_nrm/r0_nrm) << std::endl;
if( r_nrm/r0_nrm < tolerance ) return true; // Success!
}
// Compute iteration
const Scalar rho = inner(*r, *r); // <r,r> -> rho
if (iter==0) V_V(p.ptr(), *r); // r -> p (iter == 0)
else Vp_V( p.ptr(), *r, rho/rho_old ); // r+(rho/rho_old)*p -> p (iter > 0)
apply<Scalar>(A, NOTRANS, *p, q.ptr()); // A*p -> q
const Scalar alpha = rho/inner(*p, *q); // rho/<p,q> -> alpha
Vp_StV( x, +alpha, *p ); // +alpha*p + x -> x
Vp_StV( r.ptr(), -alpha, *q ); // -alpha*q + r -> r
rho_old = rho; // rho -> rho_old (for next iter)
}
return false; // Failure
} // end sillyCgSolve
开发者ID:00liujj,项目名称:trilinos,代码行数:60,代码来源:sillyCgSolve.hpp
示例3: ostab
void ImplicitBDFStepperRampingStepControl<Scalar>::initialize(
const StepperBase<Scalar>& stepper)
{
// Initialize can be called from the stepper when setInitialCondition
// is called.
using Teuchos::as;
typedef Teuchos::ScalarTraits<Scalar> ST;
using Thyra::createMember;
// Set initial time:
TimeRange<Scalar> stepperRange = stepper.getTimeRange();
TEUCHOS_TEST_FOR_EXCEPTION(
!stepperRange.isValid(),
std::logic_error,
"Error, Stepper does not have valid time range for initialization "
"of ImplicitBDFStepperRampingStepControl!\n");
if (is_null(parameterList_)) {
RCP<Teuchos::ParameterList> emptyParameterList =
Teuchos::rcp(new Teuchos::ParameterList);
this->setParameterList(emptyParameterList);
}
if (is_null(errWtVecCalc_)) {
RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc =
rcp(new ImplicitBDFStepperErrWtVecCalc<Scalar>());
errWtVecCalc_ = IBDFErrWtVecCalc;
}
stepControlState_ = UNINITIALIZED;
requestedStepSize_ = Scalar(-1.0);
currentStepSize_ = initialStepSize_;
currentOrder_ = 1;
nextStepSize_ = initialStepSize_;
nextOrder_ = 1;
numberOfSteps_ = 0;
totalNumberOfFailedSteps_ = 0;
countOfConstantStepsAfterFailure_ = 0;
if (is_null(delta_)) {
delta_ = createMember(stepper.get_x_space());
}
if (is_null(errWtVec_)) {
errWtVec_ = createMember(stepper.get_x_space());
}
V_S(delta_.ptr(),ST::zero());
if ( doOutput_(Teuchos::VERB_HIGH) ) {
RCP<Teuchos::FancyOStream> out = this->getOStream();
Teuchos::OSTab ostab(out,1,"initialize");
*out << "currentOrder_ = " << currentOrder_ << std::endl;
*out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
}
setStepControlState_(BEFORE_FIRST_STEP);
}
开发者ID:gitter-badger,项目名称:quinoa,代码行数:58,代码来源:Rythmos_ImplicitBDFStepperRampingStepControl_def.hpp
示例4: createModel
const Teuchos::RCP<TriKota::DiagonalROME<Scalar> >
createModel(
const int globalDim,
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType &g_offset
)
{
using Teuchos::RCP;
const RCP<const Teuchos::Comm<Thyra::Ordinal> > comm =
Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
const int numProcs = comm->getSize();
TEUCHOS_TEST_FOR_EXCEPT_MSG( numProcs > globalDim,
"Error, the number of processors can not be greater than the global"
" dimension of the vectors!." );
const int localDim = globalDim / numProcs;
const int localDimRemainder = globalDim % numProcs;
TEUCHOS_TEST_FOR_EXCEPT_MSG( localDimRemainder != 0,
"Error, the number of processors must divide into the global number"
" of elements exactly for now!." );
const RCP<TriKota::DiagonalROME<Scalar> > model =
Teuchos::rcp(new TriKota::DiagonalROME<Scalar>(localDim));
const RCP<const Thyra::VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
const RCP<Thyra::VectorBase<Scalar> > ps = createMember(p_space);
const Scalar ps_val = 2.0;
Thyra::V_S(ps.ptr(), ps_val);
model->setSolutionVector(ps);
model->setScalarOffset(g_offset);
return model;
}
开发者ID:00liujj,项目名称:trilinos,代码行数:32,代码来源:Diagonal_ThyraROME_def.hpp
示例5:
void DefaultPolyLineSearchPointEvaluator<Scalar>::computePoint( const ScalarMag &alpha,
const Ptr<Thyra::VectorBase<Scalar> > &p
) const
{
typedef ScalarTraits<Scalar> ST;
using Teuchos::as;
using Thyra::V_V;
using Thyra::Vp_StV;
V_V( p, *vecs_[0] );
if (alpha != ST::zero()) {
ScalarMag alpha_i = alpha;
const int n = vecs_.size();
for (int i = 1; i < n; ++i, alpha_i *= alpha) {
Vp_StV(p, alpha_i, *vecs_[i]);
}
}
}
开发者ID:haripandey,项目名称:trilinos,代码行数:17,代码来源:OptiPack_DefaultPolyLineSearchPointEvaluator.hpp
示例6: buildInverseMassMatrix
void ExplicitModelEvaluator<Scalar>::
buildInverseMassMatrix() const
{
typedef Thyra::ModelEvaluatorBase MEB;
using Teuchos::RCP;
using Thyra::createMember;
RCP<const Thyra::ModelEvaluator<Scalar> > me = this->getUnderlyingModel();
// first allocate space for the mass matrix
RCP<Thyra::LinearOpBase<Scalar> > mass = me->create_W_op();
// intialize a zero to get rid of the x-dot
if(zero_==Teuchos::null) {
zero_ = Thyra::createMember(*me->get_x_space());
Thyra::assign(zero_.ptr(),0.0);
}
// request only the mass matrix from the physics
// Model evaluator builds: alpha*u_dot + beta*F(u) = 0
MEB::InArgs<Scalar> inArgs = me->createInArgs();
inArgs.set_x(createMember(me->get_x_space()));
inArgs.set_x_dot(zero_);
inArgs.set_alpha(-1.0);
inArgs.set_beta(0.0);
// set the one time beta to ensure dirichlet conditions
// are correctly included in the mass matrix: do it for
// both epetra and Tpetra. If a panzer model evaluator has
// not been passed in...oh well you get what you asked for!
if(panzerModel_!=Teuchos::null)
panzerModel_->setOneTimeDirichletBeta(-1.0);
else if(panzerEpetraModel_!=Teuchos::null)
panzerEpetraModel_->setOneTimeDirichletBeta(-1.0);
// set only the mass matrix
MEB::OutArgs<Scalar> outArgs = me->createOutArgs();
outArgs.set_W_op(mass);
// this will fill the mass matrix operator
me->evalModel(inArgs,outArgs);
if(!massLumping_) {
invMassMatrix_ = Thyra::inverse<Scalar>(*me->get_W_factory(),mass);
}
else {
// build lumped mass matrix (assumes all positive mass entries, does a simple sum)
Teuchos::RCP<Thyra::VectorBase<Scalar> > ones = Thyra::createMember(*mass->domain());
Thyra::assign(ones.ptr(),1.0);
RCP<Thyra::VectorBase<Scalar> > invLumpMass = Thyra::createMember(*mass->range());
Thyra::apply(*mass,Thyra::NOTRANS,*ones,invLumpMass.ptr());
Thyra::reciprocal(*invLumpMass,invLumpMass.ptr());
invMassMatrix_ = Thyra::diagonal(invLumpMass);
}
}
开发者ID:gitter-badger,项目名称:quinoa,代码行数:57,代码来源:Panzer_ExplicitModelEvaluator_impl.hpp
示例7: Ng_
DiagonalROME<Scalar>::DiagonalROME(
const int localDim,
const RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm
)
:Np_(1), Ng_(1), comm_(comm), localDim_(localDim),
nonlinearTermFactor_(0.0), g_offset_(0.0)
{
typedef Teuchos::ScalarTraits<Scalar> ST;
using Thyra::createMember;
TEUCHOS_ASSERT( localDim > 0 );
// Get the comm
if (is_null(comm_)) {
comm_ = Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
}
// Locally replicated space for g
g_space_ = Thyra::locallyReplicatedDefaultSpmdVectorSpace<Scalar>(comm_, 1);
// Distributed space for p
p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, localDim, -1);
// Default solution
const RCP<Thyra::VectorBase<Scalar> > ps = createMember<Scalar>(p_space_);
V_S(ps.ptr(), ST::zero());
ps_ = ps;
// Default diagonal
const RCP<Thyra::VectorBase<Scalar> > diag = createMember<Scalar>(p_space_);
V_S(diag.ptr(), ST::one());
diag_ = diag;
diag_bar_ = diag;
// Default inner product
const RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
V_S(s_bar.ptr(), ST::one());
s_bar_ = s_bar;
// Default response offset
g_offset_ = ST::zero();
}
开发者ID:00liujj,项目名称:trilinos,代码行数:44,代码来源:Diagonal_ThyraROME_def.hpp
示例8: assertBlockFillIsActive
bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::opSupportedImpl(
EOpTransp M_trans
) const
{
using Thyra::opSupported;
assertBlockFillIsActive(false);
for ( int k = 0; k < numDiagBlocks_; ++k ) {
if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
return false;
}
return true;
// ToDo: To be safe we really should do a collective reduction of all
// clusters of processes. However, for the typical use case, every block
// will return the same info and we should be safe!
}
开发者ID:00liujj,项目名称:trilinos,代码行数:15,代码来源:Thyra_DefaultBlockedTriangularLinearOpWithSolve_def.hpp
示例9: x_space_
Simple2DModelEvaluator<Scalar>::Simple2DModelEvaluator()
: x_space_(Thyra::defaultSpmdVectorSpace<Scalar>(2)),
f_space_(x_space_),
W_factory_(Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>()),
d_(0.0),
p_(x_space_->dim(), Teuchos::ScalarTraits<Scalar>::zero()),
showGetInvalidArg_(false)
{
using Teuchos::RCP;
using Thyra::VectorBase;
using Thyra::createMember;
typedef Thyra::ModelEvaluatorBase MEB;
typedef Teuchos::ScalarTraits<Scalar> ST;
MEB::InArgsSetup<Scalar> inArgs;
inArgs.setModelEvalDescription(this->description());
inArgs.setSupports(MEB::IN_ARG_x);
prototypeInArgs_ = inArgs;
MEB::OutArgsSetup<Scalar> outArgs;
outArgs.setModelEvalDescription(this->description());
outArgs.setSupports(MEB::OUT_ARG_f);
outArgs.setSupports(MEB::OUT_ARG_W_op);
outArgs.setSupports(MEB::OUT_ARG_W_prec);
prototypeOutArgs_ = outArgs;
nominalValues_ = inArgs;
x0_ = createMember(x_space_);
V_S(x0_.ptr(), ST::zero());
nominalValues_.set_x(x0_);
set_d(10.0);
set_p(Teuchos::tuple<Scalar>(2.0, 0.0)());
set_x0(Teuchos::tuple<Scalar>(1.0, 1.0)());
}
开发者ID:00liujj,项目名称:trilinos,代码行数:37,代码来源:Thyra_Simple2DModelEvaluator_def.hpp
示例10: tab
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);
}
//.........这里部分代码省略.........
开发者ID:gitter-badger,项目名称:quinoa,代码行数:101,代码来源:OptiPack_NonlinearCG_def.hpp
示例11: runCgSolveExample
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";
//.........这里部分代码省略.........
开发者ID:00liujj,项目名称:trilinos,代码行数:101,代码来源:sillyCgSolve_serial.cpp
示例12: main
int main(int argc, char *argv[])
{
using std::endl;
typedef double Scalar;
typedef double ScalarMag;
using Teuchos::describe;
using Teuchos::RCP;
using Teuchos::rcp;
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;
using Thyra::createMember;
using Thyra::createMembers;
bool 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 {
//
// A) Read commandline options
//
CommandLineProcessor clp;
clp.throwExceptions(false);
clp.addOutputSetupOptions(true);
std::string paramsFileName = "";
clp.setOption( "params-file", ¶msFileName,
"File name for XML parameters" );
std::string extraParamsString = "";
clp.setOption( "extra-params", &extraParamsString,
"Extra XML parameter string" );
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 );
double finalTime = 1.0;
clp.setOption( "final-time", &finalTime, "Final time (the inital time)" );
int numTimeSteps = 2;
clp.setOption( "num-time-steps", &numTimeSteps, "Number of time steps" );
bool dumpFinalSolutions = false;
clp.setOption(
"dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
"Determine if the final solutions are dumpped or not." );
double maxStateError = 1e-6;
clp.setOption( "max-state-error", &maxStateError,
"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
//.........这里部分代码省略.........
开发者ID:haripandey,项目名称:trilinos,代码行数:101,代码来源:timeDiscretizedBackwardEulerMain.cpp
示例13: run_composite_linear_ops_tests
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());
//.........这里部分代码省略.........
开发者ID:agrippa,项目名称:Trilinos,代码行数:101,代码来源:test_composite_linear_ops.cpp
示例14: TEUCHOS_TEST_FOR_EXCEPTION
void CubicSplineInterpolator<Scalar>::interpolate(
const Array<Scalar> &t_values,
typename DataStore<Scalar>::DataStoreVector_t *data_out
) const
{
using Teuchos::as;
using Teuchos::outArg;
typedef Teuchos::ScalarTraits<Scalar> ST;
TEUCHOS_TEST_FOR_EXCEPTION( nodesSet_ == false, std::logic_error,
"Error!, setNodes must be called before interpolate"
);
#ifdef HAVE_RYTHMOS_DEBUG
// Check that our nodes_ have not changed between the call to setNodes and interpolate
assertNodesUnChanged<Scalar>(*nodes_,*nodes_copy_);
// Assert that the base interpolator preconditions are satisfied
assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
#endif // HAVE_RYTHMOS_DEBUG
// Output info
const RCP<FancyOStream> out = this->getOStream();
const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
Teuchos::OSTab ostab(out,1,"CSI::interpolator");
if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
*out << "nodes_:" << std::endl;
for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
*out << "nodes_[" << i << "] = " << std::endl;
(*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
}
*out << "t_values = " << std::endl;
for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
*out << "t_values[" << i << "] = " << t_values[i] << std::endl;
}
}
data_out->clear();
// Return immediately if no time points are requested ...
if (t_values.size() == 0) {
return;
}
if ((*nodes_).size() == 1) {
// trivial case of one node. Preconditions assert that t_values[0] ==
// (*nodes_)[0].time so we can just pass it out
DataStore<Scalar> DS((*nodes_)[0]);
data_out->push_back(DS);
}
else { // (*nodes_).size() >= 2
int n = 0; // index into t_values
// Loop through all of the time interpolation points in the buffer and
// satisfiy all of the requested time points that you find. NOTE: The
// loop will be existed once all of the time points are satisified (see
// return below).
for (Teuchos::Ordinal i=0 ; i < (*nodes_).size()-1; ++i) {
const Scalar& ti = (*nodes_)[i].time;
const Scalar& tip1 = (*nodes_)[i+1].time;
const TimeRange<Scalar> range_i(ti,tip1);
// For the interpolation range of [ti,tip1], satisify all of the
// requested points in this range.
while ( range_i.isInRange(t_values[n]) ) {
// First we check for exact node matches:
if (compareTimeValues(t_values[n],ti)==0) {
DataStore<Scalar> DS((*nodes_)[i]);
data_out->push_back(DS);
}
else if (compareTimeValues(t_values[n],tip1)==0) {
DataStore<Scalar> DS((*nodes_)[i+1]);
data_out->push_back(DS);
} else {
if (!splineCoeffComputed_) {
computeCubicSplineCoeff<Scalar>(*nodes_,outArg(splineCoeff_));
splineCoeffComputed_ = true;
}
DataStore<Scalar> DS;
RCP<Thyra::VectorBase<Scalar> > x = createMember((*nodes_)[i].x->space());
evaluateCubicSpline<Scalar>( splineCoeff_, i, t_values[n], outArg(*x) );
DS.time = t_values[n];
DS.x = x;
DS.accuracy = ST::zero();
data_out->push_back(DS);
}
// Move to the next user time point to consider!
n++;
if (n == as<int>(t_values.size())) {
// WE ARE ALL DONE! MOVE OUT!
return;
}
}
// Move on the the next interpolation time range
}
} // (*nodes_).size() == 1
}
开发者ID:00liujj,项目名称:trilinos,代码行数:93,代码来源:Rythmos_CubicSplineInterpolator_def.hpp
示例15: rcp
bool tBlockJacobiPreconditionerFactory::test_initializePrec(int verbosity,std::ostream & os)
{
using Thyra::zero;
bool status = false;
bool allPassed = true;
std::string constrType[3] = {
std::string("Static"),
std::string("2x2 Static Strategy"),
std::string("3x3 Static Strategy") };
// three by three bloock diagonal
std::vector<RCP<const Thyra::LinearOpBase<double> > > invD;
invD.push_back(invF_); invD.push_back(invC_); invD.push_back(invF_);
// allocate new linear operator
const RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blkOp
= Thyra::defaultBlockedLinearOp<double>();
blkOp->beginBlockFill(3,3);
blkOp->setBlock(0,0,F_); blkOp->setBlock(0,1,Bt_);
blkOp->setBlock(1,0,B_); blkOp->setBlock(1,1,C_); blkOp->setBlock(1,2,B_);
blkOp->setBlock(2,1,Bt_); blkOp->setBlock(2,2,F_);
blkOp->endBlockFill();
const RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > invBlkOp
= Thyra::defaultBlockedLinearOp<double>();
invBlkOp->beginBlockFill(3,3);
invBlkOp->setBlock(0,0,invF_);
invBlkOp->setBlock(1,1,invC_);
invBlkOp->setBlock(2,2,invF_);
invBlkOp->endBlockFill();
// build factory array
RCP<JacobiPreconditionerFactory> fact_array[3]
= { rcp(new JacobiPreconditionerFactory(invF_,invC_)),
rcp(new JacobiPreconditionerFactory(rcp(new StaticInvDiagStrategy(invF_,invC_)))),
rcp(new JacobiPreconditionerFactory(rcp(new StaticInvDiagStrategy(invD)))) };
RCP<const Thyra::LinearOpBase<double> > A[3] = {
block2x2(F_,Bt_,B_,C_),
block2x2(F_,Bt_,B_,C_),
blkOp };
// this is what the factory should build
RCP<const Thyra::LinearOpBase<double> > invA[3] = {
block2x2(invF_,zero(Bt_->range(),Bt_->domain()),zero(B_->range(),B_->domain()),invC_),
block2x2(invF_,zero(Bt_->range(),Bt_->domain()),zero(B_->range(),B_->domain()),invC_),
invBlkOp };
// test both constructors
for(int i=0;i<3;i++) {
RCP<const Thyra::LinearOpBase<double> > op;
RCP<Thyra::PreconditionerFactoryBase<double> > fact = fact_array[i];
RCP<Thyra::PreconditionerBase<double> > prec = fact->createPrec();
// initialize the preconditioner
fact->initializePrec(Thyra::defaultLinearOpSource(A[i]), &*prec);
op = prec->getRightPrecOp();
TEST_EQUALITY(op,Teuchos::null,
std::endl << " tBlockJacobiPreconditionerFactory::test_initializePrec "
<< "using \"" << constrType[i] << "\" constructor " << toString(status)
<< ": Preconditioner \"getRightPrecOp\" is not null (it should be!)");
op = prec->getLeftPrecOp();
TEST_EQUALITY(op,Teuchos::null,
std::endl << " tBlockJacobiPreconditionerFactory::test_initializePrec "
<< "using \"" << constrType[i] << "\" constructor " << toString(status)
<< ": Preconditioner \"getLeftPrecOp\" is not null (it should be!)");
op = prec->getUnspecifiedPrecOp();
TEST_NOT_EQUAL(op,Teuchos::null,
std::endl << " tBlockJacobiPreconditionerFactory::test_initializePrec "
<< "using \"" << constrType[i] << "\" constructor " << toString(status)
<< ": Preconditioner \"getUnspecifiedPrecOp\" is null (it should not be!)");
LinearOpTester<double> tester;
tester.show_all_tests(true);
std::stringstream ss;
Teuchos::FancyOStream fos(rcpFromRef(ss)," |||");
const bool result = tester.compare( *invA[i], *op, &fos );
TEST_ASSERT(result,
std::endl << " tBlockJacobiPreconditionerFactory::test_initializePrec "
<< ": Comparing factory generated operator to correct operator");
if(not result || verbosity>=10)
os << ss.str();
}
return allPassed;
}
开发者ID:haripandey,项目名称:trilinos,代码行数:92,代码来源:tBlockJacobiPreconditionerFactory.cpp
示例16: exampleImplicitlyComposedLinearOperators
int exampleImplicitlyComposedLinearOperators(
const int n0,
const int n1,
const int n2,
Teuchos::FancyOStream &out,
const Teuchos::EVerbosityLevel verbLevel,
typename Teuchos::ScalarTraits<Scalar>::magnitudeType errorTol,
const bool testAdjoint
)
{
// Using and other declarations
typedef Teuchos::ScalarTraits<Scalar> ST;
using Teuchos::as;
using Teuchos::RCP;
using Teuchos::OSTab;
using Thyra::VectorSpaceBase;
using Thyra::VectorBase;
using Thyra::MultiVectorBase;
using Thyra::LinearOpBase;
using Thyra::defaultSpmdVectorSpace;
using Thyra::randomize;
using Thyra::identity;
using Thyra::diagonal;
using Thyra::multiply;
using Thyra::add;
using Thyra::subtract;
using Thyra::scale;
using Thyra::adjoint;
using Thyra::block1x2;
using Thyra::block2x2;
using Thyra::block2x2;
out << "\n***"
<< "\n*** Demonstrating building linear operators for scalar type "
<< ST::name()
<< "\n***\n";
OSTab tab(out);
//
// A) Set up the basic objects and other inputs to build the implicitly
// composed linear operators.
//
// Create serial vector spaces in this case
const RCP<const VectorSpaceBase<Scalar> >
space0 = defaultSpmdVectorSpace<Scalar>(n0),
space1 = defaultSpmdVectorSpace<Scalar>(n1),
space2 = defaultSpmdVectorSpace<Scalar>(n2);
// Create the component linear operators first as multi-vectors
const RCP<MultiVectorBase<Scalar> >
mvA = createMembers(space2, n0, "A"),
mvB = createMembers(space0, n2, "B"),
mvC = createMembers(space0, n0, "C"),
mvE = createMembers(space0, n1, "E"),
mvF = createMembers(space0, n1, "F"),
mvJ = createMembers(space2, n1, "J"),
mvK = createMembers(space1, n2, "K"),
mvL = createMembers(space2, n1, "L"),
mvN = createMembers(space0, n1, "N"),
mvP = createMembers(space2, n1, "P"),
mvQ = createMembers(space0, n2, "Q");
// Create the vector diagonal for D
const RCP<VectorBase<Scalar> > d = createMember(space2);
// Get the constants
const Scalar
one = 1.0,
beta = 2.0,
gamma = 3.0,
eta = 4.0;
// Randomize the values in the Multi-Vector
randomize( -one, +one, mvA.ptr() );
randomize( -one, +one, mvB.ptr() );
randomize( -one, +one, mvC.ptr() );
randomize( -one, +one, d.ptr() );
randomize( -one, +one, mvE.ptr() );
randomize( -one, +one, mvF.ptr() );
randomize( -one, +one, mvJ.ptr() );
randomize( -one, +one, mvK.ptr() );
randomize( -one, +one, mvL.ptr() );
randomize( -one, +one, mvN.ptr() );
randomize( -one, +one, mvP.ptr() );
randomize( -one, +one, mvQ.ptr() );
// Get the linear operator forms of the basic component linear operators
const RCP<const LinearOpBase<Scalar> >
A = mvA,
B = mvB,
C = mvC,
E = mvE,
F = mvF,
J = mvJ,
K = mvK,
L = mvL,
N = mvN,
//.........这里部分代码省略.........
开发者ID:haripandey,项目名称:trilinos,代码行数:101,代码来源:exampleImplicitlyComposedLinearOperators.cpp
示例17: computeCubicSplineCoeff
void computeCubicSplineCoeff(
const typename DataStore<Scalar>::DataStoreVector_t & data,
const Ptr<CubicSplineCoeff<Scalar> > & coeffPtr
)
{
using Teuchos::outArg;
typedef Teuchos::ScalarTraits<Scalar> ST;
using Thyra::createMember;
TEUCHOS_TEST_FOR_EXCEPTION(
(data.size() < 2), std::logic_error,
"Error! A minimum of two data points is required for this cubic spline."
);
// time data in the DataStoreVector should be unique and sorted
Array<Scalar> t;
Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
dataStoreVectorToVector<Scalar>( data, &t, &x_vec, &xdot_vec, NULL );
#ifdef HAVE_RYTHMOS_DEBUG
assertTimePointsAreSorted<Scalar>( t );
#endif // HAVE_RYTHMOS_DEBUG
// 11/18/08 tscoffe: Question: Should I erase everything in coeffPtr or
// re-use what I can? For now, I'll erase and create new each time.
CubicSplineCoeff<Scalar>& coeff = *coeffPtr;
// If there are only two points, then we do something special and just create
// a linear polynomial between the points
|
请发表评论