Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_EpetraModelEvaluator.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include "Thyra_EpetraModelEvaluator.hpp"
43 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
45 #include "Thyra_EpetraLinearOp.hpp"
46 #include "Thyra_DetachedMultiVectorView.hpp"
47 #include "Thyra_ModelEvaluatorDelegatorBase.hpp" // Gives verbose macros!
48 #include "EpetraExt_ModelEvaluatorScalingTools.h"
49 #include "Epetra_RowMatrix.h"
50 #include "Teuchos_Time.hpp"
51 #include "Teuchos_implicit_cast.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include "Teuchos_StandardParameterEntryValidators.hpp"
54 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
55 
56 
57 namespace {
58 
59 
60 const std::string StateFunctionScaling_name = "State Function Scaling";
63  Thyra::EpetraModelEvaluator::EStateFunctionScaling
64  >
65  >
66 stateFunctionScalingValidator;
67 const std::string StateFunctionScaling_default = "None";
68 
69 
70 // Extract out the Epetra_RowMatrix from the set W in an Epetra OutArgs object
72 get_Epetra_RowMatrix(
73  const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs
74  )
75 {
76  using Teuchos::RCP;
77  const RCP<Epetra_Operator>
78  eW = epetraOutArgs.get_W();
79  const RCP<Epetra_RowMatrix>
80  ermW = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(eW,false);
82  is_null(ermW), std::logic_error,
83  "Thyra::EpetraModelEvaluator::evalModel(...): Error, if\n"
84  "scaling is turned on, the underlying Epetra_Operator created\n"
85  "an initialized by the underlying epetra model evaluator\n"
86  "\"" << epetraOutArgs.modelEvalDescription() << "\"\n"
87  "must support the Epetra_RowMatrix interface through a dynamic cast.\n"
88  "The concrete type " << Teuchos::typeName(*eW) << " does not support\n"
89  "Epetra_RowMatrix!"
90  );
91  return ermW;
92 }
93 
94 
96 create_and_assert_W(
97  const EpetraExt::ModelEvaluator &epetraModel
98  )
99 {
100  using Teuchos::RCP;
101  RCP<Epetra_Operator>
102  eW = epetraModel.create_W();
104  is_null(eW), std::logic_error,
105  "Error, the call to create_W() returned null on the "
106  "EpetraExt::ModelEvaluator object "
107  "\"" << epetraModel.description() << "\". This may mean that "
108  "the underlying model does not support more than one copy of "
109  "W at one time!" );
110  return eW;
111 }
112 
113 
114 } // namespace
115 
116 
117 namespace Thyra {
118 
119 
120 // Constructors/initializers/accessors.
121 
122 
124  :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
125  currentInArgsOutArgs_(false), finalPointWasSolved_(false)
126 {}
127 
128 
130  const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
131  const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
132  )
133  :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
134  currentInArgsOutArgs_(false), finalPointWasSolved_(false)
135 {
136  initialize(epetraModel,W_factory);
137 }
138 
139 
141  const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
142  const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
143  )
144 {
146  //typedef ModelEvaluatorBase MEB; // unused
147  //
148  epetraModel_ = epetraModel;
149  //
150  W_factory_ = W_factory;
151  //
152  x_map_ = epetraModel_->get_x_map();
153  f_map_ = epetraModel_->get_f_map();
154  if (!is_null(x_map_)) {
155  x_space_ = create_VectorSpace(x_map_);
156  f_space_ = create_VectorSpace(f_map_);
157  }
158  //
159  EpetraExt::ModelEvaluator::InArgs inArgs = epetraModel_->createInArgs();
160  p_map_.resize(inArgs.Np()); p_space_.resize(inArgs.Np());
161  p_map_is_local_.resize(inArgs.Np(),false);
162  for( int l = 0; l < implicit_cast<int>(p_space_.size()); ++l ) {
164  p_map_l = ( p_map_[l] = epetraModel_->get_p_map(l) );
165 #ifdef TEUCHOS_DEBUG
167  is_null(p_map_l), std::logic_error,
168  "Error, the the map p["<<l<<"] for the model \""
169  <<epetraModel->description()<<"\" can not be null!");
170 #endif
171 
172  p_map_is_local_[l] = !p_map_l->DistributedGlobal();
173  p_space_[l] = create_VectorSpace(p_map_l);
174  }
175  //
176  EpetraExt::ModelEvaluator::OutArgs outArgs = epetraModel_->createOutArgs();
177  g_map_.resize(outArgs.Ng()); g_space_.resize(outArgs.Ng());
178  g_map_is_local_.resize(outArgs.Ng(),false);
179  for( int j = 0; j < implicit_cast<int>(g_space_.size()); ++j ) {
181  g_map_j = ( g_map_[j] = epetraModel_->get_g_map(j) );
182  g_map_is_local_[j] = !g_map_j->DistributedGlobal();
183  g_space_[j] = create_VectorSpace( g_map_j );
184  }
185  //
186  epetraInArgsScaling_ = epetraModel_->createInArgs();
187  epetraOutArgsScaling_ = epetraModel_->createOutArgs();
188  nominalValuesAndBoundsAreUpdated_ = false;
189  finalPointWasSolved_ = false;
190  stateFunctionScalingVec_ = Teuchos::null; // Must set new scaling!
191  //
192  currentInArgsOutArgs_ = false;
193 }
194 
195 
198 {
199  return epetraModel_;
200 }
201 
202 
204  const ModelEvaluatorBase::InArgs<double>& nominalValues
205  )
206 {
207  nominalValues_.setArgs(nominalValues);
208  // Note: These must be the scaled values so we don't need to scale!
209 }
210 
211 
213  const RCP<const Epetra_Vector> &stateVariableScalingVec
214  )
215 {
216 #ifdef TEUCHOS_DEBUG
217  typedef ModelEvaluatorBase MEB;
218  TEUCHOS_TEST_FOR_EXCEPT( !this->createInArgs().supports(MEB::IN_ARG_x) );
219 #endif
220  stateVariableScalingVec_ = stateVariableScalingVec.assert_not_null();
221  invStateVariableScalingVec_ = Teuchos::null;
222  nominalValuesAndBoundsAreUpdated_ = false;
223 }
224 
225 
228 {
229  return stateVariableScalingVec_;
230 }
231 
232 
235 {
236  updateNominalValuesAndBounds();
237  return invStateVariableScalingVec_;
238 }
239 
240 
242  const RCP<const Epetra_Vector> &stateFunctionScalingVec
243  )
244 {
245  stateFunctionScalingVec_ = stateFunctionScalingVec;
246 }
247 
248 
251 {
252  return stateFunctionScalingVec_;
253 }
254 
255 
259  )
260 {
261  if(epetraModel) *epetraModel = epetraModel_;
262  if(W_factory) *W_factory = W_factory_;
263  epetraModel_ = Teuchos::null;
264  W_factory_ = Teuchos::null;
265  stateFunctionScalingVec_ = Teuchos::null;
266  stateVariableScalingVec_ = Teuchos::null;
267  invStateVariableScalingVec_ = Teuchos::null;
268  currentInArgsOutArgs_ = false;
269 }
270 
271 
274 {
275  return finalPoint_;
276 }
277 
278 
280 {
281  return finalPointWasSolved_;
282 }
283 
284 
285 // Public functions overridden from Teuchos::Describable
286 
287 
289 {
290  std::ostringstream oss;
291  oss << "Thyra::EpetraModelEvaluator{";
292  oss << "epetraModel=";
293  if(epetraModel_.get())
294  oss << "\'"<<epetraModel_->description()<<"\'";
295  else
296  oss << "NULL";
297  oss << ",W_factory=";
298  if(W_factory_.get())
299  oss << "\'"<<W_factory_->description()<<"\'";
300  else
301  oss << "NULL";
302  oss << "}";
303  return oss.str();
304 }
305 
306 
307 // Overridden from Teuchos::ParameterListAcceptor
308 
309 
311  RCP<Teuchos::ParameterList> const& paramList
312  )
313 {
314  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
315  paramList->validateParameters(*getValidParameters(),0); // Just validate my params
316  paramList_ = paramList;
317  const EStateFunctionScaling stateFunctionScaling_old = stateFunctionScaling_;
318  stateFunctionScaling_ = stateFunctionScalingValidator->getIntegralValue(
319  *paramList_, StateFunctionScaling_name, StateFunctionScaling_default
320  );
321  if( stateFunctionScaling_ != stateFunctionScaling_old )
322  stateFunctionScalingVec_ = Teuchos::null;
323  Teuchos::readVerboseObjectSublist(&*paramList_,this);
324 #ifdef TEUCHOS_DEBUG
325  paramList_->validateParameters(*getValidParameters(),0);
326 #endif // TEUCHOS_DEBUG
327 }
328 
329 
332 {
333  return paramList_;
334 }
335 
336 
339 {
340  RCP<Teuchos::ParameterList> _paramList = paramList_;
341  paramList_ = Teuchos::null;
342  return _paramList;
343 }
344 
345 
348 {
349  return paramList_;
350 }
351 
352 
355 {
356  using Teuchos::rcp;
358  using Teuchos::tuple;
359  using Teuchos::rcp_implicit_cast;
361  static RCP<const Teuchos::ParameterList> validPL;
362  if(is_null(validPL)) {
365  stateFunctionScalingValidator = rcp(
366  new StringToIntegralParameterEntryValidator<EStateFunctionScaling>(
367  tuple<std::string>(
368  "None",
369  "Row Sum"
370  ),
371  tuple<std::string>(
372  "Do not scale the state function f(...) in this class.",
373 
374  "Scale the state function f(...) and all its derivatives\n"
375  "using the row sum scaling from the initial Jacobian\n"
376  "W=d(f)/d(x). Note, this only works with Epetra_CrsMatrix\n"
377  "currently."
378  ),
379  tuple<EStateFunctionScaling>(
380  STATE_FUNC_SCALING_NONE,
381  STATE_FUNC_SCALING_ROW_SUM
382  ),
383  StateFunctionScaling_name
384  )
385  );
386  pl->set(StateFunctionScaling_name,StateFunctionScaling_default,
387  "Determines if and how the state function f(...) and all of its\n"
388  "derivatives are scaled. The scaling is done explicitly so there should\n"
389  "be no impact on the meaning of inner products or tolerances for\n"
390  "linear solves.",
391  rcp_implicit_cast<const PEV>(stateFunctionScalingValidator)
392  );
393  Teuchos::setupVerboseObjectSublist(&*pl);
394  validPL = pl;
395  }
396  return validPL;
397 }
398 
399 
400 // Overridden from ModelEvaulator.
401 
402 
404 {
405  return p_space_.size();
406 }
407 
408 
410 {
411  return g_space_.size();
412 }
413 
414 
417 {
418  return x_space_;
419 }
420 
421 
424 {
425  return f_space_;
426 }
427 
428 
431 {
432 #ifdef TEUCHOS_DEBUG
434 #endif
435  return p_space_[l];
436 }
437 
438 
441 {
442 #ifdef TEUCHOS_DEBUG
444 #endif
445  return epetraModel_->get_p_names(l);
446 }
447 
448 
451 {
452  TEUCHOS_TEST_FOR_EXCEPT( ! ( 0 <= j && j < this->Ng() ) );
453  return g_space_[j];
454 }
455 
456 
459 {
460 #ifdef TEUCHOS_DEBUG
462 #endif
463  return epetraModel_->get_g_names(j);
464 }
465 
466 
469 {
470  updateNominalValuesAndBounds();
471  return nominalValues_;
472 }
473 
474 
477 {
478  updateNominalValuesAndBounds();
479  return lowerBounds_;
480 }
481 
482 
485 {
486  updateNominalValuesAndBounds();
487  return upperBounds_;
488 }
489 
490 
493 {
494  return this->create_epetra_W_op();
495 }
496 
497 
500 {
501  return Teuchos::null;
502 }
503 
504 
507 {
508  return W_factory_;
509 }
510 
511 
513 {
514  if (!currentInArgsOutArgs_)
515  updateInArgsOutArgs();
516  return prototypeInArgs_;
517 }
518 
519 
521  const ModelEvaluatorBase::InArgs<double> &finalPoint,
522  const bool wasSolved
523  )
524 {
525  finalPoint_ = this->createInArgs();
526  finalPoint_.setArgs(finalPoint);
527  finalPointWasSolved_ = wasSolved;
528 }
529 
530 
531 // Private functions overridden from ModelEvaulatorDefaultBase
532 
533 
535 EpetraModelEvaluator::create_DfDp_op_impl(int /* l */) const
536 {
538  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
539 }
540 
541 
542 RCP<LinearOpBase<double> >
543 EpetraModelEvaluator::create_DgDx_dot_op_impl(int /* j */) const
544 {
546  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
547 }
548 
549 
550 RCP<LinearOpBase<double> >
551 EpetraModelEvaluator::create_DgDx_op_impl(int /* j */) const
552 {
554  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
555 }
556 
557 
558 RCP<LinearOpBase<double> >
559 EpetraModelEvaluator::create_DgDp_op_impl( int /* j */, int /* l */ ) const
560 {
562  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
563 }
564 
565 
566 ModelEvaluatorBase::OutArgs<double>
567 EpetraModelEvaluator::createOutArgsImpl() const
568 {
569  if (!currentInArgsOutArgs_) {
570  updateInArgsOutArgs();
571  }
572  return prototypeOutArgs_;
573 }
574 
575 
576 void EpetraModelEvaluator::evalModelImpl(
577  const ModelEvaluatorBase::InArgs<double>& inArgs_in,
578  const ModelEvaluatorBase::OutArgs<double>& outArgs
579  ) const
580 {
581 
582  using Teuchos::rcp;
583  using Teuchos::rcp_const_cast;
584  using Teuchos::rcp_dynamic_cast;
585  using Teuchos::OSTab;
587  typedef EpetraExt::ModelEvaluator EME;
588 
589  //
590  // A) Initial setup
591  //
592 
593  // Make sure that we are fully initialized!
594  this->updateNominalValuesAndBounds();
595 
596  // Make sure we grab the initial guess first!
597  InArgs<double> inArgs = this->getNominalValues();
598  // Now, copy the parameters from the input inArgs_in object to the inArgs
599  // object. Any input objects that are set in inArgs_in will overwrite those
600  // in inArgs that will already contain the nominal values. This will insure
601  // that all input parameters are set and those that are not set by the
602  // client will be at their nominal values (as defined by the underlying
603  // EpetraExt::ModelEvaluator object). The full set of Thyra input arguments
604  // must be set before these can be translated into Epetra input arguments.
605  inArgs.setArgs(inArgs_in);
606 
607  // This is a special exception: see evalModel() in Thyra::ME
608  // documentation. If inArgs() supports x_dot but the evaluate call
609  // passes in a null value, then we need to make sure the null value
610  // gets passed on instead of the nominal value.
611  if (inArgs.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) {
612  if (is_null(inArgs_in.get_x_dot()))
613  inArgs.set_x_dot(Teuchos::null);
614  }
615 
616  // Print the header and the values of the inArgs and outArgs objects!
617  typedef double Scalar; // Needed for below macro!
618  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
619  "Thyra::EpetraModelEvaluator",inArgs,outArgs,Teuchos::null
620  );
621 
622  // State function Scaling
623  const bool firstTimeStateFuncScaling
624  = (
625  stateFunctionScaling_ != STATE_FUNC_SCALING_NONE
626  && is_null(stateFunctionScalingVec_)
627  );
628 
630  VOTSLOWSF W_factory_outputTempState(W_factory_,out,verbLevel);
631 
632  Teuchos::Time timer("");
633 
634  //
635  // B) Prepressess the InArgs and OutArgs in preparation to call
636  // the underlying EpetraExt::ModelEvaluator
637  //
638 
639  //
640  // B.1) Translate InArgs from Thyra to Epetra objects
641  //
642 
643  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
644  *out << "\nSetting-up/creating input arguments ...\n";
645  timer.start(true);
646 
647  // Unwrap input Thyra objects to get Epetra objects
648  EME::InArgs epetraScaledInArgs = epetraModel_->createInArgs();
649  convertInArgsFromThyraToEpetra( inArgs, &epetraScaledInArgs );
650 
651  // Unscale the input Epetra objects which will be passed to the underlying
652  // EpetraExt::ModelEvaluator object.
653  EME::InArgs epetraInArgs = epetraModel_->createInArgs();
654  EpetraExt::unscaleModelVars(
655  epetraScaledInArgs, epetraInArgsScaling_, &epetraInArgs,
656  out.get(), verbLevel
657  );
658 
659  timer.stop();
660  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
661  OSTab(out).o() << "\nTime to setup InArgs = "<<timer.totalElapsedTime()<<" sec\n";
662 
663  //
664  // B.2) Convert from Thyra to Epetra OutArgs
665  //
666 
667  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
668  *out << "\nSetting-up/creating output arguments ...\n";
669  timer.start(true);
670 
671  // The unscaled Epetra OutArgs that will be passed to the
672  // underlying EpetraExt::ModelEvaluator object
673  EME::OutArgs epetraUnscaledOutArgs = epetraModel_->createOutArgs();
674 
675  // Various objects that are needed later (see documentation in
676  // the function convertOutArgsFromThyraToEpetra(...)
677  RCP<LinearOpBase<double> > W_op;
678  RCP<EpetraLinearOp> efwdW;
679  RCP<Epetra_Operator> eW;
680 
681  // Convert from Thyra to Epetra OutArgs and grap some of the intermediate
682  // objects accessed along the way that are needed later.
683  convertOutArgsFromThyraToEpetra(
684  outArgs,
685  &epetraUnscaledOutArgs,
686  &W_op, &efwdW, &eW
687  );
688 
689  //
690  // B.3) Setup OutArgs to computing scaling if needed
691  //
692 
693  if (firstTimeStateFuncScaling) {
694  preEvalScalingSetup(&epetraInArgs,&epetraUnscaledOutArgs,out,verbLevel);
695  }
696 
697  timer.stop();
698  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
699  OSTab(out).o()
700  << "\nTime to setup OutArgs = "
701  << timer.totalElapsedTime() <<" sec\n";
702 
703  //
704  // C) Evaluate the underlying EpetraExt model to compute the Epetra outputs
705  //
706 
707  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
708  *out << "\nEvaluating the Epetra output functions ...\n";
709  timer.start(true);
710 
711  epetraModel_->evalModel(epetraInArgs,epetraUnscaledOutArgs);
712 
713  timer.stop();
714  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
715  OSTab(out).o()
716  << "\nTime to evaluate Epetra output functions = "
717  << timer.totalElapsedTime() <<" sec\n";
718 
719  //
720  // D) Postprocess the output objects
721  //
722 
723  //
724  // D.1) Compute the scaling factors if needed
725  //
726 
727  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
728  *out << "\nCompute scale factors if needed ...\n";
729  timer.start(true);
730 
731  if (firstTimeStateFuncScaling) {
732  postEvalScalingSetup(epetraUnscaledOutArgs,out,verbLevel);
733  }
734 
735  timer.stop();
736  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
737  OSTab(out).o()
738  << "\nTime to compute scale factors = "
739  << timer.totalElapsedTime() <<" sec\n";
740 
741  //
742  // D.2) Scale the output Epetra objects
743  //
744 
745  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
746  *out << "\nScale the output objects ...\n";
747  timer.start(true);
748 
749  EME::OutArgs epetraOutArgs = epetraModel_->createOutArgs();
750  bool allFuncsWhereScaled = false;
751  EpetraExt::scaleModelFuncs(
752  epetraUnscaledOutArgs, epetraInArgsScaling_, epetraOutArgsScaling_,
753  &epetraOutArgs, &allFuncsWhereScaled,
754  out.get(), verbLevel
755  );
756  // AGS: This test precluded use of matrix-free Operators (vs. RowMatrix)
757  if (stateFunctionScaling_ != STATE_FUNC_SCALING_NONE)
759  !allFuncsWhereScaled, std::logic_error,
760  "Error, we can not currently handle epetra output objects that could not be"
761  " scaled. Special code will have to be added to handle this (i.e. using"
762  " implicit diagonal and multiplied linear operators to implicitly do"
763  " the scaling."
764  );
765 
766  timer.stop();
767  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
768  OSTab(out).o()
769  << "\nTime to scale the output objects = "
770  << timer.totalElapsedTime() << " sec\n";
771 
772  //
773  // D.3) Convert any Epetra objects to Thyra OutArgs objects that still need to
774  // be converted
775  //
776 
777  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
778  *out << "\nFinish processing and wrapping the output objects ...\n";
779  timer.start(true);
780 
781  finishConvertingOutArgsFromEpetraToThyra(
782  epetraOutArgs, W_op, efwdW, eW,
783  outArgs
784  );
785 
786  timer.stop();
787  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
788  OSTab(out).o()
789  << "\nTime to finish processing and wrapping the output objects = "
790  << timer.totalElapsedTime() <<" sec\n";
791 
792  //
793  // E) Print footer to end the function
794  //
795 
796  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
797 
798 }
799 
800 
801 // private
802 
803 
804 void EpetraModelEvaluator::convertInArgsFromEpetraToThyra(
805  const EpetraExt::ModelEvaluator::InArgs &epetraInArgs,
806  ModelEvaluatorBase::InArgs<double> *inArgs
807  ) const
808 {
809 
811  typedef ModelEvaluatorBase MEB;
812 
813  TEUCHOS_TEST_FOR_EXCEPT(!inArgs);
814 
815  if(inArgs->supports(MEB::IN_ARG_x)) {
816  inArgs->set_x( create_Vector( epetraInArgs.get_x(), x_space_ ) );
817  }
818 
819  if(inArgs->supports(MEB::IN_ARG_x_dot)) {
820  inArgs->set_x_dot( create_Vector( epetraInArgs.get_x_dot(), x_space_ ) );
821  }
822 
823  if(inArgs->supports(MEB::IN_ARG_x_mp)) {
824  inArgs->set_x_mp( epetraInArgs.get_x_mp() );
825  }
826 
827  if(inArgs->supports(MEB::IN_ARG_x_dot_mp)) {
828  inArgs->set_x_dot_mp( epetraInArgs.get_x_dot_mp() );
829  }
830 
831  const int l_Np = inArgs->Np();
832  for( int l = 0; l < l_Np; ++l ) {
833  inArgs->set_p( l, create_Vector( epetraInArgs.get_p(l), p_space_[l] ) );
834  }
835  for( int l = 0; l < l_Np; ++l ) {
836  if(inArgs->supports(MEB::IN_ARG_p_mp, l))
837  inArgs->set_p_mp( l, epetraInArgs.get_p_mp(l) );
838  }
839 
840  if(inArgs->supports(MEB::IN_ARG_t)) {
841  inArgs->set_t(epetraInArgs.get_t());
842  }
843 
844 }
845 
846 
847 void EpetraModelEvaluator::convertInArgsFromThyraToEpetra(
848  const ModelEvaluatorBase::InArgs<double> &inArgs,
849  EpetraExt::ModelEvaluator::InArgs *epetraInArgs
850  ) const
851 {
852 
853  using Teuchos::rcp;
854  using Teuchos::rcp_const_cast;
855 #ifdef HAVE_THYRA_ME_POLYNOMIAL
856  using Teuchos::Polynomial;
857 #endif // HAVE_THYRA_ME_POLYNOMIAL
858 
859 
860  TEUCHOS_TEST_FOR_EXCEPT(0==epetraInArgs);
861 
862  RCP<const VectorBase<double> > x_dot;
863  if( inArgs.supports(IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() ) {
864  RCP<const Epetra_Vector> e_x_dot = get_Epetra_Vector(*x_map_,x_dot);
865  epetraInArgs->set_x_dot(e_x_dot);
866  }
867  RCP<const Stokhos::ProductEpetraVector > x_dot_mp;
868  if( inArgs.supports(IN_ARG_x_dot_mp) && (x_dot_mp = inArgs.get_x_dot_mp()).get() ) {
869  epetraInArgs->set_x_dot_mp(x_dot_mp);
870  }
871 
872  RCP<const VectorBase<double> > x;
873  if( inArgs.supports(IN_ARG_x) && (x = inArgs.get_x()).get() ) {
874  RCP<const Epetra_Vector> e_x = get_Epetra_Vector(*x_map_,x);
875  epetraInArgs->set_x(e_x);
876  }
877  RCP<const Stokhos::ProductEpetraVector > x_mp;
878  if( inArgs.supports(IN_ARG_x_mp) && (x_mp = inArgs.get_x_mp()).get() ) {
879  epetraInArgs->set_x_mp(x_mp);
880  }
881 
882  RCP<const VectorBase<double> > p_l;
883  for(int l = 0; l < inArgs.Np(); ++l ) {
884  p_l = inArgs.get_p(l);
885  if(p_l.get()) epetraInArgs->set_p(l,get_Epetra_Vector(*p_map_[l],p_l));
886  }
887  RCP<const Stokhos::ProductEpetraVector > p_mp_l;
888  for(int l = 0; l < inArgs.Np(); ++l ) {
889  if (inArgs.supports(IN_ARG_p_mp,l)) {
890  p_mp_l = inArgs.get_p_mp(l);
891  if(p_mp_l.get()) epetraInArgs->set_p_mp(l,p_mp_l);
892  }
893  }
894 
895 #ifdef HAVE_THYRA_ME_POLYNOMIAL
896 
897  RCP<const Polynomial< VectorBase<double> > > x_dot_poly;
898  RCP<Epetra_Vector> epetra_ptr;
899  if(
900  inArgs.supports(IN_ARG_x_dot_poly)
901  && (x_dot_poly = inArgs.get_x_dot_poly()).get()
902  )
903  {
904  RCP<Polynomial<Epetra_Vector> > epetra_x_dot_poly =
905  rcp(new Polynomial<Epetra_Vector>(x_dot_poly->degree()));
906  for (unsigned int i=0; i<=x_dot_poly->degree(); i++) {
907  epetra_ptr = rcp_const_cast<Epetra_Vector>(
908  get_Epetra_Vector(*x_map_, x_dot_poly->getCoefficient(i)) );
909  epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
910  }
911  epetraInArgs->set_x_dot_poly(epetra_x_dot_poly);
912  }
913 
914  RCP<const Polynomial< VectorBase<double> > > x_poly;
915  if(
916  inArgs.supports(IN_ARG_x_poly)
917  && (x_poly = inArgs.get_x_poly()).get()
918  )
919  {
920  RCP<Polynomial<Epetra_Vector> > epetra_x_poly =
921  rcp(new Polynomial<Epetra_Vector>(x_poly->degree()));
922  for (unsigned int i=0; i<=x_poly->degree(); i++) {
923  epetra_ptr = rcp_const_cast<Epetra_Vector>(
924  get_Epetra_Vector(*x_map_, x_poly->getCoefficient(i)) );
925  epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
926  }
927  epetraInArgs->set_x_poly(epetra_x_poly);
928  }
929 
930 #endif // HAVE_THYRA_ME_POLYNOMIAL
931 
932  if( inArgs.supports(IN_ARG_t) )
933  epetraInArgs->set_t(inArgs.get_t());
934 
935  if( inArgs.supports(IN_ARG_alpha) )
936  epetraInArgs->set_alpha(inArgs.get_alpha());
937 
938  if( inArgs.supports(IN_ARG_beta) )
939  epetraInArgs->set_beta(inArgs.get_beta());
940 
941  if( inArgs.supports(IN_ARG_step_size) )
942  epetraInArgs->set_step_size(inArgs.get_step_size());
943 
944  if( inArgs.supports(IN_ARG_stage_number) )
945  epetraInArgs->set_stage_number(inArgs.get_stage_number());
946 }
947 
948 
949 void EpetraModelEvaluator::convertOutArgsFromThyraToEpetra(
950  const ModelEvaluatorBase::OutArgs<double> &outArgs,
951  EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
952  RCP<LinearOpBase<double> > *W_op_out,
953  RCP<EpetraLinearOp> *efwdW_out,
954  RCP<Epetra_Operator> *eW_out
955  ) const
956 {
957 
958  using Teuchos::rcp;
959  using Teuchos::rcp_const_cast;
960  using Teuchos::rcp_dynamic_cast;
961  using Teuchos::OSTab;
964  //typedef EpetraExt::ModelEvaluator EME; // unused
965 
966  // Assert input
967 #ifdef TEUCHOS_DEBUG
968  TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
969  TEUCHOS_ASSERT(W_op_out);
970  TEUCHOS_ASSERT(efwdW_out);
971  TEUCHOS_ASSERT(eW_out);
972 #endif
973 
974  // Create easy to use references
975  EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
976  RCP<LinearOpBase<double> > &W_op = *W_op_out;
977  RCP<EpetraLinearOp> &efwdW = *efwdW_out;
978  RCP<Epetra_Operator> &eW = *eW_out;
979 
980  // f
981  {
982  RCP<VectorBase<double> > f;
983  if( outArgs.supports(OUT_ARG_f) && (f = outArgs.get_f()).get() )
984  epetraUnscaledOutArgs.set_f(get_Epetra_Vector(*f_map_,f));
985  RCP<Stokhos::ProductEpetraVector > f_mp;
986  if( outArgs.supports(OUT_ARG_f_mp) && (f_mp = outArgs.get_f_mp()).get() )
987  epetraUnscaledOutArgs.set_f_mp(f_mp);
988  }
989 
990  // g
991  {
992  RCP<VectorBase<double> > g_j;
993  for(int j = 0; j < outArgs.Ng(); ++j ) {
994  g_j = outArgs.get_g(j);
995  if(g_j.get()) epetraUnscaledOutArgs.set_g(j,get_Epetra_Vector(*g_map_[j],g_j));
996  }
997  RCP<Stokhos::ProductEpetraVector > g_mp_j;
998  for(int j = 0; j < outArgs.Ng(); ++j ) {
999  if (outArgs.supports(OUT_ARG_g_mp,j)) {
1000  g_mp_j = outArgs.get_g_mp(j);
1001  if(g_mp_j.get()) epetraUnscaledOutArgs.set_g_mp(j,g_mp_j);
1002  }
1003  }
1004  }
1005 
1006  // W_op
1007  {
1008  RCP<Stokhos::ProductEpetraOperator > W_mp;
1009  if( outArgs.supports(OUT_ARG_W_mp) && (W_mp = outArgs.get_W_mp()).get() )
1010  epetraUnscaledOutArgs.set_W_mp(W_mp);
1011  }
1012 
1013  // W_op
1014  {
1015 
1016  if (outArgs.supports(OUT_ARG_W_op) && nonnull(W_op = outArgs.get_W_op())) {
1017  if (nonnull(W_op) && is_null(efwdW)) {
1018  efwdW = rcp_const_cast<EpetraLinearOp>(
1019  rcp_dynamic_cast<const EpetraLinearOp>(W_op, true));
1020  }
1021  }
1022 
1023  if (nonnull(efwdW)) {
1024  // By the time we get here, if we have an object in efwdW, then it
1025  // should already be embeadded with an underlying Epetra_Operator object
1026  // that was allocated by the EpetraExt::ModelEvaluator object.
1027  // Therefore, we should just have to grab this object and be on our way.
1028  eW = efwdW->epetra_op();
1029  epetraUnscaledOutArgs.set_W(eW);
1030  }
1031 
1032  // Note: The following derivative objects update in place!
1033 
1034  }
1035 
1036  // DfDp
1037  {
1038  Derivative<double> DfDp_l;
1039  for(int l = 0; l < outArgs.Np(); ++l ) {
1040  if( !outArgs.supports(OUT_ARG_DfDp,l).none()
1041  && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
1042  {
1043  epetraUnscaledOutArgs.set_DfDp(l,convert(DfDp_l,f_map_,p_map_[l]));
1044  }
1045  }
1046  MPDerivative DfDp_mp_l;
1047  for(int l = 0; l < outArgs.Np(); ++l ) {
1048  if( !outArgs.supports(OUT_ARG_DfDp_mp,l).none()
1049  && !(DfDp_mp_l = outArgs.get_DfDp_mp(l)).isEmpty() )
1050  {
1051  epetraUnscaledOutArgs.set_DfDp_mp(l,convert(DfDp_mp_l,f_map_,p_map_[l]));
1052  }
1053  }
1054  }
1055 
1056  // DgDx_dot
1057  {
1058  Derivative<double> DgDx_dot_j;
1059  for(int j = 0; j < outArgs.Ng(); ++j ) {
1060  if( !outArgs.supports(OUT_ARG_DgDx_dot,j).none()
1061  && !(DgDx_dot_j = outArgs.get_DgDx_dot(j)).isEmpty() )
1062  {
1063  epetraUnscaledOutArgs.set_DgDx_dot(j,convert(DgDx_dot_j,g_map_[j],x_map_));
1064  }
1065  }
1066  MPDerivative DgDx_dot_mp_j;
1067  for(int j = 0; j < outArgs.Ng(); ++j ) {
1068  if( !outArgs.supports(OUT_ARG_DgDx_dot_mp,j).none()
1069  && !(DgDx_dot_mp_j = outArgs.get_DgDx_dot_mp(j)).isEmpty() )
1070  {
1071  epetraUnscaledOutArgs.set_DgDx_dot_mp(j,convert(DgDx_dot_mp_j,g_map_[j],x_map_));
1072  }
1073  }
1074  }
1075 
1076  // DgDx
1077  {
1078  Derivative<double> DgDx_j;
1079  for(int j = 0; j < outArgs.Ng(); ++j ) {
1080  if( !outArgs.supports(OUT_ARG_DgDx,j).none()
1081  && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
1082  {
1083  epetraUnscaledOutArgs.set_DgDx(j,convert(DgDx_j,g_map_[j],x_map_));
1084  }
1085  }
1086  MPDerivative DgDx_mp_j;
1087  for(int j = 0; j < outArgs.Ng(); ++j ) {
1088  if( !outArgs.supports(OUT_ARG_DgDx_mp,j).none()
1089  && !(DgDx_mp_j = outArgs.get_DgDx_mp(j)).isEmpty() )
1090  {
1091  epetraUnscaledOutArgs.set_DgDx_mp(j,convert(DgDx_mp_j,g_map_[j],x_map_));
1092  }
1093  }
1094  }
1095 
1096  // DgDp
1097  {
1098  DerivativeSupport DgDp_j_l_support;
1099  Derivative<double> DgDp_j_l;
1100  for (int j = 0; j < outArgs.Ng(); ++j ) {
1101  for (int l = 0; l < outArgs.Np(); ++l ) {
1102  if (!(DgDp_j_l_support = outArgs.supports(OUT_ARG_DgDp,j,l)).none()
1103  && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
1104  {
1105  epetraUnscaledOutArgs.set_DgDp(j,l,convert(DgDp_j_l,g_map_[j],p_map_[l]));
1106  }
1107  }
1108  }
1109  DerivativeSupport DgDp_mp_j_l_support;
1110  MPDerivative DgDp_mp_j_l;
1111  for (int j = 0; j < outArgs.Ng(); ++j ) {
1112  for (int l = 0; l < outArgs.Np(); ++l ) {
1113  if (!(DgDp_mp_j_l_support = outArgs.supports(OUT_ARG_DgDp_mp,j,l)).none()
1114  && !(DgDp_mp_j_l = outArgs.get_DgDp_mp(j,l)).isEmpty() )
1115  {
1116  epetraUnscaledOutArgs.set_DgDp_mp(j,l,convert(DgDp_mp_j_l,g_map_[j],p_map_[l]));
1117  }
1118  }
1119  }
1120  }
1121 
1122 #ifdef HAVE_THYRA_ME_POLYNOMIAL
1123 
1124  // f_poly
1125  RCP<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
1126  if( outArgs.supports(OUT_ARG_f_poly) && (f_poly = outArgs.get_f_poly()).get() )
1127  {
1128  RCP<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly =
1129  Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(f_poly->degree()));
1130  for (unsigned int i=0; i<=f_poly->degree(); i++) {
1131  RCP<Epetra_Vector> epetra_ptr
1132  = Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*f_map_,
1133  f_poly->getCoefficient(i)));
1134  epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
1135  }
1136  epetraUnscaledOutArgs.set_f_poly(epetra_f_poly);
1137  }
1138 
1139 #endif // HAVE_THYRA_ME_POLYNOMIAL
1140 
1141 }
1142 
1143 
1144 void EpetraModelEvaluator::preEvalScalingSetup(
1145  EpetraExt::ModelEvaluator::InArgs *epetraInArgs_inout,
1146  EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
1147  const RCP<Teuchos::FancyOStream> &out,
1148  const Teuchos::EVerbosityLevel verbLevel
1149  ) const
1150 {
1151 
1152  typedef EpetraExt::ModelEvaluator EME;
1153 
1154 #ifdef TEUCHOS_DEBUG
1155  TEUCHOS_ASSERT(epetraInArgs_inout);
1156  TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
1157 #endif
1158 
1159  EpetraExt::ModelEvaluator::InArgs
1160  &epetraInArgs = *epetraInArgs_inout;
1161  EpetraExt::ModelEvaluator::OutArgs
1162  &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
1163 
1164  if (
1165  ( stateFunctionScaling_ == STATE_FUNC_SCALING_ROW_SUM )
1166  &&
1167  (
1168  epetraUnscaledOutArgs.supports(EME::OUT_ARG_f)
1169  &&
1170  epetraUnscaledOutArgs.funcOrDerivesAreSet(EME::OUT_ARG_f)
1171  )
1172  &&
1173  (
1174  epetraUnscaledOutArgs.supports(EME::OUT_ARG_W)
1175  &&
1176  is_null(epetraUnscaledOutArgs.get_W())
1177  )
1178  )
1179  {
1180  // This is the first pass through with scaling turned on and the client
1181  // turned on automatic scaling but did not ask for W. We must compute W
1182  // in order to compute the scale factors so we must allocate a temporary W
1183  // just to compute the scale factors and then throw it away. If the
1184  // client wants to evaluate W at the same point, then it should have
1185  // passed W in but that is not our problem here. The ModelEvaluator
1186  // relies on the client to set up the calls to allow for efficient
1187  // evaluation.
1188 
1189  if(out.get() && verbLevel >= Teuchos::VERB_LOW)
1190  *out
1191  << "\nCreating a temporary Epetra W to compute scale factors"
1192  << " for f(...) ...\n";
1193  epetraUnscaledOutArgs.set_W(create_and_assert_W(*epetraModel_));
1194  if( epetraInArgs.supports(EME::IN_ARG_beta) )
1195  epetraInArgs.set_beta(1.0);
1196  if( epetraInArgs.supports(EME::IN_ARG_alpha) )
1197  epetraInArgs.set_alpha(0.0);
1198  if( epetraInArgs.supports(EME::IN_ARG_step_size) )
1199  epetraInArgs.set_step_size(0.0);
1200  if( epetraInArgs.supports(EME::IN_ARG_stage_number) )
1201  epetraInArgs.set_stage_number(1.0);
1202  }
1203 
1204 }
1205 
1206 
1207 void EpetraModelEvaluator::postEvalScalingSetup(
1208  const EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs,
1209  const RCP<Teuchos::FancyOStream> &out,
1210  const Teuchos::EVerbosityLevel verbLevel
1211  ) const
1212 {
1213 
1214  using Teuchos::OSTab;
1215  using Teuchos::rcp;
1216  using Teuchos::rcp_const_cast;
1218 
1219  // Compute the scale factors for the state function f(...)
1220  switch(stateFunctionScaling_) {
1221 
1222  case STATE_FUNC_SCALING_ROW_SUM: {
1223 
1224  // Compute the inverse row-sum scaling from W
1225 
1226  const RCP<Epetra_RowMatrix>
1227  ermW = get_Epetra_RowMatrix(epetraUnscaledOutArgs);
1228  // Note: Above, we get the Epetra W object directly from the Epetra
1229  // OutArgs object since this might be a temporary matrix just to
1230  // compute scaling factors. In this case, the stack funtion variable
1231  // eW might be empty!
1232 
1233  RCP<Epetra_Vector>
1234  invRowSums = rcp(new Epetra_Vector(ermW->OperatorRangeMap()));
1235  // Above: From the documentation is seems that the RangeMap should be
1236  // okay but who knows for sure!
1237 
1238  ermW->InvRowSums(*invRowSums);
1239 
1240  if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW)) {
1241  *out
1242  << "\nComputed inverse row sum scaling from W that"
1243  " will be used to scale f(...) and its derivatives:\n";
1244  double minVal = 0, maxVal = 0, avgVal = 0;
1245  invRowSums->MinValue(&minVal);
1246  invRowSums->MaxValue(&maxVal);
1247  invRowSums->MeanValue(&avgVal);
1248  OSTab tab(out);
1249  *out
1250  << "min(invRowSums) = " << minVal << "\n"
1251  << "max(invRowSums) = " << maxVal << "\n"
1252  << "avg(invRowSums) = " << avgVal << "\n";
1253  }
1254 
1255  stateFunctionScalingVec_ = invRowSums;
1256 
1257  break;
1258 
1259  }
1260 
1261  default:
1262  TEUCHOS_TEST_FOR_EXCEPT("Should never get here!");
1263 
1264  }
1265 
1266  epetraOutArgsScaling_ = epetraModel_->createOutArgs();
1267 
1268  epetraOutArgsScaling_.set_f(
1269  rcp_const_cast<Epetra_Vector>(stateFunctionScalingVec_) );
1270 
1271 }
1272 
1273 
1274 void EpetraModelEvaluator::finishConvertingOutArgsFromEpetraToThyra(
1275  const EpetraExt::ModelEvaluator::OutArgs &/* epetraOutArgs */,
1276  RCP<LinearOpBase<double> > &W_op,
1277  RCP<EpetraLinearOp> &efwdW,
1278  RCP<Epetra_Operator> &/* eW */,
1279  const ModelEvaluatorBase::OutArgs<double> &/* outArgs */
1280  ) const
1281 {
1282 
1283  using Teuchos::rcp_dynamic_cast;
1284  //typedef EpetraExt::ModelEvaluator EME; // unused
1285 
1286  if (nonnull(efwdW)) {
1287  efwdW->setFullyInitialized(true);
1288  // NOTE: Above will directly update W_op also if W.get()==NULL!
1289  }
1290 
1291  if (nonnull(W_op)) {
1292  if (W_op.shares_resource(efwdW)) {
1293  // W_op was already updated above since *efwdW is the same object as *W_op
1294  }
1295  else {
1296  rcp_dynamic_cast<EpetraLinearOp>(W_op, true)->setFullyInitialized(true);
1297  }
1298  }
1299 
1300 }
1301 
1302 
1303 void EpetraModelEvaluator::updateNominalValuesAndBounds() const
1304 {
1305 
1306  using Teuchos::rcp;
1307  using Teuchos::implicit_cast;
1308  //typedef ModelEvaluatorBase MEB; // unused
1309  typedef EpetraExt::ModelEvaluator EME;
1310 
1311  if( !nominalValuesAndBoundsAreUpdated_ ) {
1312 
1313  // Gather the nominal values and bounds into Epetra InArgs objects
1314 
1315  EME::InArgs epetraOrigNominalValues;
1316  EpetraExt::gatherModelNominalValues(
1317  *epetraModel_, &epetraOrigNominalValues );
1318 
1319  EME::InArgs epetraOrigLowerBounds;
1320  EME::InArgs epetraOrigUpperBounds;
1321  EpetraExt::gatherModelBounds(
1322  *epetraModel_, &epetraOrigLowerBounds, &epetraOrigUpperBounds );
1323 
1324  // Set up Epetra InArgs scaling object
1325 
1326  epetraInArgsScaling_ = epetraModel_->createInArgs();
1327 
1328  if( !is_null(stateVariableScalingVec_) ) {
1329  invStateVariableScalingVec_
1330  = EpetraExt::createInverseModelScalingVector(stateVariableScalingVec_);
1331  if( epetraOrigNominalValues.supports(EME::IN_ARG_x_dot) ) {
1332  epetraInArgsScaling_.set_x_dot(invStateVariableScalingVec_);
1333  }
1334  if( epetraOrigNominalValues.supports(EME::IN_ARG_x) ) {
1335  epetraInArgsScaling_.set_x(invStateVariableScalingVec_);
1336  }
1337  }
1338 
1339  // Scale the original variables and bounds
1340 
1341  EME::InArgs epetraScaledNominalValues = epetraModel_->createInArgs();
1342  EpetraExt::scaleModelVars(
1343  epetraOrigNominalValues, epetraInArgsScaling_, &epetraScaledNominalValues
1344  );
1345 
1346  EME::InArgs epetraScaledLowerBounds = epetraModel_->createInArgs();
1347  EME::InArgs epetraScaledUpperBounds = epetraModel_->createInArgs();
1348  EpetraExt::scaleModelBounds(
1349  epetraOrigLowerBounds, epetraOrigUpperBounds, epetraModel_->getInfBound(),
1350  epetraInArgsScaling_,
1351  &epetraScaledLowerBounds, &epetraScaledUpperBounds
1352  );
1353 
1354  // Wrap the scaled epetra InArgs objects as Thyra InArgs objects!
1355 
1356  nominalValues_ = this->createInArgs();
1357  lowerBounds_ = this->createInArgs();
1358  upperBounds_ = this->createInArgs();
1359  convertInArgsFromEpetraToThyra(epetraScaledNominalValues, &nominalValues_);
1360  convertInArgsFromEpetraToThyra(epetraScaledLowerBounds, &lowerBounds_);
1361  convertInArgsFromEpetraToThyra(epetraScaledUpperBounds, &upperBounds_);
1362 
1363  nominalValuesAndBoundsAreUpdated_ = true;
1364 
1365  }
1366  else {
1367 
1368  // The nominal values and bounds should already be updated an should have
1369  // the currect scaling!
1370 
1371  }
1372 
1373 }
1374 
1375 
1376 void EpetraModelEvaluator::updateInArgsOutArgs() const
1377 {
1378 
1379  typedef EpetraExt::ModelEvaluator EME;
1380 
1381  const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
1382  EME::InArgs epetraInArgs = epetraModel.createInArgs();
1383  EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
1384  const int l_Np = epetraOutArgs.Np();
1385  const int l_Ng = epetraOutArgs.Ng();
1386 
1387  //
1388  // InArgs
1389  //
1390 
1391  InArgsSetup<double> inArgs;
1392  inArgs.setModelEvalDescription(this->description());
1393  inArgs.set_Np(epetraInArgs.Np());
1394  inArgs.setSupports(IN_ARG_x_dot, epetraInArgs.supports(EME::IN_ARG_x_dot));
1395  inArgs.setSupports(IN_ARG_x, epetraInArgs.supports(EME::IN_ARG_x));
1396  inArgs.setSupports(IN_ARG_x_dot_mp, epetraInArgs.supports(EME::IN_ARG_x_dot_mp));
1397  inArgs.setSupports(IN_ARG_x_mp, epetraInArgs.supports(EME::IN_ARG_x_mp));
1398 #ifdef HAVE_THYRA_ME_POLYNOMIAL
1399  inArgs.setSupports(IN_ARG_x_dot_poly,
1400  epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
1401  inArgs.setSupports(IN_ARG_x_poly, epetraInArgs.supports(EME::IN_ARG_x_poly));
1402 #endif // HAVE_THYRA_ME_POLYNOMIAL
1403  inArgs.setSupports(IN_ARG_t, epetraInArgs.supports(EME::IN_ARG_t));
1404  inArgs.setSupports(IN_ARG_alpha, epetraInArgs.supports(EME::IN_ARG_alpha));
1405  inArgs.setSupports(IN_ARG_beta, epetraInArgs.supports(EME::IN_ARG_beta));
1406  inArgs.setSupports(IN_ARG_step_size, epetraInArgs.supports(EME::IN_ARG_step_size));
1407  inArgs.setSupports(IN_ARG_stage_number, epetraInArgs.supports(EME::IN_ARG_stage_number));
1408  for(int l=0; l<l_Np; ++l) {
1409  inArgs.setSupports(IN_ARG_p_mp, l, epetraInArgs.supports(EME::IN_ARG_p_mp, l));
1410  }
1411  prototypeInArgs_ = inArgs;
1412 
1413  //
1414  // OutArgs
1415  //
1416 
1417  OutArgsSetup<double> outArgs;
1418  outArgs.setModelEvalDescription(this->description());
1419  outArgs.set_Np_Ng(l_Np, l_Ng);
1420  // f
1421  outArgs.setSupports(OUT_ARG_f, epetraOutArgs.supports(EME::OUT_ARG_f));
1422  outArgs.setSupports(OUT_ARG_f_mp, epetraOutArgs.supports(EME::OUT_ARG_f_mp));
1423  if (outArgs.supports(OUT_ARG_f)) {
1424  // W_op
1425  outArgs.setSupports(OUT_ARG_W_op, epetraOutArgs.supports(EME::OUT_ARG_W));
1426  outArgs.set_W_properties(convert(epetraOutArgs.get_W_properties()));
1427  // DfDp
1428  for(int l=0; l<l_Np; ++l) {
1429  outArgs.setSupports(OUT_ARG_DfDp, l,
1430  convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp, l)));
1431  if(!outArgs.supports(OUT_ARG_DfDp, l).none())
1432  outArgs.set_DfDp_properties(l,
1433  convert(epetraOutArgs.get_DfDp_properties(l)));
1434  if (outArgs.supports(OUT_ARG_f_mp))
1435  {
1436  outArgs.setSupports(OUT_ARG_DfDp_mp, l,
1437  convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp_mp, l)));
1438  if(!outArgs.supports(OUT_ARG_DfDp_mp, l).none())
1439  outArgs.set_DfDp_mp_properties(l,
1440  convert(epetraOutArgs.get_DfDp_mp_properties(l)));
1441  }
1442  }
1443  }
1444  // DgDx_dot and DgDx
1445  for(int j=0; j<l_Ng; ++j) {
1446  if (inArgs.supports(IN_ARG_x_dot))
1447  outArgs.setSupports(OUT_ARG_DgDx_dot, j,
1448  convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot, j)));
1449  if(!outArgs.supports(OUT_ARG_DgDx_dot, j).none())
1450  outArgs.set_DgDx_dot_properties(j,
1451  convert(epetraOutArgs.get_DgDx_dot_properties(j)));
1452  if (inArgs.supports(IN_ARG_x))
1453  outArgs.setSupports(OUT_ARG_DgDx, j,
1454  convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx, j)));
1455  if(!outArgs.supports(OUT_ARG_DgDx, j).none())
1456  outArgs.set_DgDx_properties(j,
1457  convert(epetraOutArgs.get_DgDx_properties(j)));
1458  if (inArgs.supports(IN_ARG_x_dot_mp))
1459  outArgs.setSupports(OUT_ARG_DgDx_dot_mp, j,
1460  convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot_mp, j)));
1461  if(!outArgs.supports(OUT_ARG_DgDx_dot_mp, j).none())
1462  outArgs.set_DgDx_dot_mp_properties(j,
1463  convert(epetraOutArgs.get_DgDx_dot_mp_properties(j)));
1464  if (inArgs.supports(IN_ARG_x_mp))
1465  outArgs.setSupports(OUT_ARG_DgDx_mp, j,
1466  convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_mp, j)));
1467  if(!outArgs.supports(OUT_ARG_DgDx_mp, j).none())
1468  outArgs.set_DgDx_mp_properties(j,
1469  convert(epetraOutArgs.get_DgDx_mp_properties(j)));
1470  }
1471  // DgDp
1472  for(int j=0; j < l_Ng; ++j) for(int l=0; l < l_Np; ++l) {
1473  const EME::DerivativeSupport epetra_DgDp_j_l_support =
1474  epetraOutArgs.supports(EME::OUT_ARG_DgDp, j, l);
1475  outArgs.setSupports(OUT_ARG_DgDp, j, l,
1476  convert(epetra_DgDp_j_l_support));
1477  if(!outArgs.supports(OUT_ARG_DgDp, j, l).none())
1478  outArgs.set_DgDp_properties(j, l,
1479  convert(epetraOutArgs.get_DgDp_properties(j, l)));
1480  const EME::DerivativeSupport epetra_DgDp_mp_j_l_support =
1481  epetraOutArgs.supports(EME::OUT_ARG_DgDp_mp, j, l);
1482  outArgs.setSupports(OUT_ARG_DgDp_mp, j, l,
1483  convert(epetra_DgDp_mp_j_l_support));
1484  if(!outArgs.supports(OUT_ARG_DgDp_mp, j, l).none())
1485  outArgs.set_DgDp_mp_properties(j, l,
1486  convert(epetraOutArgs.get_DgDp_mp_properties(j, l)));
1487  }
1488  for(int j=0; j<l_Ng; ++j) {
1489  outArgs.setSupports(OUT_ARG_g_mp, j, epetraOutArgs.supports(EME::OUT_ARG_g_mp, j));
1490  }
1491 #ifdef HAVE_THYRA_ME_POLYNOMIAL
1492  outArgs.setSupports(OUT_ARG_f_poly,
1493  epetraOutArgs.supports(EME::OUT_ARG_f_poly));
1494 #endif // HAVE_THYRA_ME_POLYNOMIAL
1495  prototypeOutArgs_ = outArgs;
1496 
1497  // We are current!
1498  currentInArgsOutArgs_ = true;
1499 
1500 }
1501 
1502 
1503 RCP<EpetraLinearOp>
1504 EpetraModelEvaluator::create_epetra_W_op() const
1505 {
1506  return Thyra::partialNonconstEpetraLinearOp(
1507  this->get_f_space(), this->get_x_space(),
1508  create_and_assert_W(*epetraModel_)
1509  );
1510 }
1511 
1512 
1513 } // namespace Thyra
1514 
1515 
1516 //
1517 // Non-member utility functions
1518 //
1519 
1520 
1522 Thyra::epetraModelEvaluator(
1523  const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
1524  const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
1525  )
1526 {
1527  return Teuchos::rcp(new EpetraModelEvaluator(epetraModel,W_factory));
1528 }
1529 
1530 
1532 Thyra::convert(
1533  const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation
1534  )
1535 {
1536  switch(mvOrientation) {
1537  case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
1538  return ModelEvaluatorBase::DERIV_MV_BY_COL;
1539  case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
1540  return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
1541  default:
1543  }
1544  TEUCHOS_UNREACHABLE_RETURN(ModelEvaluatorBase::DERIV_MV_BY_COL);
1545 }
1546 
1547 
1548 EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
1549 Thyra::convert(
1550  const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation
1551  )
1552 {
1553  switch(mvOrientation) {
1554  case ModelEvaluatorBase::DERIV_MV_BY_COL :
1555  return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
1556  case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
1557  return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
1558  default:
1560  }
1561  TEUCHOS_UNREACHABLE_RETURN(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL);
1562 }
1563 
1564 
1566 Thyra::convert(
1567  const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties
1568  )
1569 {
1570  ModelEvaluatorBase::EDerivativeLinearity linearity;
1571  switch(derivativeProperties.linearity) {
1572  case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
1573  linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
1574  break;
1575  case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
1576  linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
1577  break;
1578  case EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
1579  linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
1580  break;
1581  default:
1583  }
1584  ModelEvaluatorBase::ERankStatus rank;
1585  switch(derivativeProperties.rank) {
1586  case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
1587  rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
1588  break;
1589  case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
1590  rank = ModelEvaluatorBase::DERIV_RANK_FULL;
1591  break;
1592  case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
1593  rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
1594  break;
1595  default:
1597  }
1598  return ModelEvaluatorBase::DerivativeProperties(
1599  linearity,rank,derivativeProperties.supportsAdjoint);
1600 }
1601 
1602 
1604 Thyra::convert(
1605  const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport
1606  )
1607 {
1608  ModelEvaluatorBase::DerivativeSupport ds;
1609  if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
1610  ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
1611  if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
1612  ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
1613  if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
1614  ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
1615  return ds;
1616 }
1617 
1618 
1619 EpetraExt::ModelEvaluator::Derivative
1620 Thyra::convert(
1621  const ModelEvaluatorBase::Derivative<double> &derivative,
1622  const RCP<const Epetra_Map> &fnc_map,
1623  const RCP<const Epetra_Map> &var_map
1624  )
1625 {
1626  typedef ModelEvaluatorBase MEB;
1627  if(derivative.getLinearOp().get()) {
1628  return EpetraExt::ModelEvaluator::Derivative(
1629  Teuchos::rcp_const_cast<Epetra_Operator>(
1630  Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
1631  )
1632  );
1633  }
1634  else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
1635  return EpetraExt::ModelEvaluator::Derivative(
1636  EpetraExt::ModelEvaluator::DerivativeMultiVector(
1638  ( derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
1639  ? *fnc_map
1640  : *var_map
1641  )
1642  ,derivative.getDerivativeMultiVector().getMultiVector()
1643  )
1644  ,convert(derivative.getDerivativeMultiVector().getOrientation())
1645  )
1646  );
1647  }
1648  return EpetraExt::ModelEvaluator::Derivative();
1649 }
1650 EpetraExt::ModelEvaluator::MPDerivative
1651 Thyra::convert(
1652  const ModelEvaluatorBase::MPDerivative &derivative,
1653  const RCP<const Epetra_Map> &/* fnc_map */,
1654  const RCP<const Epetra_Map> &/* var_map */
1655  )
1656 {
1657  //typedef ModelEvaluatorBase MEB; // unused
1658  if(derivative.getLinearOp().get()) {
1659  return EpetraExt::ModelEvaluator::MPDerivative(
1660  derivative.getLinearOp()
1661  );
1662  }
1663  else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
1664  return EpetraExt::ModelEvaluator::MPDerivative(
1665  EpetraExt::ModelEvaluator::MPDerivativeMultiVector(
1666  derivative.getDerivativeMultiVector().getMultiVector()
1667  ,convert(derivative.getDerivativeMultiVector().getOrientation())
1668  )
1669  );
1670  }
1671  return EpetraExt::ModelEvaluator::MPDerivative();
1672 }
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible...
const ModelEvaluatorBase::InArgs< double > & getFinalPoint() const
ModelEvaluatorBase::InArgs< double > getLowerBounds() const
void setNominalValues(const ModelEvaluatorBase::InArgs< double > &nominalValues)
Set the nominal values.
bool DistributedGlobal() const
RCP< const EpetraExt::ModelEvaluator > getEpetraModel() const
void reportFinalPoint(const ModelEvaluatorBase::InArgs< double > &finalPoint, const bool wasSolved)
bool is_null(const boost::shared_ptr< T > &p)
basic_OSTab< char > OSTab
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< LinearOpBase< double > > create_W_op() const
ModelEvaluatorBase::InArgs< double > getNominalValues() const
RCP< Teuchos::ParameterList > getNonconstParameterList()
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< PreconditionerBase< double > > create_W_prec() const
Returns null currently.
void setStateVariableScalingVec(const RCP< const Epetra_Vector > &stateVariableScalingVec)
Set the state variable scaling vector s_x (see above).
ModelEvaluatorBase::InArgs< double > createInArgs() const
RCP< const Epetra_Vector > getStateVariableScalingVec() const
Get the inverse state variable scaling vector inv_s_x (see above).
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
RCP< const Teuchos::ParameterList > getParameterList() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RCP< const LinearOpWithSolveFactoryBase< double > > get_W_factory() const
RCP< const Epetra_Vector > getStateVariableInvScalingVec() const
Get the state variable scaling vector s_x (see above).
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
void setArgs(const InArgs< Scalar > &inArgs, bool ignoreUnsupported=false, bool cloneObjects=false)
Set non-null arguments (does not overwrite non-NULLs with NULLs) .
ModelEvaluatorBase::InArgs< double > getUpperBounds() const
void setStateFunctionScalingVec(const RCP< const Epetra_Vector > &stateFunctionScalingVec)
Set the state function scaling vector s_f (see above).
RCP< const Epetra_Vector > getStateFunctionScalingVec() const
Get the state function scaling vector s_f (see above).
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
void initialize(const RCP< const EpetraExt::ModelEvaluator > &epetraModel, const RCP< LinearOpWithSolveFactoryBase< double > > &W_factory)
RCP< const VectorSpaceBase< double > > get_f_space() const
TypeTo implicit_cast(const TypeFrom &t)
RCP< const VectorSpaceBase< double > > get_x_space() const
RCP< const VectorSpaceBase< double > > get_g_space(int j) const
Simple public strict containing properties of a derivative object.
RCP< const VectorSpaceBase< double > > get_p_space(int l) const
void resize(size_type new_size, const value_type &x=value_type())
RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)
const RCP< T > & assert_not_null() const
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Base subclass for ModelEvaluator that defines some basic types.
ModelEvaluatorBase::EDerivativeMultiVectorOrientation convert(const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation)
void uninitialize(RCP< const EpetraExt::ModelEvaluator > *epetraModel=NULL, RCP< LinearOpWithSolveFactoryBase< double > > *W_factory=NULL)
bool nonnull(const boost::shared_ptr< T > &p)
Teuchos::ArrayView< const std::string > get_g_names(int j) const
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
size_type size() const
KOKKOS_INLINE_FUNCTION constexpr unsigned rank(const View< D, P...> &V)
Determines the forms of a general derivative that are supported.
#define TEUCHOS_ASSERT(assertion_test)
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
RCP< Teuchos::ParameterList > unsetParameterList()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string typeName(const T &t)