本文整理汇总了C++中Teuchos::rcp方法的典型用法代码示例。如果您正苦于以下问题:C++ Teuchos::rcp方法的具体用法?C++ Teuchos::rcp怎么用?C++ Teuchos::rcp使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类Teuchos
的用法示例。
在下文中一共展示了Teuchos::rcp方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: if
TEUCHOS_UNIT_TEST(tUniqueGlobalIndexer_Utilities,ArrayToFieldVector)
{
typedef Intrepid::FieldContainer<int> IntFieldContainer;
// build global (or serial communicator)
#ifdef HAVE_MPI
RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
#else
RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
#endif
int myRank = eComm->MyPID();
int numProcs = eComm->NumProc();
TEUCHOS_ASSERT(numProcs==2);
RCP<panzer::UniqueGlobalIndexer<short,int> > globalIndexer
= rcp(new panzer::unit_test::UniqueGlobalIndexer(myRank,numProcs));
panzer::ArrayToFieldVector<short,int,KokkosClassic::DefaultNode::DefaultNodeType> atfv(globalIndexer);
std::map<std::string,IntFieldContainer> dataU, dataT;
fillFieldContainer(globalIndexer->getFieldNum("U"),"block_0",*globalIndexer,dataU["block_0"]);
fillFieldContainer(globalIndexer->getFieldNum("U"),"block_1",*globalIndexer,dataU["block_1"]);
fillFieldContainer(globalIndexer->getFieldNum("T"),"block_0",*globalIndexer,dataT["block_0"]);
Teuchos::RCP<Tpetra::MultiVector<int,int,int> > reducedUDataVector = atfv.getDataVector<int>("U",dataU);
Teuchos::RCP<Tpetra::MultiVector<int,int,int> > reducedTDataVector = atfv.getDataVector<int>("T",dataT);
std::vector<int> fields_u(reducedUDataVector->getLocalLength());
std::vector<int> fields_t(reducedTDataVector->getLocalLength());
reducedUDataVector->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(fields_u));
reducedTDataVector->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(fields_t));
if(myRank==0) {
TEST_EQUALITY(reducedUDataVector->getLocalLength(),6);
TEST_EQUALITY(reducedTDataVector->getLocalLength(),5);
TEST_EQUALITY(fields_u[0], 6);
TEST_EQUALITY(fields_u[1], 0);
TEST_EQUALITY(fields_u[2], 2);
TEST_EQUALITY(fields_u[3], 8);
TEST_EQUALITY(fields_u[4],10);
TEST_EQUALITY(fields_u[5],13);
TEST_EQUALITY(fields_t[0], 7);
TEST_EQUALITY(fields_t[1], 1);
TEST_EQUALITY(fields_t[2], 3);
TEST_EQUALITY(fields_t[3], 9);
TEST_EQUALITY(fields_t[4],11);
}
else if(myRank==1) {
TEST_EQUALITY(reducedUDataVector->getLocalLength(),4);
TEST_EQUALITY(reducedTDataVector->getLocalLength(),1);
TEST_EQUALITY(fields_u[0], 4);
TEST_EQUALITY(fields_u[1],12);
TEST_EQUALITY(fields_u[2],15);
TEST_EQUALITY(fields_u[3],14);
TEST_EQUALITY(fields_t[0], 5);
}
else
TEUCHOS_ASSERT(false);
}
示例2: if
//----------------------------------------------------------------------------
Teuchos::RCP<AAdapt::AbstractAdapter>
AAdapt::AdaptationFactory::createAdapter() {
using Teuchos::rcp;
Teuchos::RCP<AAdapt::AbstractAdapter> strategy;
std::string& method = adapt_params_->get("Method", "");
#if defined(ALBANY_STK)
if(method == "Copy Remesh") {
strategy = rcp(new AAdapt::CopyRemesh(adapt_params_,
param_lib_,
state_mgr_,
commT_));
}
#if defined(ALBANY_LCM) && defined(ALBANY_BGL)
else if(method == "Topmod") {
strategy = rcp(new AAdapt::TopologyMod(adapt_params_,
param_lib_,
state_mgr_,
commT_));
}
#endif
#if defined(ALBANY_LCM) && defined(LCM_SPECULATIVE)
else if(method == "Random") {
strategy = rcp(new AAdapt::RandomFracture(adapt_params_,
param_lib_,
state_mgr_,
commT_));
}
#endif
#endif
#if 0
#if defined(ALBANY_LCM) && defined(ALBANY_STK_PERCEPT)
else if(method == "Unif Size") {
strategy = rcp(new AAdapt::STKAdapt<AAdapt::STKUnifRefineField>(adapt_params_,
param_lib_,
state_mgr_,
epetra_comm_));
}
#endif
#endif
#if defined(ALBANY_STK)
else
#endif
TEUCHOS_TEST_FOR_EXCEPTION(true,
Teuchos::Exceptions::InvalidParameter,
std::endl <<
"Error! Unknown adaptivity method requested:"
<< method <<
" !" << std::endl
<< "Supplied parameter list is " <<
std::endl << *adapt_params_);
return strategy;
}
示例3: 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 verbose = false, debug = false, 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 maxiters = -1; // maximum number of iterations allowed per linear system
int maxsubspace = 50; // maximum number of blocks the solver can use for the subspace
int maxrestarts = 15; // number of restarts allowed
std::string filename("orsirr1.hb");
MT tol = 1.0e-5; // relative residual tolerance
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("debug","nondebug",&debug,"Print debugging information from solver.");
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 -1;
}
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->Seed();
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);
}
//
// ********Other information used by block solver***********
// *****************(can be user specified)******************
//
const int NumGlobalElements = B->GlobalLength();
if (maxiters == -1)
maxiters = NumGlobalElements/blocksize - 1; // maximum number of iterations to run
//
ParameterList belosList;
belosList.set( "Num Blocks", maxsubspace); // Maximum number of blocks in Krylov factorization
belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
int verbosity = Belos::Errors + Belos::Warnings;
if (verbose) {
verbosity += Belos::TimingDetails + Belos::StatusTestDetails;
if (frequency > 0)
belosList.set( "Output Frequency", frequency );
}
if (debug) {
verbosity += Belos::Debug;
}
belosList.set( "Verbosity", verbosity );
//
// Construct an unpreconditioned linear problem instance.
//.........这里部分代码省略.........
示例4: if
/*! \brief Create a mesh of approximately the desired size.
*
* We want 3 dimensions close to equal in length.
*/
const RCP<tMVector_t> getMeshCoordinates(
const RCP<const Teuchos::Comm<int> > & comm,
zgno_t numGlobalCoords)
{
int rank = comm->getRank();
int nprocs = comm->getSize();
double k = log(numGlobalCoords) / 3;
double xdimf = exp(k) + 0.5;
ssize_t xdim = static_cast<ssize_t>(floor(xdimf));
ssize_t ydim = xdim;
ssize_t zdim = numGlobalCoords / (xdim*ydim);
ssize_t num=xdim*ydim*zdim;
ssize_t diff = numGlobalCoords - num;
ssize_t newdiff = 0;
while (diff > 0){
if (zdim > xdim && zdim > ydim){
zdim++;
newdiff = diff - (xdim*ydim);
if (newdiff < 0)
if (diff < -newdiff)
zdim--;
}
else if (ydim > xdim && ydim > zdim){
ydim++;
newdiff = diff - (xdim*zdim);
if (newdiff < 0)
if (diff < -newdiff)
ydim--;
}
else{
xdim++;
newdiff = diff - (ydim*zdim);
if (newdiff < 0)
if (diff < -newdiff)
xdim--;
}
diff = newdiff;
}
num=xdim*ydim*zdim;
diff = numGlobalCoords - num;
if (diff < 0)
diff /= -numGlobalCoords;
else
diff /= numGlobalCoords;
if (rank == 0){
if (diff > .01)
cout << "Warning: Difference " << diff*100 << " percent" << endl;
cout << "Mesh size: " << xdim << "x" << ydim << "x" <<
zdim << ", " << num << " vertices." << endl;
}
// Divide coordinates.
ssize_t numLocalCoords = num / nprocs;
ssize_t leftOver = num % nprocs;
ssize_t gid0 = 0;
if (rank <= leftOver)
gid0 = zgno_t(rank) * (numLocalCoords+1);
else
gid0 = (leftOver * (numLocalCoords+1)) +
((zgno_t(rank) - leftOver) * numLocalCoords);
if (rank < leftOver)
numLocalCoords++;
ssize_t gid1 = gid0 + numLocalCoords;
zgno_t *ids = new zgno_t [numLocalCoords];
if (!ids)
throw bad_alloc();
ArrayRCP<zgno_t> idArray(ids, 0, numLocalCoords, true);
for (ssize_t i=gid0; i < gid1; i++)
*ids++ = zgno_t(i);
RCP<const tMap_t> idMap = rcp(
new tMap_t(num, idArray.view(0, numLocalCoords), 0, comm));
// Create a Tpetra::MultiVector of coordinates.
zscalar_t *x = new zscalar_t [numLocalCoords*3];
if (!x)
throw bad_alloc();
ArrayRCP<zscalar_t> coordArray(x, 0, numLocalCoords*3, true);
zscalar_t *y = x + numLocalCoords;
zscalar_t *z = y + numLocalCoords;
zgno_t xStart = 0;
zgno_t yStart = 0;
//.........这里部分代码省略.........
示例5: 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());
//.........这里部分代码省略.........
示例6: createResponseFunction
void
Albany::ResponseFactory::
createResponseFunction(
const std::string& name,
Teuchos::ParameterList& responseParams,
Teuchos::Array< Teuchos::RCP<AbstractResponseFunction> >& responses) const
{
using std::string;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ParameterList;
using Teuchos::Array;
RCP<const Epetra_Comm> comm = app->getComm();
if (name == "Solution Average") {
responses.push_back(rcp(new Albany::SolutionAverageResponseFunction(comm)));
}
else if (name == "Solution Two Norm") {
responses.push_back(rcp(new Albany::SolutionTwoNormResponseFunction(comm)));
}
else if (name == "Solution Values") {
responses.push_back(rcp(new Albany::SolutionValuesResponseFunction(app, responseParams)));
}
else if (name == "Solution Max Value") {
int eq = responseParams.get("Equation", 0);
int neq = responseParams.get("Num Equations", 1);
bool inor = responseParams.get("Interleaved Ordering", true);
responses.push_back(
rcp(new Albany::SolutionMaxValueResponseFunction(comm, neq, eq, inor)));
}
else if (name == "Solution Two Norm File") {
responses.push_back(
rcp(new Albany::SolutionFileResponseFunction<Albany::NormTwo>(comm)));
}
else if (name == "Solution Inf Norm File") {
responses.push_back(
rcp(new Albany::SolutionFileResponseFunction<Albany::NormInf>(comm)));
}
else if (name == "Aggregated") {
int num_aggregate_responses = responseParams.get<int>("Number");
Array< RCP<AbstractResponseFunction> > aggregated_responses;
Array< RCP<ScalarResponseFunction> > scalar_responses;
for (int i=0; i<num_aggregate_responses; i++) {
std::string id = Albany::strint("Response",i);
std::string name = responseParams.get<std::string>(id);
std::string sublist_name = Albany::strint("ResponseParams",i);
createResponseFunction(name, responseParams.sublist(sublist_name),
aggregated_responses);
}
scalar_responses.resize(aggregated_responses.size());
for (int i=0; i<aggregated_responses.size(); i++) {
TEUCHOS_TEST_FOR_EXCEPTION(
aggregated_responses[i]->isScalarResponse() != true, std::logic_error,
"Response function " << i << " is not a scalar response function." <<
std::endl <<
"The aggregated response can only aggregate scalar response " <<
"functions!");
scalar_responses[i] =
Teuchos::rcp_dynamic_cast<ScalarResponseFunction>(
aggregated_responses[i]);
}
responses.push_back(
rcp(new Albany::AggregateScalarResponseFunction(comm, scalar_responses)));
}
else if (name == "Field Integral" ||
name == "Field Value" ||
name == "Field Average" ||
name == "Surface Velocity Mismatch" ||
name == "Center Of Mass" ||
name == "Save Field" ||
name == "Region Boundary" ||
name == "Element Size Field" ||
name == "IP to Nodal Field" ||
name == "Save Nodal Fields" ||
name == "PHAL Field Integral") {
responseParams.set("Name", name);
for (int i=0; i<meshSpecs.size(); i++) {
responses.push_back(
rcp(new Albany::FieldManagerScalarResponseFunction(
app, prob, meshSpecs[i], stateMgr, responseParams)));
}
}
else if (name == "Solution") {
responses.push_back(
rcp(new Albany::SolutionResponseFunction(app, responseParams)));
}
else if (name == "KL") {
Array< RCP<AbstractResponseFunction> > base_responses;
//.........这里部分代码省略.........
示例7: Compute
int Ifpack_SORa::Compute(){
if(!IsInitialized_) Initialize();
Epetra_Map *RowMap=const_cast<Epetra_Map*>(&A_->RowMatrixRowMap());
Epetra_Vector Adiag(*RowMap);
Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
double *vals_s,*vals_h;
bool RowMatrixMode=(Acrs_==Teuchos::null);
// Label
sprintf(Label_, "IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_);
if(RowMatrixMode){
if(!A_->Comm().MyPID()) cout<<"SORa: RowMatrix mode enabled"<<endl;
// RowMatrix mode, build Acrs_ the hard way.
Epetra_RowMatrix *Arow=&*A_;
Epetra_Map *ColMap=const_cast<Epetra_Map*>(&A_->RowMatrixColMap());
int Nmax=A_->MaxNumEntries();
int length;
vector<int> indices(Nmax);
vector<double> values(Nmax);
Epetra_CrsMatrix *Acrs=new Epetra_CrsMatrix(Copy,*RowMap,Nmax);
for(int i=0;i<Arow->NumMyRows();i++){
Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices[0]);
for(int j=0;j<length;j++)
indices[j]=ColMap->GID(indices[j]);
Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices[0]);
}
Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
Acrs_=rcp(Acrs,true);
}
// Create Askew2
// Note: Here I let EpetraExt do the thinking for me. Since it gets all the maps correct for the E + F^T stencil.
// There are probably more efficient ways to do this but this method has the bonus of allowing code reuse.
IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,-1,Askew2));
Askew2->FillComplete();
// Create Aherm2
IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,1,Aherm2));
Aherm2->FillComplete();
int nnz2=Askew2->NumMyNonzeros();
int N=Askew2->NumMyRows();
// Grab pointers
IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
// Sanity checking: Make sure the sparsity patterns are the same
#define SANITY_CHECK
#ifdef SANITY_CHECK
for(int i=0;i<N;i++)
if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
for(int i=0;i<nnz2;i++)
if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
#endif
// Dirichlet Detection & Nuking of Aherm2 and Askew2
// Note: Relies on Aherm2/Askew2 having identical sparsity patterns (see sanity check above)
if(HaveOAZBoundaries_){
int numBCRows;
int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
delete [] dirRows;
delete dirCols;
}
// Grab diagonal of A
A_->ExtractDiagonalCopy(Adiag);
// Allocate the diagonal for W
Epetra_Vector *Wdiag = new Epetra_Vector(*RowMap);
// Build the W matrix (lower triangle only)
// Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
int maxentries=Askew2->MaxNumEntries();
int* gids=new int [maxentries];
double* newvals=new double[maxentries];
W=new Epetra_CrsMatrix(Copy,*RowMap,0);
for(int i=0;i<N;i++){
// Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
int rowgid=Acrs_->GRID(i);
double c_data=0.0;
double ipdamp=0.0;
int idx=0;
for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){
int colgid=Askew2->GCID(colind_s[j]);
c_data+=fabs(vals_s[j]);
if(rowgid>colgid){
// Rely on the fact that off-diagonal entries are always numbered last, dropping the entry entirely.
if(colind_s[j] < N) {
gids[idx]=colgid;
newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
//.........这里部分代码省略.........
示例8: main
int main(int argc, char *argv[]) {
int status=0; // 0 = pass, failures are incremented
// Initialize MPI and timer
int Proc=0;
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
double total_time = -MPI_Wtime();
(void) MPI_Comm_rank(MPI_COMM_WORLD, &Proc);
MPI_Comm appComm = MPI_COMM_WORLD;
#else
int appComm=0;
#endif
using Teuchos::RCP;
using Teuchos::rcp;
// Command-line argument for input file
//char* defaultfile="input_1.xml";
try {
RCP<EpetraExt::ModelEvaluator> Model = rcp(new MockModelEval_A(appComm));
// Set input arguments to evalModel call
EpetraExt::ModelEvaluator::InArgs inArgs = Model->createInArgs();
RCP<Epetra_Vector> x = rcp(new Epetra_Vector(*(Model->get_x_init())));
inArgs.set_x(x);
int num_p = inArgs.Np(); // Number of *vectors* of parameters
RCP<Epetra_Vector> p1;
if (num_p > 0) {
p1 = rcp(new Epetra_Vector(*(Model->get_p_init(0))));
inArgs.set_p(0,p1);
}
int numParams = p1->MyLength(); // Number of parameters in p1 vector
// Set output arguments to evalModel call
EpetraExt::ModelEvaluator::OutArgs outArgs = Model->createOutArgs();
RCP<Epetra_Vector> f = rcp(new Epetra_Vector(x->Map()));
outArgs.set_f(f);
int num_g = outArgs.Ng(); // Number of *vectors* of responses
RCP<Epetra_Vector> g1;
if (num_g > 0) {
g1 = rcp(new Epetra_Vector(*(Model->get_g_map(0))));
outArgs.set_g(0,g1);
}
RCP<Epetra_Operator> W_op = Model->create_W();
outArgs.set_W(W_op);
RCP<Epetra_MultiVector> dfdp = rcp(new Epetra_MultiVector(
*(Model->get_x_map()), numParams));
outArgs.set_DfDp(0, dfdp);
RCP<Epetra_MultiVector> dgdp = rcp(new Epetra_MultiVector(g1->Map(), numParams));
outArgs.set_DgDp(0, 0, dgdp);
RCP<Epetra_MultiVector> dgdx = rcp(new Epetra_MultiVector(x->Map(), g1->MyLength()));
outArgs.set_DgDx(0, dgdx);
// Now, evaluate the model!
Model->evalModel(inArgs, outArgs);
// Print out everything
if (Proc == 0)
cout << "Finished Model Evaluation: Printing everything {Exact in brackets}"
<< "\n-----------------------------------------------------------------"
<< std::setprecision(9) << endl;
x->Print(cout << "\nSolution vector! {3,3,3,3}\n");
if (num_p>0) p1->Print(cout << "\nParameters! {1,1}\n");
f->Print(cout << "\nResidual! {8,5,0,-7}\n");
if (num_g>0) g1->Print(cout << "\nResponses! {2}\n");
RCP<Epetra_CrsMatrix> W = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_op, true);
W->Print(cout << "\nJacobian! {6 on diags}\n");
dfdp->Print(cout << "\nDfDp sensitivity MultiVector! {-1,0,0,0}{0,-4,-6,-8}\n");
dgdp->Print(cout << "\nDgDp response sensitivity MultiVector!{2,2}\n");
dgdx->Print(cout << "\nDgDx^T response gradient MultiVector! {-2,-2,-2,-2}\n");
if (Proc == 0)
cout <<
"\n-----------------------------------------------------------------\n";
}
catch (std::exception& e) {
cout << e.what() << endl;
status = 10;
}
catch (string& s) {
cout << s << endl;
status = 20;
}
catch (char *s) {
cout << s << endl;
status = 30;
}
//.........这里部分代码省略.........
示例9: TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(initial_condition_control, control)
{
using Teuchos::RCP;
using Teuchos::rcp;
Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
// setup mesh
/////////////////////////////////////////////
RCP<panzer_stk_classic::STK_Interface> mesh;
{
RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
pl->set<int>("X Elements",2);
pl->set<int>("Y Elements",2);
pl->set<int>("X Blocks",2);
pl->set<int>("Y Blocks",1);
panzer_stk_classic::SquareQuadMeshFactory mesh_factory;
mesh_factory.setParameterList(pl);
mesh = mesh_factory.buildMesh(MPI_COMM_WORLD);
mesh->writeToExodus("test.exo");
}
RCP<const shards::CellTopology> ct = mesh->getCellTopology("eblock-0_0");
panzer::CellData cellData(4,ct);
RCP<panzer::IntegrationRule> int_rule = rcp(new panzer::IntegrationRule(2,cellData));
ICFieldDescriptor densDesc;
densDesc.fieldName = "DENSITY";
densDesc.basisName = "Const";
densDesc.basisOrder = 0;
ICFieldDescriptor condDesc;
condDesc.fieldName = "CONDUCTIVITY";
condDesc.basisName = "HGrad";
condDesc.basisOrder = 1;
RCP<panzer::PureBasis> const_basis = rcp(new panzer::PureBasis(densDesc.basisName,densDesc.basisOrder,cellData));
RCP<panzer::PureBasis> hgrad_basis = rcp(new panzer::PureBasis(condDesc.basisName,condDesc.basisOrder,cellData));
RCP<const panzer::FieldPattern> constFP = rcp(new panzer::Intrepid2FieldPattern(const_basis->getIntrepid2Basis()));
RCP<const panzer::FieldPattern> hgradFP = rcp(new panzer::Intrepid2FieldPattern(hgrad_basis->getIntrepid2Basis()));
// setup DOF manager
/////////////////////////////////////////////
RCP<panzer::ConnManager<int,int> > conn_manager
= Teuchos::rcp(new panzer_stk_classic::STKConnManager<int>(mesh));
RCP<panzer::DOFManager<int,int> > dofManager
= rcp(new panzer::DOFManager<int,int>(conn_manager,MPI_COMM_WORLD));
dofManager->addField(densDesc.fieldName, constFP);
dofManager->addField(condDesc.fieldName, hgradFP);
dofManager->buildGlobalUnknowns();
Teuchos::RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > elof
= Teuchos::rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(comm.getConst(),dofManager));
Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof = elof;
// setup worksets
/////////////////////////////////////////////
std::map<std::string,panzer::WorksetNeeds> needs;
needs["eblock-0_0"].cellData = cellData;
needs["eblock-0_0"].int_rules.push_back(int_rule);
needs["eblock-0_0"].bases = { const_basis, hgrad_basis};
needs["eblock-0_0"].rep_field_name = { densDesc.fieldName, condDesc.fieldName};
needs["eblock-1_0"].cellData = cellData;
needs["eblock-1_0"].int_rules.push_back(int_rule);
needs["eblock-1_0"].bases = { const_basis, hgrad_basis};
needs["eblock-1_0"].rep_field_name = { densDesc.fieldName, condDesc.fieldName};
Teuchos::RCP<panzer_stk_classic::WorksetFactory> wkstFactory
= Teuchos::rcp(new panzer_stk_classic::WorksetFactory(mesh)); // build STK workset factory
Teuchos::RCP<panzer::WorksetContainer> wkstContainer // attach it to a workset container (uses lazy evaluation)
= Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,needs));
// setup field manager builder
/////////////////////////////////////////////
// Add in the application specific closure model factory
panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory;
user_app::STKModelFactory_TemplateBuilder cm_builder;
cm_factory.buildObjects(cm_builder);
Teuchos::ParameterList user_data("User Data");
Teuchos::ParameterList ic_closure_models("Initial Conditions");
ic_closure_models.sublist("eblock-0_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
ic_closure_models.sublist("eblock-0_0").sublist(condDesc.fieldName).set<double>("Value",9.0);
ic_closure_models.sublist("eblock-1_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
ic_closure_models.sublist("eblock-1_0").sublist(condDesc.fieldName).set<double>("Value",9.0);
std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
block_ids_to_cell_topo["eblock-0_0"] = mesh->getCellTopology("eblock-0_0");
block_ids_to_cell_topo["eblock-1_0"] = mesh->getCellTopology("eblock-1_0");
//.........这里部分代码省略.........
示例10: main
// =============================================================================
int main (int argc, char *argv[])
{
Teuchos::GlobalMPISession(&argc, &argv, NULL);
// Create a communicator for Epetra objects
#ifdef HAVE_MPI
Epetra_MpiComm eComm(MPI_COMM_WORLD);
#else
Epetra_SerialComm eComm();
#endif
const RCP<Teuchos::FancyOStream> out =
Teuchos::VerboseObjectBase::getDefaultOStream();
bool success = true;
try {
// ===========================================================================
// handle command line arguments
Teuchos::CommandLineProcessor My_CLP;
My_CLP.setDocString("Linear solver testbed for the 1D Poisson matrix.\n");
std::string action("matvec");
My_CLP.setOption("action",
&action,
"Which action to perform with the operator (matvec, solve_cg, solve_minres, solve_gmres)"
);
std::string solver("cg");
// My_CLP.setOption("solver", &solver, "Krylov subspace method (cg, minres, gmres)");
// Operator op = JAC;
// Operator allOpts[] = {JAC, KEO, KEOREG, POISSON1D};
// std::string allOptNames[] = {"jac", "keo", "keoreg", "poisson1d"};
// My_CLP.setOption("operator", &op, 4, allOpts, allOptNames);
bool verbose = true;
My_CLP.setOption("verbose", "quiet",
&verbose,
"Print messages and results.");
int frequency = 10;
My_CLP.setOption("frequency",
&frequency,
"Solvers frequency for printing residuals (#iters).");
// Make sure this value is large enough to keep the cores busy for a while.
int n = 1000;
My_CLP.setOption("size",
&n,
"Size of the equation system (default: 1000).");
// print warning for unrecognized arguments
My_CLP.recogniseAllOptions(true);
My_CLP.throwExceptions(false);
// finally, parse the command line
TEUCHOS_ASSERT_EQUALITY(My_CLP.parse (argc, argv),
Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL
);
// =========================================================================
// Construct Epetra matrix.
RCP<Teuchos::Time> matrixConstructTime =
Teuchos::TimeMonitor::getNewTimer("Epetra matrix construction");
RCP<Epetra_CrsMatrix> A;
{
Teuchos::TimeMonitor tm(*matrixConstructTime);
A = contructEpetraMatrix(n, eComm);
}
// A->Print(std::cout);
// create initial guess and right-hand side
RCP<Epetra_Vector> x =
rcp(new Epetra_Vector(A->OperatorDomainMap()));
RCP<Epetra_MultiVector> b =
rcp(new Epetra_Vector(A->OperatorRangeMap()));
// b->Random();
TEUCHOS_ASSERT_EQUALITY(0, b->PutScalar(1.0));
if (action.compare("matvec") == 0) {
TEUCHOS_ASSERT_EQUALITY(0, x->PutScalar(1.0));
RCP<Teuchos::Time> mvTime = Teuchos::TimeMonitor::getNewTimer("Epetra operator apply");
{
Teuchos::TimeMonitor tm(*mvTime);
// Don't TEUCHOS_ASSERT_EQUALITY() here for speed.
A->Apply(*x, *b);
}
// print timing data
Teuchos::TimeMonitor::summarize();
} else {
// -----------------------------------------------------------------------
// Belos part
Teuchos::ParameterList belosList;
// Relative convergence tolerance requested
belosList.set("Convergence Tolerance", 1.0e-12);
if (verbose) {
belosList.set("Verbosity",
Belos::Errors +
//.........这里部分代码省略.........
示例11: main
int main(int argc, char *argv[])
{
#ifndef HAVE_TPETRA_COMPLEX_DOUBLE
# error "Anasazi: This test requires Scalar = std::complex<double> to be enabled in Tpetra."
#else
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::tuple;
using std::cout;
using std::endl;
typedef std::complex<double> ST;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Tpetra::MultiVector<ST> MV;
typedef MV::global_ordinal_type GO;
typedef Tpetra::Operator<ST> OP;
typedef Anasazi::MultiVecTraits<ST,MV> MVT;
typedef Anasazi::OperatorTraits<ST,MV,OP> OPT;
Teuchos::GlobalMPISession mpisess (&argc, &argv, &std::cout);
bool success = false;
const ST ONE = SCT::one ();
int info = 0;
RCP<const Teuchos::Comm<int> > comm =
Tpetra::DefaultPlatform::getDefaultPlatform ().getComm ();
const int MyPID = comm->getRank ();
bool verbose = false;
bool debug = false;
bool insitu = false;
bool herm = false;
std::string which("LM");
std::string filename;
int nev = 4;
int blockSize = 4;
MT tol = 1.0e-6;
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
cmdp.setOption("insitu","exsitu",&insitu,"Perform in situ restarting.");
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
cmdp.setOption("herm","nonherm",&herm,"Solve Hermitian or non-Hermitian problem.");
cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix (assumes non-Hermitian unless specified otherwise).");
cmdp.setOption("nev",&nev,"Number of eigenvalues to compute.");
cmdp.setOption("blockSize",&blockSize,"Block size for the algorithm.");
cmdp.setOption("tol",&tol,"Tolerance for convergence.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
if (debug) verbose = true;
if (filename == "") {
// get default based on herm
if (herm) {
filename = "mhd1280b.cua";
}
else {
filename = "mhd1280a.cua";
}
}
if (MyPID == 0) {
cout << Anasazi::Anasazi_Version() << endl << endl;
}
// Get the data from the HB file
int dim,dim2,nnz;
int rnnzmax;
double *dvals;
int *colptr,*rowind;
nnz = -1;
if (MyPID == 0) {
info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,&colptr,&rowind,&dvals);
// find maximum NNZ over all rows
vector<int> rnnz(dim,0);
for (int *ri=rowind; ri<rowind+nnz; ++ri) {
++rnnz[*ri-1];
}
rnnzmax = *std::max_element(rnnz.begin(),rnnz.end());
}
else {
// address uninitialized data warnings
dvals = NULL;
colptr = NULL;
rowind = NULL;
}
Teuchos::broadcast(*comm,0,&info);
Teuchos::broadcast(*comm,0,&nnz);
Teuchos::broadcast(*comm,0,&dim);
Teuchos::broadcast(*comm,0,&rnnzmax);
if (info == 0 || nnz < 0) {
if (MyPID == 0) {
cout << "Error reading '" << filename << "'" << endl
<< "End Result: TEST FAILED" << endl;
//.........这里部分代码省略.........
示例12: Initialize
void IncSVDPOD::Initialize( const Teuchos::RCP< Teuchos::ParameterList >& params,
const Teuchos::RCP< const Epetra_MultiVector >& ss,
const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > >& fileio ) {
using Teuchos::rcp;
// Get the "Reduced Basis Method" sublist.
Teuchos::ParameterList rbmethod_params = params->sublist( "Reduced Basis Method" );
// Get the maximum basis size
maxBasisSize_ = rbmethod_params.get<int>("Max Basis Size");
TEUCHOS_TEST_FOR_EXCEPTION(maxBasisSize_ < 2,std::invalid_argument,"\"Max Basis Size\" must be at least 2.");
// Get a filter
filter_ = rbmethod_params.get<Teuchos::RCP<Filter<double> > >("Filter",Teuchos::null);
if (filter_ == Teuchos::null) {
int k = rbmethod_params.get("Rank",(int)5);
filter_ = rcp( new RangeFilter<double>(LARGEST,k,k) );
}
// Get convergence tolerance
tol_ = rbmethod_params.get<double>("Convergence Tolerance",tol_);
// Get debugging flag
debug_ = rbmethod_params.get<bool>("IncSVD Debug",debug_);
// Get verbosity level
verbLevel_ = rbmethod_params.get<int>("IncSVD Verbosity Level",verbLevel_);
// Get an Anasazi orthomanager
if (rbmethod_params.isType<
Teuchos::RCP< Anasazi::OrthoManager<double,Epetra_MultiVector> >
>("Ortho Manager")
)
{
ortho_ = rbmethod_params.get<
Teuchos::RCP<Anasazi::OrthoManager<double,Epetra_MultiVector> >
>("Ortho Manager");
TEUCHOS_TEST_FOR_EXCEPTION(ortho_ == Teuchos::null,std::invalid_argument,"User specified null ortho manager.");
}
else {
std::string omstr = rbmethod_params.get("Ortho Manager","DGKS");
if (omstr == "SVQB") {
ortho_ = rcp( new Anasazi::SVQBOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
}
else { // if omstr == "DGKS"
ortho_ = rcp( new Anasazi::BasicOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
}
}
// Lmin,Lmax,Kstart
lmin_ = rbmethod_params.get("Min Update Size",1);
TEUCHOS_TEST_FOR_EXCEPTION(lmin_ < 1 || lmin_ >= maxBasisSize_,std::invalid_argument,
"Method requires 1 <= min update size < max basis size.");
lmax_ = rbmethod_params.get("Max Update Size",maxBasisSize_);
TEUCHOS_TEST_FOR_EXCEPTION(lmin_ > lmax_,std::invalid_argument,"Max update size must be >= min update size.");
startRank_ = rbmethod_params.get("Start Rank",lmin_);
TEUCHOS_TEST_FOR_EXCEPTION(startRank_ < 1 || startRank_ > maxBasisSize_,std::invalid_argument,
"Starting rank must be in [1,maxBasisSize_)");
// MaxNumPasses
maxNumPasses_ = rbmethod_params.get("Maximum Number Passes",maxNumPasses_);
TEUCHOS_TEST_FOR_EXCEPTION(maxNumPasses_ != -1 && maxNumPasses_ <= 0, std::invalid_argument,
"Maximum number passes must be -1 or > 0.");
// Save the pointer to the snapshot matrix
TEUCHOS_TEST_FOR_EXCEPTION(ss == Teuchos::null,std::invalid_argument,"Input snapshot matrix cannot be null.");
A_ = ss;
// MaxNumCols
maxNumCols_ = A_->NumVectors();
maxNumCols_ = rbmethod_params.get("Maximum Number Columns",maxNumCols_);
TEUCHOS_TEST_FOR_EXCEPTION(maxNumCols_ < A_->NumVectors(), std::invalid_argument,
"Maximum number of columns must be at least as many as in the initializing data set.");
// V locally replicated or globally distributed
// this must be true for now
// Vlocal_ = rbmethod_params.get("V Local",Vlocal_);
// Allocate space for the factorization
sigma_.reserve( maxBasisSize_ );
U_ = Teuchos::null;
V_ = Teuchos::null;
U_ = rcp( new Epetra_MultiVector(ss->Map(),maxBasisSize_,false) );
if (Vlocal_) {
Epetra_LocalMap lclmap(maxNumCols_,0,ss->Comm());
V_ = rcp( new Epetra_MultiVector(lclmap,maxBasisSize_,false) );
}
else {
Epetra_Map gblmap(maxNumCols_,0,ss->Comm());
V_ = rcp( new Epetra_MultiVector(gblmap,maxBasisSize_,false) );
}
B_ = rcp( new Epetra_SerialDenseMatrix(maxBasisSize_,maxBasisSize_) );
resNorms_.reserve(maxBasisSize_);
// clear counters
numProc_ = 0;
curNumPasses_ = 0;
// we are now initialized, albeit with null rank
isInitialized_ = true;
//.........这里部分代码省略.........
示例13: main
int main(int argc, char* argv[])
{
info("Desired number of eigenvalues: %d.", NUMBER_OF_EIGENVALUES);
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("domain.mesh", &mesh);
// Perform initial mesh refinements (optional).
for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_dirichlet(BDY_ALL);
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_zero(BDY_ALL);
// Create an H1 space with default shapeset.
H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
int ndof = Space::get_num_dofs(&space);
info("ndof: %d.", ndof);
// Initialize the weak formulation for the left hand side, i.e., H.
WeakForm wf_left, wf_right;
wf_left.add_matrix_form(callback(bilinear_form_left));
wf_right.add_matrix_form(callback(bilinear_form_right));
// Initialize matrices.
RCP<SparseMatrix> matrix_left = rcp(new CSCMatrix());
RCP<SparseMatrix> matrix_right = rcp(new CSCMatrix());
// Assemble the matrices.
bool is_linear = true;
DiscreteProblem dp_left(&wf_left, &space, is_linear);
dp_left.assemble(matrix_left.get());
DiscreteProblem dp_right(&wf_right, &space, is_linear);
dp_right.assemble(matrix_right.get());
EigenSolver es(matrix_left, matrix_right);
info("Calling Pysparse...");
es.solve(NUMBER_OF_EIGENVALUES, TARGET_VALUE, TOL, MAX_ITER);
info("Pysparse finished.");
es.print_eigenvalues();
// Initializing solution vector, solution and ScalarView.
double* coeff_vec;
Solution sln;
ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
// Reading solution vectors and visualizing.
double* eigenval = new double[NUMBER_OF_EIGENVALUES];
int neig = es.get_n_eigs();
if (neig != NUMBER_OF_EIGENVALUES) error("Mismatched number of eigenvectors in the eigensolver output file.");
for (int ieig = 0; ieig < neig; ieig++) {
eigenval[ieig] = es.get_eigenvalue(ieig);
int n;
es.get_eigenvector(ieig, &coeff_vec, &n);
// Convert coefficient vector into a Solution.
Solution::vector_to_solution(coeff_vec, &space, &sln);
// Visualize the solution.
char title[100];
sprintf(title, "Solution %d, val = %g", ieig, eigenval[ieig]);
view.set_title(title);
view.show(&sln);
// Wait for keypress.
View::wait(HERMES_WAIT_KEYPRESS);
}
return 0;
};
示例14: MyOp
//
// Constructor
//
MyOp(const GlobalOrdinal n, const RCP< const Teuchos::Comm<int> > comm)
{
//
// Construct a map for our block row distribution
//
opMap_ = rcp( new Map(n, 0, comm) );
//
// Get the rank of this process and the number of processes
// We're going to have to do something special with the first and last processes
//
myRank_ = comm->getRank();
numProcs_ = comm->getSize();
//
// Get the local number of rows
//
LocalOrdinal nlocal = opMap_->getNodeNumElements();
//
// Define the distribution that you need for the matvec.
// When you define this for your own operator, it is helpful to draw pictures
// on a sheet of paper to keep track of who needs to receive which elements.
// Here, each process needs to receive one element from each of its neighbors.
//
// All processes but the first will receive one entry from the
// previous process.
if (myRank_ > 0) {
nlocal++;
}
// All processes but the last will receive one entry from the next
// process.
if (myRank_ < numProcs_-1) {
nlocal++;
}
// Construct a list of columns where this process has nonzero elements.
// For our tridiagonal matrix, this is firstRowItOwns-1:lastRowItOwns+1.
std::vector<GlobalOrdinal> indices;
indices.reserve (nlocal);
// The first process is a special case...
if (myRank_ > 0) {
indices.push_back (opMap_->getMinGlobalIndex () - 1);
}
for (GlobalOrdinal i = opMap_->getMinGlobalIndex ();
i <= opMap_->getMaxGlobalIndex (); ++i) {
indices.push_back (i);
}
// So is the last process...
if (myRank_ < numProcs_-1) {
indices.push_back (opMap_->getMaxGlobalIndex () + 1);
}
// Wrap our vector in an array view, which is like a raw pointer.
Teuchos::ArrayView<const GlobalOrdinal> elementList (indices);
// Make a Map for handling the redistribution. There will be some
// redundancies (i.e., some of the entries will be owned by
// multiple MPI processes).
GlobalOrdinal numGlobalElements = n + 2*(numProcs_-1);
redistMap_ = rcp (new Map (numGlobalElements, elementList, 0, comm));
// Make an Import object that describes how data will be
// redistributed. It takes a Map describing who owns what
// originally, and a Map that describes who you WANT to own what.
importer_= rcp (new Import (opMap_, redistMap_));
}
示例15: run
bool run( const Teuchos::RCP<const Teuchos::Comm<int> > & comm ,
const CMD & cmd)
{
typedef typename Kokkos::Compat::KokkosDeviceWrapperNode<Device> NodeType;
bool success = true;
try {
const int comm_rank = comm->getRank();
// Create Tpetra Node -- do this first as it initializes host/device
Teuchos::RCP<NodeType> node = createKokkosNode<NodeType>( cmd , *comm );
// Set up stochastic discretization
using Teuchos::Array;
using Teuchos::RCP;
using Teuchos::rcp;
typedef Stokhos::OneDOrthogPolyBasis<int,double> one_d_basis;
typedef Stokhos::LegendreBasis<int,double> legendre_basis;
typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > order_type;
typedef Stokhos::TotalOrderBasis<int,double,order_type> product_basis;
typedef Stokhos::Sparse3Tensor<int,double> Cijk;
const int dim = cmd.CMD_USE_UQ_DIM;
const int order = cmd.CMD_USE_UQ_ORDER ;
Array< RCP<const one_d_basis> > bases(dim);
for (int i=0; i<dim; i++)
bases[i] = rcp(new legendre_basis(order, true));
RCP<const product_basis> basis = rcp(new product_basis(bases));
RCP<Cijk> cijk = basis->computeTripleProductTensor();
typedef Stokhos::DynamicStorage<int,double,Device> Storage;
typedef Sacado::UQ::PCE<Storage> Scalar;
typename Scalar::cijk_type kokkos_cijk =
Stokhos::create_product_tensor<Device>(*basis, *cijk);
Kokkos::setGlobalCijkTensor(kokkos_cijk);
// typedef Stokhos::TensorProductQuadrature<int,double> quadrature;
// RCP<const quadrature> quad = rcp(new quadrature(basis));
// const int num_quad_points = quad->size();
// const Array<double>& quad_weights = quad->getQuadWeights();
// const Array< Array<double> >& quad_points = quad->getQuadPoints();
// const Array< Array<double> >& quad_values = quad->getBasisAtQuadPoints();
// Print output headers
const std::vector< size_t > widths =
print_headers( std::cout , cmd , comm_rank );
using Kokkos::Example::FENL::TrivialManufacturedSolution;
using Kokkos::Example::FENL::ElementComputationKLCoefficient;
using Kokkos::Example::BoxElemPart;
using Kokkos::Example::FENL::fenl;
using Kokkos::Example::FENL::Perf;
const double bc_lower_value = 1 ;
const double bc_upper_value = 2 ;
const TrivialManufacturedSolution manufactured_solution;
int nelem[3] = { cmd.CMD_USE_FIXTURE_X ,
cmd.CMD_USE_FIXTURE_Y ,
cmd.CMD_USE_FIXTURE_Z };
// Create KL diffusion coefficient
const double kl_mean = cmd.CMD_USE_MEAN;
const double kl_variance = cmd.CMD_USE_VAR;
const double kl_correlation = cmd.CMD_USE_COR;
typedef ElementComputationKLCoefficient< Scalar, double, Device > KL;
KL diffusion_coefficient( kl_mean, kl_variance, kl_correlation, dim );
typedef typename KL::RandomVariableView RV;
typedef typename RV::HostMirror HRV;
RV rv = diffusion_coefficient.getRandomVariables();
HRV hrv = Kokkos::create_mirror_view(rv);
// Set random variables
// ith random variable \xi_i = \psi_I(\xi) / \psi_I(1.0)
// where I is determined by the basis ordering (since the component basis
// functions have unit two-norm, \psi_I(1.0) might not be 1.0). We compute
// this by finding the index of the multivariate term that is first order in
// the ith slot, all other orders 0
Teuchos::Array<double> point(dim, 1.0);
Teuchos::Array<double> basis_vals(basis->size());
basis->evaluateBases(point, basis_vals);
for (int i=0; i<dim; ++i) {
Stokhos::MultiIndex<int> term(dim, 0);
term[i] = 1;
int index = basis->index(term);
hrv(i).fastAccessCoeff(index) = 1.0 / basis_vals[index];
}
Kokkos::deep_copy( rv, hrv );
// Compute stochastic response using stochastic Galerkin method
Scalar response = 0;
Perf perf;
if ( cmd.CMD_USE_FIXTURE_QUADRATIC )
perf = fenl< Scalar , Device , BoxElemPart::ElemQuadratic >
( comm , node , cmd.CMD_PRINT , cmd.CMD_USE_TRIALS ,
cmd.CMD_USE_ATOMIC , cmd.CMD_USE_BELOS , cmd.CMD_USE_MUELU ,
cmd.CMD_USE_MEANBASED ,
nelem , diffusion_coefficient , manufactured_solution ,
bc_lower_value , bc_upper_value ,
false , response);
else
//.........这里部分代码省略.........