46#include "Teuchos_as.hpp"
49 Teuchos::RCP<EpetraExt::ModelEvaluator> underlyingME_,
50 const Teuchos::RCP<EpetraExt::MultiComm> &globalComm_,
51 const std::vector<Epetra_Vector*> initGuessVec_,
52 Teuchos::RCP<std::vector< Teuchos::RCP<Epetra_Vector> > > q_vec_,
53 Teuchos::RCP<std::vector< Teuchos::RCP<Epetra_Vector> > > matching_vec_
55 underlyingME(underlyingME_),
56 globalComm(globalComm_),
59 timeStepsOnTimeDomain(globalComm_->NumTimeStepsOnDomain()),
60 numTimeDomains(globalComm_->NumSubDomains()),
61 timeDomain(globalComm_->SubDomainRank()),
62#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
65#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
68#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
71#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
74 matching_vec(matching_vec_)
78 std::cout <<
"----------MultiPoint Partition Info------------"
80 <<
"\n\tSpatial Decomposition = " <<
globalComm->SubDomainComm().NumProc()
83 <<
"\n\tTotal Number of Steps = " <<
globalComm->NumTimeSteps();
84 std::cout <<
"\n-----------------------------------------------" << std::endl;
92#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
93 if(
split_W->RowMatrixRowMap().GlobalIndicesInt()) {
98 (*rowStencil_int)[i].push_back(0);
99 (*rowIndex_int).push_back(i +
globalComm->FirstTimeStepOnDomain());
106#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
107 if(
split_W->RowMatrixRowMap().GlobalIndicesInt()) {
112 (*rowStencil_LL)[i].push_back(0);
113 (*rowIndex_LL).push_back(i +
globalComm->FirstTimeStepOnDomain());
120 throw "EpetraExt::MultiPointModelEvaluator::MultiPointModelEvaluator: Global indices unknown";
134 TEUCHOS_TEST_FOR_EXCEPT(underlyingOutArgs.
Np()!=2);
138 num_p0 = underlyingME_->get_p_map(0)->NumMyElements();
186#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
192#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
205 TEUCHOS_TEST_FOR_EXCEPT(!(*
matching_vec)[0]->Map().SameAs(*(underlyingME_->get_g_map(0))));
206 TEUCHOS_TEST_FOR_EXCEPT(
num_g0 != 1);
215 if (underlyingNg)
delete block_DgDx;
216#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
217 delete rowStencil_int;
220#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
221 delete rowStencil_LL;
237 return Teuchos::rcp(&(block_W->OperatorDomainMap()),
false);
247 return underlyingME->get_p_map(l);
252 return underlyingME->get_g_map(j);
257 return solution_init;
262 return underlyingME->get_p_init(l);
290 DERIV_LINEARITY_NONCONST
295 outArgs.
setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
298 DERIV_LINEARITY_CONST
299 ,DERIV_RANK_DEFICIENT
305 outArgs.
setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
308 DERIV_LINEARITY_NONCONST
309 ,DERIV_RANK_DEFICIENT
313 outArgs.
setSupports(OUT_ARG_DgDp,0,0, orientation_DgDp);
316 DERIV_LINEARITY_NONCONST
317 ,DERIV_RANK_DEFICIENT
339 Teuchos::RCP<const Epetra_Vector> p_in = inArgs.
get_p(0);
340 if (p_in.get()) underlyingInArgs.
set_p(0, p_in);
342 Teuchos::RCP<const Epetra_Vector> x_in = inArgs.
get_x();
343 block_x->Epetra_Vector::operator=(*x_in);
346 Teuchos::RCP<Epetra_Vector> f_out = outArgs.
get_f();
348 Teuchos::RCP<Epetra_Operator> W_out = outArgs.
get_W();
349 Teuchos::RCP<EpetraExt::BlockCrsMatrix> W_block =
350 Teuchos::rcp_dynamic_cast<EpetraExt::BlockCrsMatrix>(W_out);
352 Teuchos::RCP<Epetra_Vector> g_out;
353 if (underlyingNg) g_out = outArgs.
get_g(0);
354 if (g_out.get()) g_out->PutScalar(0.0);
369 bool need_g = g_out.get();
375 for (
int i=0; i < timeStepsOnTimeDomain; i++) {
378 underlyingInArgs.
set_p(1, (*q_vec)[i]);
382#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
383 block_x->ExtractBlockValues(*split_x, (*rowIndex_LL)[i]);
387#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
388 block_x->ExtractBlockValues(*split_x, (*rowIndex_int)[i]);
391 underlyingInArgs.
set_x(split_x);
394 if (f_out.get()) underlyingOutArgs.
set_f(split_f);
396 if (need_g) underlyingOutArgs.
set_g(0, split_g);
398 if (W_out.get()) underlyingOutArgs.
set_W(split_W);
404 if (!DgDp_out.
isEmpty()) underlyingOutArgs.
set_DgDp(0, 0, *deriv_DgDp);
407 underlyingME->evalModel(underlyingInArgs, underlyingOutArgs);
411 if (matchingProblem) {
413 double diff = (*split_g)[0] - (*(*matching_vec)[i])[0];
414 double nrmlz = fabs((*(*matching_vec)[i])[0]) + 1.0e-6;
415 (*split_g)[0] = 0.5 * diff * diff/(nrmlz*nrmlz);
416 if (!DgDx_out.
isEmpty()) split_DgDx->Scale(diff/(nrmlz*nrmlz));
417 if (!DgDp_out.
isEmpty()) split_DgDp->Scale(diff/(nrmlz*nrmlz));
423#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
424 if (f_out.get()) block_f->LoadBlockValues(*split_f, (*rowIndex_LL)[i]);
425 if (W_out.get()) W_block->LoadBlock(*split_W, i, 0);
427 if (!DfDp_out.
isEmpty()) block_DfDp->LoadBlockValues(*split_DfDp, (*rowIndex_LL)[i]);
428 if (!DgDx_out.
isEmpty()) block_DgDx->LoadBlockValues(*split_DgDx, (*rowIndex_LL)[i]);
432#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
433 if (f_out.get()) block_f->LoadBlockValues(*split_f, (*rowIndex_int)[i]);
434 if (W_out.get()) W_block->LoadBlock(*split_W, i, 0);
436 if (!DfDp_out.
isEmpty()) block_DfDp->LoadBlockValues(*split_DfDp, (*rowIndex_int)[i]);
437 if (!DgDx_out.
isEmpty()) block_DgDx->LoadBlockValues(*split_DgDx, (*rowIndex_int)[i]);
442 if (g_out.get()) g_out->Update(1.0, *split_g, 1.0);
450 if (f_out.get()) f_out->operator=(*block_f);
457 if (numTimeDomains > 1) {
458 double factorToZeroOutCopies = 0.0;
459 if (globalComm->SubDomainComm().MyPID()==0) factorToZeroOutCopies = 1.0;
461 (*g_out).Scale(factorToZeroOutCopies);
462 double* vPtr = &((*g_out)[0]);
464 globalComm->SumAll( &(tmp[0]), vPtr, num_g0);
470 globalComm->SumAll(tmp[0], mvPtr, num_dg0dp0);
Simple aggregate class for a derivative object represented as a column-wise multi-vector or its trans...
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Teuchos::RCP< Epetra_MultiVector > getMultiVector() const
void setSupports(EInArgsMembers arg, bool supports=true)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
void set_x(const Teuchos::RCP< const Epetra_Vector > &x)
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
void set_p(int l, const Teuchos::RCP< const Epetra_Vector > &p_l)
void setModelEvalDescription(const std::string &modelEvalDescription)
void set_W_properties(const DerivativeProperties &properties)
void set_DfDp_properties(int l, const DerivativeProperties &properties)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
void set_DgDx_properties(int j, const DerivativeProperties &properties)
void set_Np_Ng(int Np, int Ng)
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > get_W() const
void set_g(int j, const Evaluation< Epetra_Vector > &g_j)
Set g(j) where 0 <= j && j < this->Ng().
bool supports(EOutArgsMembers arg) const
Derivative get_DfDp(int l) const
Derivative get_DgDx(int j) const
void set_DfDp(int l, const Derivative &DfDp_l)
void set_DgDp(int j, int l, const Derivative &DgDp_j_l)
void set_DgDx(int j, const Derivative &DgDx_j)
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
Evaluation< Epetra_Vector > get_f() const
void set_f(const Evaluation< Epetra_Vector > &f)
Derivative get_DgDp(int j, int l) const
void set_W(const Teuchos::RCP< Epetra_Operator > &W)
EDerivativeMultiVectorOrientation orientation_DgDp
Some local data.
Teuchos::RCP< EpetraExt::MultiComm > globalComm
Pointer to the global (full XYZT) communicator.
std::vector< std::vector< long long > > * rowStencil_LL
Teuchos::RCP< const Epetra_Map > get_x_map() const
EpetraExt::BlockVector * block_x
Pointer to global multipoint solution vector – local storage.
EpetraExt::ModelEvaluator::DerivativeMultiVector * derivMV_DfDp
Teuchos::RCP< Epetra_RowMatrix > split_W
Pointer to split (spatial) Jacobian matrix.
EpetraExt::ModelEvaluator::Derivative * deriv_DfDp
Teuchos::RCP< Epetra_Vector > split_f
Split (spatial) residual vector – local storage.
Teuchos::RCP< const Epetra_Map > get_f_map() const
std::vector< int > * rowIndex_int
Set of indices into global XYZT Jacobian matrix.
~MultiPointModelEvaluator()
Teuchos::RCP< EpetraExt::BlockVector > solution_init
Pointer to initial multipoint solution vector.
OutArgs createOutArgs() const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< Epetra_Vector > split_x
Split (spatial) input vector – local storage.
int timeStepsOnTimeDomain
Number of time steps computed on each time domain.
Teuchos::RCP< Epetra_MultiVector > split_DfDp
Split sensitivity vector – local storage.
EpetraExt::ModelEvaluator::Derivative * deriv_DgDp
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
\breif .
EpetraExt::BlockVector * block_f
Pointer to global multipoint residual vector – local storage.
Teuchos::RCP< std::vector< Teuchos::RCP< Epetra_Vector > > > matching_vec
Array of vectors that have data for g-matching optimization problem.
EpetraExt::BlockMultiVector * block_DgDx
Pointer to global multipoint DfDp multi vector – local storage.
Teuchos::RCP< Epetra_MultiVector > split_DgDp
std::vector< long long > * rowIndex_LL
Teuchos::RCP< EpetraExt::ModelEvaluator > underlyingME
InArgs createInArgs() const
EpetraExt::ModelEvaluator::DerivativeMultiVector * derivMV_DgDx
Teuchos::RCP< Epetra_MultiVector > split_DgDx
Split sensitivity vector – local storage.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
MultiPointModelEvaluator(Teuchos::RCP< EpetraExt::ModelEvaluator > underlyingME_, const Teuchos::RCP< EpetraExt::MultiComm > &globalComm_, const std::vector< Epetra_Vector * > initGuessVec, Teuchos::RCP< std::vector< Teuchos::RCP< Epetra_Vector > > > qvec, Teuchos::RCP< std::vector< Teuchos::RCP< Epetra_Vector > > > matching_vec=Teuchos::null)
EpetraExt::ModelEvaluator::DerivativeMultiVector * derivMV_DgDp
EpetraExt::ModelEvaluator::Derivative * deriv_DgDx
std::vector< std::vector< int > > * rowStencil_int
Stencil for each row of global XYZT Jacobian matrix.
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< Epetra_Vector > split_g
Split vector of response functions – local storage.
int numTimeDomains
Total number of time step domains.
Teuchos::RCP< EpetraExt::BlockCrsMatrix > block_W
Pointer to global XYZT Jacobian matrix.
EpetraExt::BlockMultiVector * block_DfDp
Pointer to global multipoint DfDp multi vector – local storage.
int underlyingNg
Number of g vectors supported by underlyingME, often used as a bool.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
\breif .
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const