本文整理汇总了C++中Teuchos类的典型用法代码示例。如果您正苦于以下问题:C++ Teuchos类的具体用法?C++ Teuchos怎么用?C++ Teuchos使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Teuchos类的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc, char *argv[])
{
using std::endl;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ParameterList;
using Teuchos::CommandLineProcessor;
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);
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 );
bool dumpFinalSolutions = false;
clp.setOption(
"dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
"Determine if the final solutions are dumpped or not." );
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 );
//
// Get the base parameter list that all other parameter lists will be read
// from.
//
RCP<ParameterList>
paramList = Teuchos::parameterList();
//
// Create the underlying EpetraExt::ModelEvaluator
//
RCP<LorenzModel>
lorenzModel = rcp(new LorenzModel( epetra_comm, *paramList ));
//
// Create the Thyra-wrapped ModelEvaluator
//
RCP<Thyra::ModelEvaluator<double> >
thyraLorenzModel = Thyra::epetraModelEvaluator(lorenzModel,Teuchos::null);
//
// Create the Rythmos GAASP ErrorEstimator
//
RCP<Rythmos::GAASPErrorEstimator> gaaspEE = rcp(new Rythmos::GAASPErrorEstimator);
gaaspEE->setModel(thyraLorenzModel);
//gaaspEE->setQuantityOfInterest( AVERAGE_ERROR_QTY ); // Not passed through yet.
Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
pl->sublist("GAASP Interface Parameters").set<double>("eTime",1.0);
pl->sublist("GAASP Interface Parameters").set<double>("timeStep",0.1);
//pl->sublist("GAASP Interface Parameters").set<double>("timeStep",0.5);
gaaspEE->setParameterList(pl);
//RCP<const Rythmos::ErrorEstimateBase<double> > error = gaaspEE->getErrorEstimate();
//double uTOL = 1.0e-8;
double uTOL = 1.0e-2;
RCP<const Rythmos::ErrorEstimateBase<double> > error = gaaspEE->controlGlobalError(uTOL);
double err = error->getTotalError();
out->precision(15);
*out << "err = " << err << std::endl;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
//.........这里部分代码省略.........
示例2: TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Belos_Hypre, Laplace2D){
const double tol = 1E-7;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ParameterList;
typedef Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> LinearProblem;
//
// Create Laplace2D
//
#ifdef HAVE_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm();
#endif
Teuchos::ParameterList GaleriList;
int nx = 10 * Comm.NumProc();
int ny = 10 * Comm.NumProc();
GaleriList.set("nx", nx);
GaleriList.set("ny", ny);
Epetra_Map Map(nx*ny,0,Comm);
RCP<Epetra_CrsMatrix> Crs_Matrix = rcp(Galeri::CreateCrsMatrix("Laplace2D", &Map, GaleriList));
int NumProc = Crs_Matrix->Comm().NumProc();
//
// Create the hypre preconditioner
//
RCP<Ifpack_Hypre> preconditioner = rcp(new Ifpack_Hypre(Crs_Matrix.get()));
TEST_EQUALITY(preconditioner->Initialize(),0);
TEST_EQUALITY(preconditioner->SetParameter(Preconditioner, ParaSails),0); // Use a Euclid Preconditioner (but not really used)
TEST_EQUALITY(preconditioner->SetParameter(Preconditioner),0); // Solve the problem
TEST_EQUALITY(preconditioner->Compute(),0);
//
// Create the solution vector and rhs
//
int numVec = 1;
RCP<Epetra_MultiVector> X = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorDomainMap(), numVec));
RCP<Epetra_MultiVector> KnownX = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorDomainMap(), numVec));
KnownX->Random();
RCP<Epetra_MultiVector> B = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorRangeMap(), numVec));
Crs_Matrix->Apply(*KnownX, *B);
//
// Test the EpetraExt wrapper
// amk November 24, 2015: Should we deprecate this?
//
// RCP<ParameterList> pl = rcp(new ParameterList());
// TEST_EQUALITY(X->PutScalar(0.0),0);
// HYPRE_IJMatrix hypre_mat = preconditioner->HypreMatrix();
// RCP<EpetraExt_HypreIJMatrix> Hyp_Matrix = rcp(new EpetraExt_HypreIJMatrix(hypre_mat));
// TEST_EQUALITY(Hyp_Matrix->SetParameter(Preconditioner, ParaSails),0);
// TEST_EQUALITY(Hyp_Matrix->SetParameter(Preconditioner),0);
// TEST_EQUALITY(EquivalentMatrices(*Hyp_Matrix, *Crs_Matrix, tol), true);
// RCP<LinearProblem> problem1 = rcp(new LinearProblem(Crs_Matrix,X,B));
// problem1->setLeftPrec(Hyp_Matrix);
// TEST_EQUALITY(problem1->setProblem(),true);
// Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator> solMgr1(problem1,pl);
// Belos::ReturnType rv1 = solMgr1.solve(); // TEST_EQUALITY(solMgr2.solve(),Belos::Converged);
// TEST_EQUALITY(rv1,Belos::Converged);
// TEST_EQUALITY(EquivalentVectors(*X, *KnownX, tol*10*pow(10.0,NumProc)), true);
//
// Test the Ifpack hypre interface
//
RCP<ParameterList> pl2 = rcp(new ParameterList());
RCP<Epetra_Operator> invOp = rcp(new Epetra_InvOperator(preconditioner.get()));
TEST_EQUALITY(X->PutScalar(0.0),0);
RCP<LinearProblem> problem2 = rcp(new LinearProblem(Crs_Matrix,X,B));
problem2->setLeftPrec(invOp);
TEST_EQUALITY(problem2->setProblem(),true);
Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator> solMgr2(problem2,pl2);
Belos::ReturnType rv2 = solMgr2.solve(); // TEST_EQUALITY(solMgr2.solve(),Belos::Converged);
TEST_EQUALITY(rv2,Belos::Converged);
TEST_EQUALITY(EquivalentVectors(*X, *KnownX, tol*10*pow(10.0,NumProc)), true);
}
示例3: TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(tIterativePreconditionerFactory, parameter_list_init)
{
// build global (or serial communicator)
#ifdef HAVE_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
Teko::LinearOp A = build2x2(Comm,1,2,3,4);
Thyra::LinearOpTester<double> tester;
tester.show_all_tests(true);
{
RCP<Teuchos::ParameterList> pl = buildLibPL(4,"Amesos");
RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());
RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));
try {
precFact->initializeFromParameterList(*pl);
out << "Passed correct parameter list" << std::endl;
Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
}
catch(...) {
success = false;
out << "Failed correct parameter list" << std::endl;
}
}
{
Teuchos::ParameterList pl;
pl.set("Preconditioner Type","Amesos");
RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());
RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));
try {
precFact->initializeFromParameterList(pl);
out << "Passed iteration count" << std::endl;
Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
}
catch(...) {
out << "Failed iteration count" << std::endl;
}
}
{
Teuchos::ParameterList pl;
pl.set("Iteration Count",4);
pl.set("Precondiioner Type","Amesos");
RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());
try {
precFact->initializeFromParameterList(pl);
success = false;
out << "Failed preconditioner type" << std::endl;
// these should not be executed
RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));
Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
}
catch(const std::exception & exp) {
out << "Passed preconditioner type" << std::endl;
}
}
}
示例4: main
int main(int argc, char *argv[]) {
Teuchos::GlobalMPISession mpiSession(&argc,&argv);
typedef double Scalar;
typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
typedef Tpetra::Map<>::local_ordinal_type LO;
typedef Tpetra::Map<>::global_ordinal_type GO;
typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;
typedef Tpetra::CrsMatrix<Scalar,LO,GO> MAT;
typedef Tpetra::MultiVector<Scalar,LO,GO> MV;
using Tpetra::global_size_t;
using Teuchos::tuple;
using Teuchos::RCP;
using Teuchos::rcp;
Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
Teuchos::RCP<const Teuchos::Comm<int> > comm = platform.getComm();
size_t myRank = comm->getRank();
RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
*fos << Amesos2::version() << std::endl << std::endl;
bool printMatrix = false;
bool printSolution = false;
bool printTiming = false;
bool printResidual = false;
bool printLUStats = false;
bool verbose = false;
std::string solver_name = "SuperLU";
std::string filedir = "../test/matrices/";
std::string filename = "arc130.mtx";
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("filedir",&filedir,"Directory where matrix-market files are located");
cmdp.setOption("filename",&filename,"Filename for Matrix-Market test matrix.");
cmdp.setOption("print-matrix","no-print-matrix",&printMatrix,"Print the full matrix after reading it.");
cmdp.setOption("print-solution","no-print-solution",&printSolution,"Print solution vector after solve.");
cmdp.setOption("print-timing","no-print-timing",&printTiming,"Print solver timing statistics");
cmdp.setOption("print-residual","no-print-residual",&printResidual,"Print solution residual");
cmdp.setOption("print-lu-stats","no-print-lu-stats",&printLUStats,"Print nnz in L and U factors");
cmdp.setOption("solver", &solver_name, "Which TPL solver library to use.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
// Before we do anything, check that the solver is enabled
if( !Amesos2::query(solver_name) ){
std::cerr << solver_name << " not enabled. Exiting..." << std::endl;
return EXIT_SUCCESS; // Otherwise CTest will pick it up as
// failure, which it isn't really
}
const size_t numVectors = 1;
// create a Map
global_size_t nrows = 6;
RCP<Tpetra::Map<LO,GO> > map
= rcp( new Tpetra::Map<LO,GO>(nrows,0,comm) );
std::string mat_pathname = filedir + filename;
RCP<MAT> A = Tpetra::MatrixMarket::Reader<MAT>::readSparseFile(mat_pathname,comm);
if( printMatrix ){
A->describe(*fos, Teuchos::VERB_EXTREME);
}
else if( verbose && myRank==0 ){
*fos << std::endl << A->description() << std::endl << std::endl;
}
// get the maps
RCP<const Tpetra::Map<LO,GO> > dmnmap = A->getDomainMap();
RCP<const Tpetra::Map<LO,GO> > rngmap = A->getRangeMap();
// Create random X
RCP<MV> Xhat = rcp( new MV(dmnmap,numVectors) );
RCP<MV> X = rcp( new MV(dmnmap,numVectors) );
X->setObjectLabel("X");
Xhat->setObjectLabel("Xhat");
X->randomize();
RCP<MV> B = rcp(new MV(rngmap,numVectors));
A->apply(*X, *B);
// Constructor from Factory
RCP<Amesos2::Solver<MAT,MV> > solver;
try{
solver = Amesos2::create<MAT,MV>(solver_name, A, Xhat, B);
} catch (std::invalid_argument e){
*fos << e.what() << std::endl;
return 0;
}
#ifdef SHYLU_NODEBASKER
if( Amesos2::query("Basker") ) {
Teuchos::ParameterList amesos2_params("Amesos2");
amesos2_params.sublist(solver_name).set("num_threads", 1, "Number of threads");
solver->setParameters( Teuchos::rcpFromRef(amesos2_params) );
//.........这里部分代码省略.........
示例5: main
int main(int argc, char *argv[]) {
//
int MyPID = 0;
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
MyPID = Comm.MyPID();
#else
Epetra_SerialComm Comm;
#endif
//
typedef double ST;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
typedef Belos::MultiVecTraits<ST,MV> MVT;
typedef Belos::OperatorTraits<ST,MV,OP> OPT;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
bool success = false;
bool verbose = false;
try {
bool proc_verbose = false;
int frequency = -1; // frequency of status test output.
int blocksize = 1; // blocksize
int numrhs = 1; // number of right-hand sides to solve for
int maxrestarts = 15; // maximum number of restarts allowed
int maxiters = -1; // maximum number of iterations allowed per linear system
int maxsubspace = 25; // maximum number of blocks the solver can use for the subspace
std::string filename("orsirr1.hb");
MT tol = 1.0e-5; // relative residual tolerance
// Specify whether to use RHS as initial guess. If false, use zero.
bool useRHSAsInitialGuess = false;
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("use-rhs","use-zero",&useRHSAsInitialGuess,"Use RHS as initial guess.");
cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
cmdp.setOption("filename",&filename,"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return EXIT_FAILURE;
}
if (!verbose)
frequency = -1; // reset frequency if test is not verbose
//
// *************Get the problem*********************
//
RCP<Epetra_Map> Map;
RCP<Epetra_CrsMatrix> A;
RCP<Epetra_MultiVector> B, X;
RCP<Epetra_Vector> vecB, vecX;
EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB);
A->OptimizeStorage();
proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
// Check to see if the number of right-hand sides is the same as requested.
if (numrhs>1) {
X = rcp( new Epetra_MultiVector( *Map, numrhs ) );
B = rcp( new Epetra_MultiVector( *Map, numrhs ) );
X->Random();
OPT::Apply( *A, *X, *B );
X->PutScalar( 0.0 );
}
else {
X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
}
// If requested, use a copy of B as initial guess
if (useRHSAsInitialGuess)
{
X->Update(1.0, *B, 0.0);
}
//
// ************Construct preconditioner*************
//
ParameterList ifpackList;
// allocates an IFPACK factory. No data is associated
// to this object (only method Create()).
Ifpack Factory;
// create the preconditioner. For valid PrecType values,
// please check the documentation
std::string PrecType = "ILU"; // incomplete LU
//.........这里部分代码省略.........
示例6: main
int main(int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>
typedef Tpetra::Vector<SC,LO,GO,NO> TVEC;
typedef Tpetra::MultiVector<SC,LO,GO,NO> TMV;
typedef Tpetra::CrsMatrix<SC,LO,GO,NO,LMO> TCRS;
typedef Xpetra::CrsMatrix<SC,LO,GO,NO,LMO> XCRS;
typedef Xpetra::TpetraCrsMatrix<SC,LO,GO,NO,LMO> XTCRS;
typedef Xpetra::Matrix<SC,LO,GO,NO,LMO> XMAT;
typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO,LMO> XWRAP;
typedef Belos::OperatorT<TMV> TOP;
typedef Belos::OperatorTraits<SC,TMV,TOP> TOPT;
typedef Belos::MultiVecTraits<SC,TMV> TMVT;
typedef Belos::LinearProblem<SC,TMV,TOP> TProblem;
typedef Belos::SolverManager<SC,TMV,TOP> TBelosSolver;
typedef Belos::BlockGmresSolMgr<SC,TMV,TOP> TBelosGMRES;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::TimeMonitor;
Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
Teuchos::CommandLineProcessor clp(false);
GO nx,ny,nz;
nx=100; ny=100; nz=100;
double stretchx, stretchy, stretchz, h, delta;
stretchx=1.0; stretchy=1.0; stretchz=1.0;
h=0.01; delta=2.0;
int PMLXL, PMLXR, PMLYL, PMLYR, PMLZL, PMLZR;
PMLXL=10; PMLXR=10; PMLYL=10; PMLYR=10; PMLZL=10; PMLZR=10;
double omega, shift;
omega=20.0*M_PI;
shift=0.5;
Galeri::Xpetra::Parameters<GO> matrixParameters(clp, nx, ny, nz, "Helmholtz1D", 0, stretchx, stretchy, stretchz,
h, delta, PMLXL, PMLXR, PMLYL, PMLYR, PMLZL, PMLZR, omega, shift);
Xpetra::Parameters xpetraParameters(clp);
RCP<TimeMonitor> globalTimeMonitor = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time")));
RCP<TimeMonitor> tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));
Teuchos::ParameterList pl = matrixParameters.GetParameterList();
RCP<MultiVector> coordinates;
Teuchos::ParameterList galeriList;
galeriList.set("nx", pl.get("nx", nx));
galeriList.set("ny", pl.get("ny", ny));
galeriList.set("nz", pl.get("nz", nz));
RCP<const Map> map;
if (matrixParameters.GetMatrixType() == "Helmholtz1D") {
map = MapFactory::Build(xpetraParameters.GetLib(), matrixParameters.GetNumGlobalElements(), 0, comm);
coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("1D", map, matrixParameters.GetParameterList());
}
else if (matrixParameters.GetMatrixType() == "Helmholtz2D") {
map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian2D", comm, galeriList);
coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("2D", map, matrixParameters.GetParameterList());
}
else if (matrixParameters.GetMatrixType() == "Helmholtz3D") {
map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian3D", comm, galeriList);
coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("3D", map, matrixParameters.GetParameterList());
}
RCP<const Tpetra::Map<LO, GO, NO> > tmap = Xpetra::toTpetra(map);
Teuchos::ParameterList matrixParams = matrixParameters.GetParameterList();
// Build problem
RCP<Galeri::Xpetra::Problem_Helmholtz<Map,CrsMatrixWrap,MultiVector> > Pr =
Galeri::Xpetra::BuildProblem_Helmholtz<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>(matrixParameters.GetMatrixType(), map, matrixParams);
RCP<Matrix> A = Pr->BuildMatrix();
RCP<MultiVector> nullspace = MultiVectorFactory::Build(map,1);
nullspace->putScalar( (SC) 1.0);
comm->barrier();
tm = Teuchos::null;
// Construct a multigrid preconditioner
tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 2 - MueLu Setup")));
// Multigrid Hierarchy
RCP<Hierarchy> H = rcp(new Hierarchy(A));
H->GetLevel(0)->Set("Nullspace",nullspace);
FactoryManager Manager;
H->Setup(Manager, 0, 5);
//H->Write(-1,-1);
tm = Teuchos::null;
// Solve Ax = b
tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 3 - LHS and RHS initialization")));
RCP<TVEC> X = Tpetra::createVector<SC,LO,GO,NO>(tmap);
RCP<TVEC> B = Tpetra::createVector<SC,LO,GO,NO>(tmap);
X->putScalar((SC) 0.0);
B->putScalar((SC) 0.0);
if(comm->getRank()==0) {
//.........这里部分代码省略.........
示例7: main
int main(int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ArrayRCP;
using Teuchos::TimeMonitor;
using Teuchos::ParameterList;
// =========================================================================
// MPI initialization using Teuchos
// =========================================================================
Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
// =========================================================================
// Convenient definitions
// =========================================================================
typedef Teuchos::ScalarTraits<SC> STS;
SC zero = STS::zero(), one = STS::one();
// =========================================================================
// Parameters initialization
// =========================================================================
Teuchos::CommandLineProcessor clp(false);
GO nx = 100, ny = 100, nz = 100;
Galeri::Xpetra::Parameters<GO> galeriParameters(clp, nx, ny, nz, "Laplace2D"); // manage parameters of the test case
Xpetra::Parameters xpetraParameters(clp); // manage parameters of Xpetra
std::string xmlFileName = "scalingTest.xml"; clp.setOption("xml", &xmlFileName, "read parameters from a file [default = 'scalingTest.xml']");
bool printTimings = true; clp.setOption("timings", "notimings", &printTimings, "print timings to screen");
int writeMatricesOPT = -2; clp.setOption("write", &writeMatricesOPT, "write matrices to file (-1 means all; i>=0 means level i)");
std::string dsolveType = "cg", solveType; clp.setOption("solver", &dsolveType, "solve type: (none | cg | gmres | standalone)");
double dtol = 1e-12, tol; clp.setOption("tol", &dtol, "solver convergence tolerance");
std::string mapFile; clp.setOption("map", &mapFile, "map data file");
std::string matrixFile; clp.setOption("matrix", &matrixFile, "matrix data file");
std::string coordFile; clp.setOption("coords", &coordFile, "coordinates data file");
int numRebuilds = 0; clp.setOption("rebuild", &numRebuilds, "#times to rebuild hierarchy");
int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations");
bool scaleResidualHistory = true; clp.setOption("scale", "noscale", &scaleResidualHistory, "scaled Krylov residual history");
switch (clp.parse(argc, argv)) {
case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED: return EXIT_SUCCESS;
case Teuchos::CommandLineProcessor::PARSE_ERROR:
case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE;
case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL: break;
}
Xpetra::UnderlyingLib lib = xpetraParameters.GetLib();
ParameterList paramList;
Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<ParameterList>(¶mList), *comm);
bool isDriver = paramList.isSublist("Run1");
if (isDriver) {
// update galeriParameters with the values from the XML file
ParameterList& realParams = galeriParameters.GetParameterList();
for (ParameterList::ConstIterator it = realParams.begin(); it != realParams.end(); it++) {
const std::string& name = realParams.name(it);
if (paramList.isParameter(name))
realParams.setEntry(name, paramList.getEntry(name));
}
}
// Retrieve matrix parameters (they may have been changed on the command line)
// [for instance, if we changed matrix type from 2D to 3D we need to update nz]
ParameterList galeriList = galeriParameters.GetParameterList();
// =========================================================================
// Problem construction
// =========================================================================
std::ostringstream galeriStream;
comm->barrier();
RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time")));
RCP<TimeMonitor> tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));
RCP<Matrix> A;
RCP<const Map> map;
RCP<MultiVector> coordinates;
RCP<MultiVector> nullspace;
if (matrixFile.empty()) {
galeriStream << "========================================================\n" << xpetraParameters << galeriParameters;
// Galeri will attempt to create a square-as-possible distribution of subdomains di, e.g.,
// d1 d2 d3
// d4 d5 d6
// d7 d8 d9
// d10 d11 d12
// A perfect distribution is only possible when the #processors is a perfect square.
// This *will* result in "strip" distribution if the #processors is a prime number or if the factors are very different in
// size. For example, np=14 will give a 7-by-2 distribution.
// If you don't want Galeri to do this, specify mx or my on the galeriList.
std::string matrixType = galeriParameters.GetMatrixType();
// Create map and coordinates
// In the future, we hope to be able to first create a Galeri problem, and then request map and coordinates from it
// At the moment, however, things are fragile as we hope that the Problem uses same map and coordinates inside
if (matrixType == "Laplace1D") {
//.........这里部分代码省略.........
示例8: valuestemp
void
SupportGraph<MatrixType>::findSupport ()
{
typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
global_ordinal_type, node_type> crs_matrix_type;
typedef Tpetra::Vector<scalar_type, local_ordinal_type,
global_ordinal_type, node_type> vec_type;
typedef std::pair<int, int> E;
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
boost::no_property,
boost::property<boost::edge_weight_t,
magnitude_type> > graph_type;
typedef typename boost::graph_traits<graph_type>::edge_descriptor edge_type;
typedef typename boost::graph_traits<graph_type>::vertex_descriptor
vertex_type;
const scalar_type zero = STS::zero();
const scalar_type one = STS::one();
//Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cout));
size_t num_verts = A_local_->getNodeNumRows();
size_t num_edges
= (A_local_->getNodeNumEntries() - A_local_->getNodeNumDiags())/2;
// Create data structures for the BGL code
// and temp data structures for extraction
E *edge_array = new E[num_edges];
magnitude_type *weights = new magnitude_type[num_edges];
size_t num_entries;
size_t max_num_entries = A_local_->getNodeMaxNumRowEntries();
std::vector<scalar_type> valuestemp (max_num_entries);
std::vector<local_ordinal_type> indicestemp (max_num_entries);
std::vector<magnitude_type> diagonal (num_verts);
Tpetra::ArrayView<scalar_type> values (valuestemp);
Tpetra::ArrayView<local_ordinal_type> indices (indicestemp);
// Extract from the tpetra matrix keeping only one edge per pair
// (assume symmetric)
size_t offDiagCount = 0;
for (size_t row = 0; row < num_verts; ++row) {
A_local_->getLocalRowCopy (row, indices, values, num_entries);
for (size_t colIndex = 0; colIndex < num_entries; ++colIndex) {
if(row == Teuchos::as<size_t>(indices[colIndex])) {
diagonal[row] = values[colIndex];
}
if((row < Teuchos::as<size_t>(indices[colIndex]))
&& (values[colIndex] < zero)) {
edge_array[offDiagCount] = E(row, indices[colIndex]);
weights[offDiagCount] = values[colIndex];
if (Randomize_) {
// Add small random pertubation.
weights[offDiagCount] *= one +
STS::magnitude(STS::rmin() * STS::random());
}
offDiagCount++;
}
}
}
// Create BGL graph
graph_type g(edge_array, edge_array + num_edges, weights, num_verts);
typedef typename boost::property_map
<graph_type, boost::edge_weight_t>::type type;
type weight = get (boost::edge_weight, g);
std::vector<edge_type> spanning_tree;
// Run Kruskal, actually maximal weight ST since edges are negative
boost::kruskal_minimum_spanning_tree(g, std::back_inserter (spanning_tree));
// Create array to store the exact number of non-zeros per row
Teuchos::ArrayRCP<size_t> NumNz (num_verts, 1);
typedef typename std::vector<edge_type>::iterator edge_iterator_type;
// Find the degree of all the vertices
for (edge_iterator_type ei = spanning_tree.begin(); ei != spanning_tree.end();
++ei) {
local_ordinal_type localsource = source(*ei,g);
local_ordinal_type localtarget = target(*ei,g);
// We only want upper triangular entries, might need to swap
if (localsource > localtarget) {
localsource = target(*ei, g);
localtarget = source(*ei, g);
}
NumNz[localsource] += 1;
}
// Create an stl vector of stl vectors to hold indices and values
std::vector<std::vector<local_ordinal_type> > Indices (num_verts);
//.........这里部分代码省略.........
示例9: apply
void
SupportGraph<MatrixType>::
apply (const Tpetra::MultiVector<scalar_type,
local_ordinal_type,
global_ordinal_type,
node_type>& X,
Tpetra::MultiVector<scalar_type,
local_ordinal_type,
global_ordinal_type,
node_type>& Y,
Teuchos::ETransp mode,
scalar_type alpha,
scalar_type beta) const
{
using Teuchos::FancyOStream;
using Teuchos::getFancyOStream;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::rcpFromRef;
using Teuchos::Time;
using Teuchos::TimeMonitor;
typedef scalar_type DomainScalar;
typedef scalar_type RangeScalar;
typedef Tpetra::MultiVector<DomainScalar, local_ordinal_type,
global_ordinal_type, node_type> MV;
RCP<FancyOStream> out = getFancyOStream(rcpFromRef(std::cout));
// Create a timer for this method, if it doesn't exist already.
// TimeMonitor::getNewCounter registers the timer, so that
// TimeMonitor's class methods like summarize() will report the
// total time spent in successful calls to this method.
const std::string timerName ("Ifpack2::SupportGraph::apply");
RCP<Time> timer = TimeMonitor::lookupCounter(timerName);
if (timer.is_null()) {
timer = TimeMonitor::getNewCounter(timerName);
}
{ // Start timing here.
Teuchos::TimeMonitor timeMon (*timer);
TEUCHOS_TEST_FOR_EXCEPTION(
! isComputed(), std::runtime_error,
"Ifpack2::SupportGraph::apply: You must call compute() to compute the "
"incomplete factorization, before calling apply().");
TEUCHOS_TEST_FOR_EXCEPTION(
X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
"Ifpack2::SupportGraph::apply: X and Y must have the same number of "
"columns. X has " << X.getNumVectors() << " columns, but Y has "
<< Y.getNumVectors() << " columns.");
TEUCHOS_TEST_FOR_EXCEPTION(
beta != STS::zero(), std::logic_error,
"Ifpack2::SupportGraph::apply: This method does not currently work when "
"beta != 0.");
// If X and Y are pointing to the same memory location,
// we need to create an auxiliary vector, Xcopy
RCP<const MV> Xcopy;
if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
Xcopy = rcp (new MV(X));
}
else {
Xcopy = rcpFromRef(X);
}
if (alpha != STS::one()) {
Y.scale(alpha);
}
RCP<MV> Ycopy = rcpFromRef(Y);
solver_->setB(Xcopy);
solver_->setX(Ycopy);
solver_->solve ();
} // Stop timing here.
++NumApply_;
// timer->totalElapsedTime() returns the total time over all timer
// calls. Thus, we use = instead of +=.
ApplyTime_ = timer->totalElapsedTime();
}
示例10: main
int main(int argc, char *argv[]) {
//
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init(&argc,&argv);
#endif
//
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
bool success = false;
bool verbose = false;
try {
//
// Get test parameters from command-line processor
//
bool proc_verbose = false;
int frequency = -1; // frequency of status test output.
std::string filename("gr_30_30.hb"); // default input filename
double tol = 1.0e-10; // relative residual tolerance
int numBlocks = 30; // maximum number of blocks the solver can use for the Krylov subspace
int recycleBlocks = 3; // maximum number of blocks the solver can use for the recycle space
int numrhs = 1; // number of right-hand sides to solve for
int maxiters = -1; // maximum number of iterations allowed per linear system
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("filename",&filename,"Filename for test matrix.");
cmdp.setOption("tol",&tol,"Relative residual tolerance used by the RCG solver.");
cmdp.setOption("max-subspace",&numBlocks,"Maximum number of vectors in search space (not including recycle space).");
cmdp.setOption("recycle",&recycleBlocks,"Number of vectors in recycle space.");
cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
if (!verbose)
frequency = -1; // reset frequency if test is not verbose
//
// Get the problem
//
int MyPID;
RCP<Epetra_CrsMatrix> A;
RCP<Epetra_MultiVector> B, X;
int return_val =Belos::createEpetraProblem(filename,NULL,&A,&B,&X,&MyPID);
if(return_val != 0) return return_val;
proc_verbose = ( verbose && (MyPID==0) );
//
// Solve using Belos
//
typedef double ST;
typedef Epetra_Operator OP;
typedef Epetra_MultiVector MV;
typedef Belos::OperatorTraits<ST,MV,OP> OPT;
typedef Belos::MultiVecTraits<ST,MV> MVT;
//
// *****Construct initial guess and right-hand sides *****
//
if (numrhs != 1) {
X = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
X->Random();
B = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
OPT::Apply( *A, *X, *B );
MVT::MvInit( *X, 0.0 );
}
else { // initialize exact solution to be vector of ones
MVT::MvInit( *X, 1.0 );
OPT::Apply( *A, *X, *B );
MVT::MvInit( *X, 0.0 );
}
//
// ********Other information used by block solver***********
// *****************(can be user specified)******************
//
const int NumGlobalElements = B->GlobalLength();
if (maxiters == -1)
maxiters = NumGlobalElements - 1; // maximum number of iterations to run
//
ParameterList belosList;
belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
belosList.set( "Num Blocks", numBlocks); // Maximum number of blocks in Krylov space
belosList.set( "Num Recycled Blocks", recycleBlocks ); // Number of vectors in recycle space
belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
if (verbose) {
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails );
if (frequency > 0)
belosList.set( "Output Frequency", frequency );
}
else
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
//
// Construct an unpreconditioned linear problem instance.
//
Belos::LinearProblem<double,MV,OP> problem( A, X, B );
bool set = problem.setProblem();
if (set == false) {
if (proc_verbose)
//.........这里部分代码省略.........
示例11: main
int main(int argc, char *argv[]) {
//Check number of arguments
if (argc < 4) {
std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
std::cout <<"Usage:\n\n";
std::cout <<" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n";
std::cout <<" where \n";
std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
exit(1);
}
// This little trick lets us print to std::cout only if
// a (dummy) command-line argument is provided.
int iprint = argc - 1;
Teuchos::RCP<std::ostream> outStream;
Teuchos::oblackholestream bhs; // outputs nothing
if (iprint > 2)
outStream = Teuchos::rcp(&std::cout, false);
else
outStream = Teuchos::rcp(&bhs, false);
// Save the format state of the original std::cout.
Teuchos::oblackholestream oldFormatState;
oldFormatState.copyfmt(std::cout);
*outStream \
<< "===============================================================================\n" \
<< "| |\n" \
<< "| Example: Apply Stiffness Matrix for |\n" \
<< "| Poisson Equation on Hexahedral Mesh |\n" \
<< "| |\n" \
<< "| Questions? Contact Pavel Bochev ([email protected]), |\n" \
<< "| Denis Ridzal ([email protected]), |\n" \
<< "| Kara Peterson ([email protected]). |\n" \
<< "| |\n" \
<< "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
<< "| Trilinos website: http://trilinos.sandia.gov |\n" \
<< "| |\n" \
<< "===============================================================================\n";
// ************************************ GET INPUTS **************************************
int deg = atoi(argv[1]); // polynomial degree to use
int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
// *********************************** CELL TOPOLOGY **********************************
// Get cell topology for base hexahedron
typedef shards::CellTopology CellTopology;
CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
// Get dimensions
int numNodesPerElem = hex_8.getNodeCount();
int spaceDim = hex_8.getDimension();
// *********************************** GENERATE MESH ************************************
*outStream << "Generating mesh ... \n\n";
*outStream << " NX" << " NY" << " NZ\n";
*outStream << std::setw(5) << NX <<
std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
// Print mesh information
int numElems = NX*NY*NZ;
int numNodes = (NX+1)*(NY+1)*(NZ+1);
*outStream << " Number of Elements: " << numElems << " \n";
*outStream << " Number of Nodes: " << numNodes << " \n\n";
// Cube
double leftX = 0.0, rightX = 1.0;
double leftY = 0.0, rightY = 1.0;
double leftZ = 0.0, rightZ = 1.0;
// Mesh spacing
double hx = (rightX-leftX)/((double)NX);
double hy = (rightY-leftY)/((double)NY);
double hz = (rightZ-leftZ)/((double)NZ);
// Get nodal coordinates
FieldContainer<double> nodeCoord(numNodes, spaceDim);
FieldContainer<int> nodeOnBoundary(numNodes);
int inode = 0;
for (int k=0; k<NZ+1; k++)
{
for (int j=0; j<NY+1; j++)
{
for (int i=0; i<NX+1; i++)
{
nodeCoord(inode,0) = leftX + (double)i*hx;
nodeCoord(inode,1) = leftY + (double)j*hy;
//.........这里部分代码省略.........
示例12: main
int main(int argc, char *argv[]) {
typedef double ST;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Tpetra::MultiVector<> MV;
typedef Tpetra::Operator<> OP;
typedef Belos::MultiVecTraits<ST,MV> MVT;
typedef Belos::OperatorTraits<ST,MV,OP> OPT;
typedef Tpetra::CrsMatrix<> CrsMatrix;
typedef Ifpack2::Preconditioner<> Prec;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
// ************************* Initialize MPI **************************
Teuchos::oblackholestream blackhole;
Teuchos::GlobalMPISession mpiSession (&argc, &argv, &blackhole);
// ************** Get the default communicator and node **************
RCP<const Teuchos::Comm<int> > comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
const int myRank = comm->getRank();
bool verbose = true;
bool success = true;
bool proc_verbose = false;
bool leftprec = true; // left preconditioning or right.
int frequency = -1; // frequency of status test output.
int numrhs = 1; // number of right-hand sides to solve for
int maxiters = -1; // maximum number of iterations allowed per linear system
std::string filename("cage4.mtx");
MT tol = 1.0e-5; // relative residual tolerance
// ***************** Read the command line arguments *****************
Teuchos::CommandLineProcessor cmdp(false,false);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("left-prec","right-prec",&leftprec,"Left preconditioning or right.");
cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
cmdp.setOption("filename",&filename,"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
if (!verbose)
frequency = -1; // reset frequency if test is not verbose
proc_verbose = verbose && (myRank==0); /* Only print on the zero processor */
// ************************* Get the problem *************************
RCP<CrsMatrix> A = Tpetra::MatrixMarket::Reader<CrsMatrix>::readSparseFile(filename,comm);
RCP<MV> B = rcp(new MV(A->getRowMap(),numrhs,false));
RCP<MV> X = rcp(new MV(A->getRowMap(),numrhs,false));
OPT::Apply(*A, *X, *B);
MVT::MvInit(*X);
MVT::MvInit(*B,1);
// ******************** Construct preconditioner *********************
Ifpack2::Factory factory;
RCP<Prec> M = factory.create("RELAXATION", A.getConst());
ParameterList ifpackParams;
ifpackParams.set("relaxation: type","Jacobi");
M->setParameters(ifpackParams);
M->initialize();
M->compute();
// ******* Create parameter list for the Belos solver manager ********
const int NumGlobalElements = MVT::GetGlobalLength(*B);
if (maxiters == -1)
maxiters = NumGlobalElements - 1; // maximum number of iterations to run
//
ParameterList belosList;
belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
if (numrhs > 1) {
belosList.set( "Show Maximum Residual Norm Only", true ); // Show only the maximum residual norm
}
if (verbose) {
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
Belos::TimingDetails + Belos::StatusTestDetails );
if (frequency > 0)
belosList.set( "Output Frequency", frequency );
}
else
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + Belos::FinalSummary );
// ************ Construct a preconditioned linear problem ************
RCP<Belos::LinearProblem<double,MV,OP> > problem
= rcp( new Belos::LinearProblem<double,MV,OP>( A, X, B ) );
if (leftprec) {
problem->setLeftPrec( M );
}
else {
problem->setRightPrec( M );
}
bool set = problem->setProblem();
if (set == false) {
if (proc_verbose)
std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
return -1;
//.........这里部分代码省略.........
示例13: sublist
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!
//.........这里部分代码省略.........
示例14: buildClosureModels
Teuchos::RCP< std::vector< Teuchos::RCP<PHX::Evaluator<panzer::Traits> > > >
user_app::MyModelFactory<EvalT>::
buildClosureModels(const std::string& model_id,
const Teuchos::ParameterList& models,
const panzer::FieldLayoutLibrary& fl,
const Teuchos::RCP<panzer::IntegrationRule>& ir,
const Teuchos::ParameterList& default_params,
const Teuchos::ParameterList& user_data,
const Teuchos::RCP<panzer::GlobalData>& global_data,
PHX::FieldManager<panzer::Traits>& fm) const
{
using std::string;
using std::vector;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ParameterList;
using PHX::Evaluator;
RCP< vector< RCP<Evaluator<panzer::Traits> > > > evaluators =
rcp(new vector< RCP<Evaluator<panzer::Traits> > > );
if (!models.isSublist(model_id)) {
models.print(std::cout);
std::stringstream msg;
msg << "Falied to find requested model, \"" << model_id
<< "\", for equation set:\n" << std::endl;
TEUCHOS_TEST_FOR_EXCEPTION(!models.isSublist(model_id), std::logic_error, msg.str());
}
std::vector<Teuchos::RCP<const panzer::PureBasis> > bases;
fl.uniqueBases(bases);
const ParameterList& my_models = models.sublist(model_id);
for (ParameterList::ConstIterator model_it = my_models.begin();
model_it != my_models.end(); ++model_it) {
bool found = false;
const std::string key = model_it->first;
ParameterList input;
const Teuchos::ParameterEntry& entry = model_it->second;
const ParameterList& plist = Teuchos::getValue<Teuchos::ParameterList>(entry);
#ifdef HAVE_STOKHOS
if (plist.isType<double>("Value") && plist.isType<double>("UQ")
&& plist.isParameter("Expansion")
&& (typeid(EvalT)==typeid(panzer::Traits::SGResidual) ||
typeid(EvalT)==typeid(panzer::Traits::SGJacobian)) ) {
{ // at IP
input.set("Name", key);
input.set("Value", plist.get<double>("Value"));
input.set("UQ", plist.get<double>("UQ"));
input.set("Expansion", plist.get<Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > >("Expansion"));
input.set("Data Layout", ir->dl_scalar);
RCP< Evaluator<panzer::Traits> > e =
rcp(new user_app::ConstantModel<EvalT,panzer::Traits>(input));
evaluators->push_back(e);
}
for (std::vector<Teuchos::RCP<const panzer::PureBasis> >::const_iterator basis_itr = bases.begin();
basis_itr != bases.end(); ++basis_itr) { // at BASIS
input.set("Name", key);
input.set("Value", plist.get<double>("Value"));
input.set("UQ", plist.get<double>("UQ"));
input.set("Expansion", plist.get<Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > >("Expansion"));
Teuchos::RCP<const panzer::BasisIRLayout> basis = basisIRLayout(*basis_itr,*ir);
input.set("Data Layout", basis->functional);
RCP< Evaluator<panzer::Traits> > e =
rcp(new user_app::ConstantModel<EvalT,panzer::Traits>(input));
evaluators->push_back(e);
}
found = true;
}
else
#endif
if (plist.isType<std::string>("Type")) {
if (plist.get<std::string>("Type") == "Parameter") {
{ // at IP
RCP< Evaluator<panzer::Traits> > e =
rcp(new panzer::Parameter<EvalT,panzer::Traits>(key,ir->dl_scalar,plist.get<double>("Value"),*global_data->pl));
evaluators->push_back(e);
}
for (std::vector<Teuchos::RCP<const panzer::PureBasis> >::const_iterator basis_itr = bases.begin();
basis_itr != bases.end(); ++basis_itr) { // at BASIS
Teuchos::RCP<const panzer::BasisIRLayout> basis = basisIRLayout(*basis_itr,*ir);
RCP< Evaluator<panzer::Traits> > e =
rcp(new panzer::Parameter<EvalT,panzer::Traits>(key,basis->functional,plist.get<double>("Value"),*global_data->pl));
evaluators->push_back(e);
}
found = true;
}
}
else if (plist.isType<double>("Value")) {
{ // at IP
//.........这里部分代码省略.........
示例15: initializePrec
void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
using Teuchos::rcp_dynamic_cast;
// we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
#ifdef HAVE_MUELU_TPETRA
typedef MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> MueTpOp;
typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpOp;
typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
#endif
// Check precondition
TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
TEUCHOS_ASSERT(prec);
// Create a copy, as we may remove some things from the list
ParameterList paramList = *paramList_;
// Retrieve wrapped concrete Xpetra matrix from FwdOp
const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
// Check whether it is Epetra/Tpetra
bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
RCP<XpMat> A = Teuchos::null;
if(bIsBlocked) {
Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
// MueLu needs a non-const object as input
RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
// wrap the forward operator as an Xpetra::Matrix that MueLu can work with
RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
RCP<const XpMap> rowmap00 = A00->getRowMap();
RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
// create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
// save blocked matrix
A = bMat;
} else {
RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
// MueLu needs a non-const object as input
RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
// wrap the forward operator as an Xpetra::Matrix that MueLu can work with
A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
}
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
// Retrieve concrete preconditioner object
const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
// extract preconditioner operator
RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
// Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy
RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;
// make a decision whether to (re)build the multigrid preconditioner or reuse the old one
// rebuild preconditioner if startingOver == true
// reuse preconditioner if startingOver == false
const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
//.........这里部分代码省略.........