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);
285 outArgs.set_Np_Ng(1, underlyingNg);
286 outArgs.setSupports(OUT_ARG_f,
true);
287 outArgs.setSupports(OUT_ARG_W,
true);
288 outArgs.set_W_properties(
290 DERIV_LINEARITY_NONCONST
295 outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
296 outArgs.set_DfDp_properties(
298 DERIV_LINEARITY_CONST
299 ,DERIV_RANK_DEFICIENT
305 outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
306 outArgs.set_DgDx_properties(
308 DERIV_LINEARITY_NONCONST
309 ,DERIV_RANK_DEFICIENT
313 outArgs.setSupports(OUT_ARG_DgDp,0,0, orientation_DgDp);
314 outArgs.set_DgDp_properties(
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);
361 DgDx_out = outArgs.get_DgDx(0);
362 DgDp_out = outArgs.get_DgDp(0,0);
363 if (!DgDx_out.isEmpty()) DgDx_out.
getMultiVector()->PutScalar(0.0);
364 if (!DgDp_out.isEmpty()) DgDp_out.getMultiVector()->PutScalar(0.0);
369 bool need_g = g_out.get();
371 if ( !DgDx_out.isEmpty() || !DgDp_out.isEmpty() ) need_g =
true;
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);
400 if (!DfDp_out.isEmpty()) underlyingOutArgs.set_DfDp(0, *deriv_DfDp);
402 if (!DgDx_out.isEmpty()) underlyingOutArgs.set_DgDx(0, *deriv_DgDx);
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);
444 if (!DgDp_out.isEmpty())
445 DgDp_out.getMultiVector()->Update(1.0, *split_DgDp, 1.0);
450 if (f_out.get()) f_out->operator=(*block_f);
451 if (!DfDp_out.isEmpty())
452 DfDp_out.getMultiVector()->operator=(*block_DfDp);
453 if (!DgDx_out.isEmpty())
454 DgDx_out.getMultiVector()->operator=(*block_DgDx);
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);
466 if (!DgDp_out.isEmpty()) {
467 DgDp_out.getMultiVector()->
Scale(factorToZeroOutCopies);
468 double* mvPtr = (*DgDp_out.getMultiVector())[0];
470 globalComm->SumAll(tmp[0], mvPtr, num_dg0dp0);