Thyra Version of the Day
Loading...
Searching...
No Matches
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
57namespace {
58
59
60const std::string StateFunctionScaling_name = "State Function Scaling";
63 Thyra::EpetraModelEvaluator::EStateFunctionScaling
64 >
65 >
66stateFunctionScalingValidator;
67const std::string StateFunctionScaling_default = "None";
68
69
70// Extract out the Epetra_RowMatrix from the set W in an Epetra OutArgs object
72get_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
96create_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
117namespace 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,
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,
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;
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
535EpetraModelEvaluator::create_DfDp_op_impl(int /* l */) const
536{
538 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
539}
540
541
542RCP<LinearOpBase<double> >
543EpetraModelEvaluator::create_DgDx_dot_op_impl(int /* j */) const
544{
546 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
547}
548
549
550RCP<LinearOpBase<double> >
551EpetraModelEvaluator::create_DgDx_op_impl(int /* j */) const
552{
554 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
555}
556
557
558RCP<LinearOpBase<double> >
559EpetraModelEvaluator::create_DgDp_op_impl( int /* j */, int /* l */ ) const
560{
562 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
563}
564
565
566ModelEvaluatorBase::OutArgs<double>
567EpetraModelEvaluator::createOutArgsImpl() const
568{
569 if (!currentInArgsOutArgs_) {
570 updateInArgsOutArgs();
571 }
572 return prototypeOutArgs_;
573}
574
575
576void 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
804void EpetraModelEvaluator::convertInArgsFromEpetraToThyra(
805 const EpetraExt::ModelEvaluator::InArgs &epetraInArgs,
806 ModelEvaluatorBase::InArgs<double> *inArgs
807 ) const
808{
809
811 typedef ModelEvaluatorBase MEB;
812
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
847void 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
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
949void 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
1144void 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
1207void 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
1274void 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
1303void EpetraModelEvaluator::updateNominalValuesAndBounds() const
1304{
1305
1306 using Teuchos::rcp;
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
1376void 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
1503RCP<EpetraLinearOp>
1504EpetraModelEvaluator::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
1522Thyra::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
1532Thyra::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
1548EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
1549Thyra::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
1566Thyra::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
1604Thyra::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
1619EpetraExt::ModelEvaluator::Derivative
1620Thyra::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}
1650EpetraExt::ModelEvaluator::MPDerivative
1651Thyra::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}
size_type size() const
void resize(size_type new_size, const value_type &x=value_type())
const RCP< T > & assert_not_null() const
T * get() const
RCP< Teuchos::ParameterList > unsetParameterList()
RCP< Teuchos::ParameterList > getNonconstParameterList()
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).
void initialize(const RCP< const EpetraExt::ModelEvaluator > &epetraModel, const RCP< LinearOpWithSolveFactoryBase< double > > &W_factory)
RCP< const VectorSpaceBase< double > > get_f_space() const
void setStateVariableScalingVec(const RCP< const Epetra_Vector > &stateVariableScalingVec)
Set the state variable scaling vector s_x (see above).
RCP< LinearOpBase< double > > create_W_op() const
ModelEvaluatorBase::InArgs< double > getUpperBounds() const
void setNominalValues(const ModelEvaluatorBase::InArgs< double > &nominalValues)
Set the nominal values.
RCP< PreconditionerBase< double > > create_W_prec() const
Returns null currently.
RCP< const LinearOpWithSolveFactoryBase< double > > get_W_factory() const
ModelEvaluatorBase::InArgs< double > getNominalValues() const
RCP< const Teuchos::ParameterList > getParameterList() const
void reportFinalPoint(const ModelEvaluatorBase::InArgs< double > &finalPoint, const bool wasSolved)
ModelEvaluatorBase::InArgs< double > createInArgs() const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
RCP< const VectorSpaceBase< double > > get_g_space(int j) const
RCP< const EpetraExt::ModelEvaluator > getEpetraModel() const
void uninitialize(RCP< const EpetraExt::ModelEvaluator > *epetraModel=NULL, RCP< LinearOpWithSolveFactoryBase< double > > *W_factory=NULL)
RCP< const VectorSpaceBase< double > > get_x_space() const
RCP< const Epetra_Vector > getStateVariableInvScalingVec() const
Get the state variable scaling vector s_x (see above).
ModelEvaluatorBase::InArgs< double > getLowerBounds() const
RCP< const Epetra_Vector > getStateVariableScalingVec() const
Get the inverse state variable scaling vector inv_s_x (see above).
RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::ArrayView< const std::string > get_g_names(int j) const
RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
RCP< const VectorSpaceBase< double > > get_p_space(int l) const
const ModelEvaluatorBase::InArgs< double > & getFinalPoint() const
ModelEvaluatorBase::EDerivativeMultiVectorOrientation convert(const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation)
Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects.
Determines the forms of a general derivative that are supported.
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
void setArgs(const InArgs< Scalar > &inArgs, bool ignoreUnsupported=false, bool cloneObjects=false)
Set non-null arguments (does not overwrite non-NULLs with NULLs) .
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
Base subclass for ModelEvaluator that defines some basic types.
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.
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
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.
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.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const std::shared_ptr< T > &p)
bool nonnull(const std::shared_ptr< T > &p)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
TypeTo implicit_cast(const TypeFrom &t)
std::string typeName(const T &t)
basic_OSTab< char > OSTab
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)
Simple public strict containing properties of a derivative object.