Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_ModelEvaluatorDefaultBase.hpp
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#ifndef THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
43#define THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
44
45#include "Thyra_VectorBase.hpp"
46
47#include "Thyra_ModelEvaluator.hpp"
48#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
49
50
51#ifdef HAVE_THYRA_ME_POLYNOMIAL
52
53
54// Define the polynomial traits class specializtaion
55// Teuchos::PolynomialTraits<VectorBase > before there is any chance of an
56// instantiation requiring the definition of this class. The Intel C++ 9.1
57// compiler is having trouble compiling Thyra_EpetraModelEvaluator.cpp because
58// Thyra_ModelEvaluatorBase_decl.hpp is only including
59// "Teuchos_Polynomial.hpp" which does not contain the needed specialization.
60// By including it here, all client code is likely to compile just fine. I am
61// going through all of the in order to avoid having to put
62// Thyra_PolynomialVectorTraits.hpp in the interfaces directory. The problem
63// with putting Thyra_PolynomialVectorTraits.hpp the interfaces directory is
64// that it contains calls to code in the support directory and creates an
65// unacceptable circular dependency that will cause problems once we move to
66// subpackages in the CMake build system.
67#include "Thyra_PolynomialVectorTraits.hpp"
68
69
70#endif // HAVE_THYRA_ME_POLYNOMIAL
71
72
73namespace Thyra {
74
75
76//
77// Undocumentated (from the user's perspective) types that are used to
78// implement ModelEvaluatorDefaultBase. These types are defined outside of
79// the templated class since they do nt depend on the scalar type.
80//
81
82
83namespace ModelEvaluatorDefaultBaseTypes {
84
85
86// Type used to determine if the ModelEvaluatorDefaultBase implementation will
87// provide a LinearOpBase wrapper for derivative object given in
88// MultiVectorBase form.
89class DefaultDerivLinearOpSupport {
90public:
91 DefaultDerivLinearOpSupport()
92 :provideDefaultLinearOp_(false),
93 mvImplOrientation_(ModelEvaluatorBase::DERIV_MV_BY_COL)
94 {}
95 DefaultDerivLinearOpSupport(
97 )
98 :provideDefaultLinearOp_(true),
99 mvImplOrientation_(mvImplOrientation_in)
100 {}
101 bool provideDefaultLinearOp() const
102 { return provideDefaultLinearOp_; }
104 { return mvImplOrientation_; }
105private:
106 bool provideDefaultLinearOp_;
108};
109
110
111// Type used to determine if the ModelEvaluatorDefaultBase implementation will
112// provide an explicit transpose copy of a derivative object given in
113// MultiVectorBase form.
114class DefaultDerivMvAdjointSupport {
115public:
116 DefaultDerivMvAdjointSupport()
117 :provideDefaultAdjoint_(false),
118 mvAdjointCopyOrientation_(ModelEvaluatorBase::DERIV_MV_BY_COL)
119 {}
120 DefaultDerivMvAdjointSupport(
121 const ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvAdjointCopyOrientation_in
122 )
123 :provideDefaultAdjoint_(true),
124 mvAdjointCopyOrientation_(mvAdjointCopyOrientation_in)
125 {}
126 bool provideDefaultAdjoint() const
127 { return provideDefaultAdjoint_; }
128 ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvAdjointCopyOrientation() const
129 { return mvAdjointCopyOrientation_; }
130private:
131 bool provideDefaultAdjoint_;
133};
134
135
136// Type used to remember a pair of transposed multi-vectors to implement a
137// adjoint copy.
138template<class Scalar>
139struct MultiVectorAdjointPair {
140 MultiVectorAdjointPair()
141 {}
142 MultiVectorAdjointPair(
143 const RCP<MultiVectorBase<Scalar> > &in_mvOuter,
144 const RCP<const MultiVectorBase<Scalar> > &in_mvImplAdjoint
145 )
146 : mvOuter(in_mvOuter),
147 mvImplAdjoint(in_mvImplAdjoint)
148 {}
149 RCP<MultiVectorBase<Scalar> > mvOuter;
150 RCP<const MultiVectorBase<Scalar> > mvImplAdjoint;
151private:
152};
153
154
155
156
157} // namespace ModelEvaluatorDefaultBaseTypes
158
159
187template<class Scalar>
188class ModelEvaluatorDefaultBase : virtual public ModelEvaluator<Scalar>
189{
190public:
191
194
196 int Np() const;
198 int Ng() const;
215 ) const;
220
221#ifdef Thyra_BUILD_HESSIAN_SUPPORT
223 virtual RCP<LinearOpBase<Scalar> > create_hess_f_xx() const;
225 virtual RCP<LinearOpBase<Scalar> > create_hess_f_xp(int l) const;
227 virtual RCP<LinearOpBase<Scalar> > create_hess_f_pp( int l1, int l2 ) const;
229 virtual RCP<LinearOpBase<Scalar> > create_hess_g_xx(int j) const;
231 virtual RCP<LinearOpBase<Scalar> > create_hess_g_xp( int j, int l ) const;
233 virtual RCP<LinearOpBase<Scalar> > create_hess_g_pp( int j, int l1, int l2 ) const;
234#endif // ifdef Thyra_BUILD_HESSIAN_SUPPORT
235
237
238protected:
239
242
252
258
260
261private:
262
265
267 virtual RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
268
270 virtual RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
271
273 virtual RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
274
276 virtual RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
277
279
282
284 virtual ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const = 0;
285
287 virtual void evalModelImpl(
290 ) const = 0;
291
293
294protected:
295
298
299private:
300
301 // //////////////////////////////
302 // Private tpyes
303
304 typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
305 DefaultDerivLinearOpSupport;
306
307 typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
308 DefaultDerivMvAdjointSupport;
309
310 typedef ModelEvaluatorDefaultBaseTypes::MultiVectorAdjointPair<Scalar>
311 MultiVectorAdjointPair;
312
313 // //////////////////////////////
314 // Private data members
315
316 bool isInitialized_;
317
318 Array<DefaultDerivLinearOpSupport> DfDp_default_op_support_;
319
320 Array<DefaultDerivLinearOpSupport> DgDx_dot_default_op_support_;
321
322 Array<DefaultDerivLinearOpSupport> DgDx_default_op_support_;
323
324 Array<Array<DefaultDerivLinearOpSupport> > DgDp_default_op_support_;
325 Array<Array<DefaultDerivMvAdjointSupport> > DgDp_default_mv_support_;
326
327 bool default_W_support_;
328
329 ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
330
331 // //////////////////////////////
332 // Private member functions
333
334 void lazyInitializeDefaultBase() const;
335
336 void assert_l(const int l) const;
337
338 void assert_j(const int j) const;
339
340 // //////////////////////////////
341 // Private static functions
342
343 static DefaultDerivLinearOpSupport
344 determineDefaultDerivLinearOpSupport(
345 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
346 );
347
349 createDefaultLinearOp(
350 const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
351 const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
352 const RCP<const VectorSpaceBase<Scalar> > &var_space
353 );
354
356 updateDefaultLinearOpSupport(
357 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
358 const DefaultDerivLinearOpSupport &defaultLinearOpSupport
359 );
360
362 getOutArgImplForDefaultLinearOpSupport(
364 const DefaultDerivLinearOpSupport &defaultLinearOpSupport
365 );
366
367 static DefaultDerivMvAdjointSupport
368 determineDefaultDerivMvAdjointSupport(
369 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
370 const VectorSpaceBase<Scalar> &fnc_space,
371 const VectorSpaceBase<Scalar> &var_space
372 );
373
375 updateDefaultDerivMvAdjointSupport(
376 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
377 const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
378 );
379
380};
381
382
383} // namespace Thyra
384
385
386//
387// Implementations
388//
389
390
391#include "Thyra_ModelEvaluatorHelpers.hpp"
392#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
393#include "Thyra_DetachedMultiVectorView.hpp"
394#include "Teuchos_Assert.hpp"
395
396
397namespace Thyra {
398
399
400// Overridden from ModelEvaluator
401
402
403template<class Scalar>
405{
406 lazyInitializeDefaultBase();
407 return prototypeOutArgs_.Np();
408}
409
410
411template<class Scalar>
413{
414 lazyInitializeDefaultBase();
415 return prototypeOutArgs_.Ng();
416}
417
418
419template<class Scalar>
422{
423 lazyInitializeDefaultBase();
424 if (default_W_support_)
425 return this->get_W_factory()->createOp();
426 return Teuchos::null;
427}
428
429
430template<class Scalar>
433{
434 lazyInitializeDefaultBase();
435#ifdef TEUCHOS_DEBUG
436 assert_l(l);
437#endif
438 const DefaultDerivLinearOpSupport
439 defaultLinearOpSupport = DfDp_default_op_support_[l];
440 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
441 return createDefaultLinearOp(
442 defaultLinearOpSupport,
443 this->get_f_space(),
444 this->get_p_space(l)
445 );
446 }
447 return this->create_DfDp_op_impl(l);
448}
449
450
451template<class Scalar>
454{
455 lazyInitializeDefaultBase();
456#ifdef TEUCHOS_DEBUG
457 assert_j(j);
458#endif
459 const DefaultDerivLinearOpSupport
460 defaultLinearOpSupport = DgDx_dot_default_op_support_[j];
461 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
462 return createDefaultLinearOp(
463 defaultLinearOpSupport,
464 this->get_g_space(j),
465 this->get_x_space()
466 );
467 }
468 return this->create_DgDx_dot_op_impl(j);
469}
470
471
472template<class Scalar>
475{
476 lazyInitializeDefaultBase();
477#ifdef TEUCHOS_DEBUG
478 assert_j(j);
479#endif
480 const DefaultDerivLinearOpSupport
481 defaultLinearOpSupport = DgDx_default_op_support_[j];
482 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
483 return createDefaultLinearOp(
484 defaultLinearOpSupport,
485 this->get_g_space(j),
486 this->get_x_space()
487 );
488 }
489 return this->create_DgDx_op_impl(j);
490}
491
492
493template<class Scalar>
496{
497 lazyInitializeDefaultBase();
498#ifdef TEUCHOS_DEBUG
499 assert_j(j);
500 assert_l(l);
501#endif
502 const DefaultDerivLinearOpSupport
503 defaultLinearOpSupport = DgDp_default_op_support_[j][l];
504 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
505 return createDefaultLinearOp(
506 defaultLinearOpSupport,
507 this->get_g_space(j),
508 this->get_p_space(l)
509 );
510 }
511 return this->create_DgDp_op_impl(j,l);
512}
513
514
515template<class Scalar>
518{
519 lazyInitializeDefaultBase();
520 return prototypeOutArgs_;
521}
522
523
524template<class Scalar>
528 ) const
529{
530
531 using Teuchos::outArg;
532 typedef ModelEvaluatorBase MEB;
533
534 lazyInitializeDefaultBase();
535
536 const int l_Np = outArgs.Np();
537 const int l_Ng = outArgs.Ng();
538
539 //
540 // A) Assert that the inArgs and outArgs object match this class!
541 //
542
543#ifdef TEUCHOS_DEBUG
544 assertInArgsEvalObjects(*this,inArgs);
545 assertOutArgsEvalObjects(*this,outArgs,&inArgs);
546#endif
547
548 //
549 // B) Setup the OutArgs object for the underlying implementation's
550 // evalModelImpl(...) function
551 //
552
553 MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
554 Array<MultiVectorAdjointPair> DgDp_temp_adjoint_copies;
555
556 {
557
558 outArgsImpl.setArgs(outArgs,true);
559
560 // DfDp(l)
561 if (outArgsImpl.supports(MEB::OUT_ARG_f)) {
562 for ( int l = 0; l < l_Np; ++l ) {
563 const DefaultDerivLinearOpSupport defaultLinearOpSupport =
564 DfDp_default_op_support_[l];
565 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
566 outArgsImpl.set_DfDp( l,
567 getOutArgImplForDefaultLinearOpSupport(
568 outArgs.get_DfDp(l), defaultLinearOpSupport
569 )
570 );
571 }
572 else {
573 // DfDp(l) already set by outArgsImpl.setArgs(...)!
574 }
575 }
576 }
577
578 // DgDx_dot(j)
579 for ( int j = 0; j < l_Ng; ++j ) {
580 const DefaultDerivLinearOpSupport defaultLinearOpSupport =
581 DgDx_dot_default_op_support_[j];
582 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
583 outArgsImpl.set_DgDx_dot( j,
584 getOutArgImplForDefaultLinearOpSupport(
585 outArgs.get_DgDx_dot(j), defaultLinearOpSupport
586 )
587 );
588 }
589 else {
590 // DgDx_dot(j) already set by outArgsImpl.setArgs(...)!
591 }
592 }
593
594 // DgDx(j)
595 for ( int j = 0; j < l_Ng; ++j ) {
596 const DefaultDerivLinearOpSupport defaultLinearOpSupport =
597 DgDx_default_op_support_[j];
598 if (defaultLinearOpSupport.provideDefaultLinearOp()) {
599 outArgsImpl.set_DgDx( j,
600 getOutArgImplForDefaultLinearOpSupport(
601 outArgs.get_DgDx(j), defaultLinearOpSupport
602 )
603 );
604 }
605 else {
606 // DgDx(j) already set by outArgsImpl.setArgs(...)!
607 }
608 }
609
610 // DgDp(j,l)
611 for ( int j = 0; j < l_Ng; ++j ) {
612 const Array<DefaultDerivLinearOpSupport> &DgDp_default_op_support_j =
613 DgDp_default_op_support_[j];
614 const Array<DefaultDerivMvAdjointSupport> &DgDp_default_mv_support_j =
615 DgDp_default_mv_support_[j];
616 for ( int l = 0; l < l_Np; ++l ) {
617 const DefaultDerivLinearOpSupport defaultLinearOpSupport =
618 DgDp_default_op_support_j[l];
619 const DefaultDerivMvAdjointSupport defaultMvAdjointSupport =
620 DgDp_default_mv_support_j[l];
621 MEB::Derivative<Scalar> DgDp_j_l;
622 if (!outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none())
623 DgDp_j_l = outArgs.get_DgDp(j,l);
624 if (
625 defaultLinearOpSupport.provideDefaultLinearOp()
626 && !is_null(DgDp_j_l.getLinearOp())
627 )
628 {
629 outArgsImpl.set_DgDp( j, l,
630 getOutArgImplForDefaultLinearOpSupport(
631 DgDp_j_l, defaultLinearOpSupport
632 )
633 );
634 }
635 else if (
636 defaultMvAdjointSupport.provideDefaultAdjoint()
637 && !is_null(DgDp_j_l.getMultiVector())
638 )
639 {
640 const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv =
641 DgDp_j_l.getMultiVector();
642 if (
643 defaultMvAdjointSupport.mvAdjointCopyOrientation()
644 ==
645 DgDp_j_l.getMultiVectorOrientation()
646 )
647 {
648 // The orientation of the multi-vector is different so we need to
649 // create a temporary copy to pass to evalModelImpl(...) and then
650 // copy it back again!
651 const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv_adj =
652 createMembers(DgDp_j_l_mv->domain(), DgDp_j_l_mv->range()->dim());
653 outArgsImpl.set_DgDp( j, l,
654 MEB::Derivative<Scalar>(
655 DgDp_j_l_mv_adj,
656 getOtherDerivativeMultiVectorOrientation(
657 defaultMvAdjointSupport.mvAdjointCopyOrientation()
658 )
659 )
660 );
661 // Remember these multi-vectors so that we can do the transpose copy
662 // back after the evaluation!
663 DgDp_temp_adjoint_copies.push_back(
664 MultiVectorAdjointPair(DgDp_j_l_mv, DgDp_j_l_mv_adj)
665 );
666 }
667 else {
668 // The form of the multi-vector is supported by evalModelImpl(..)
669 // and is already set on the outArgsImpl object.
670 }
671 }
672 else {
673 // DgDp(j,l) already set by outArgsImpl.setArgs(...)!
674 }
675 }
676 }
677
678 // W
679 {
681 if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
683 W_factory = this->get_W_factory();
684 // Extract the underlying W_op object (if it exists)
686 uninitializeOp<Scalar>(*W_factory, W.ptr(), outArg(W_op_const));
688 if (!is_null(W_op_const)) {
689 // Here we remove the const. This is perfectly safe since
690 // w.r.t. this class, we put a non-const object in there and we can
691 // expect to change that object after the fact. That is our
692 // prerogative.
693 W_op = Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(W_op_const);
694 }
695 else {
696 // The W_op object has not been initialized yet so create it. The
697 // next time through, we should not have to do this!
698 W_op = this->create_W_op();
699 }
700 outArgsImpl.set_W_op(W_op);
701 }
702 }
703 }
704
705 //
706 // C) Evaluate the underlying model implementation!
707 //
708
709 this->evalModelImpl( inArgs, outArgsImpl );
710
711 //
712 // D) Post-process the output arguments
713 //
714
715 // Do explicit transposes for DgDp(j,l) if needed
716 const int numMvAdjointCopies = DgDp_temp_adjoint_copies.size();
717 for ( int adj_copy_i = 0; adj_copy_i < numMvAdjointCopies; ++adj_copy_i ) {
718 const MultiVectorAdjointPair adjPair =
719 DgDp_temp_adjoint_copies[adj_copy_i];
720 doExplicitMultiVectorAdjoint( *adjPair.mvImplAdjoint, &*adjPair.mvOuter );
721 }
722
723 // Update W given W_op and W_factory
724 {
726 if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
728 W_factory = this->get_W_factory();
729 W_factory->setOStream(this->getOStream());
730 W_factory->setVerbLevel(this->getVerbLevel());
731 initializeOp<Scalar>(*W_factory, outArgsImpl.get_W_op().getConst(), W.ptr());
732 }
733 }
734
735}
736
737
738// protected
739
740
741// Setup functions called by subclasses
742
743template<class Scalar>
745{
746
747 typedef ModelEvaluatorBase MEB;
748
749 // In case we throw half way thorugh, set to uninitialized
750 isInitialized_ = false;
751 default_W_support_ = false;
752
753 //
754 // A) Get the InArgs and OutArgs from the subclass
755 //
756
757 const MEB::InArgs<Scalar> inArgs = this->createInArgs();
758 const MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
759
760 //
761 // B) Validate the subclasses InArgs and OutArgs
762 //
763
764#ifdef TEUCHOS_DEBUG
765 assertInArgsOutArgsSetup( this->description(), inArgs, outArgsImpl );
766#endif // TEUCHOS_DEBUG
767
768 //
769 // C) Set up support for default derivative objects and prototype OutArgs
770 //
771
772 const int l_Ng = outArgsImpl.Ng();
773 const int l_Np = outArgsImpl.Np();
774
775 // Set support for all outputs supported in the underly implementation
776 MEB::OutArgsSetup<Scalar> outArgs;
777 outArgs.setModelEvalDescription(this->description());
778 outArgs.set_Np_Ng(l_Np,l_Ng);
779 outArgs.setSupports(outArgsImpl);
780
781 // DfDp
782 DfDp_default_op_support_.clear();
783 if (outArgs.supports(MEB::OUT_ARG_f)) {
784 for ( int l = 0; l < l_Np; ++l ) {
785 const MEB::DerivativeSupport DfDp_l_impl_support =
786 outArgsImpl.supports(MEB::OUT_ARG_DfDp,l);
787 const DefaultDerivLinearOpSupport DfDp_l_op_support =
788 determineDefaultDerivLinearOpSupport(DfDp_l_impl_support);
789 DfDp_default_op_support_.push_back(DfDp_l_op_support);
790 outArgs.setSupports(
791 MEB::OUT_ARG_DfDp, l,
792 updateDefaultLinearOpSupport(
793 DfDp_l_impl_support, DfDp_l_op_support
794 )
795 );
796 }
797 }
798
799 // DgDx_dot
800 DgDx_dot_default_op_support_.clear();
801 for ( int j = 0; j < l_Ng; ++j ) {
802 const MEB::DerivativeSupport DgDx_dot_j_impl_support =
803 outArgsImpl.supports(MEB::OUT_ARG_DgDx_dot,j);
804 const DefaultDerivLinearOpSupport DgDx_dot_j_op_support =
805 determineDefaultDerivLinearOpSupport(DgDx_dot_j_impl_support);
806 DgDx_dot_default_op_support_.push_back(DgDx_dot_j_op_support);
807 outArgs.setSupports(
808 MEB::OUT_ARG_DgDx_dot, j,
809 updateDefaultLinearOpSupport(
810 DgDx_dot_j_impl_support, DgDx_dot_j_op_support
811 )
812 );
813 }
814
815 // DgDx
816 DgDx_default_op_support_.clear();
817 for ( int j = 0; j < l_Ng; ++j ) {
818 const MEB::DerivativeSupport DgDx_j_impl_support =
819 outArgsImpl.supports(MEB::OUT_ARG_DgDx,j);
820 const DefaultDerivLinearOpSupport DgDx_j_op_support =
821 determineDefaultDerivLinearOpSupport(DgDx_j_impl_support);
822 DgDx_default_op_support_.push_back(DgDx_j_op_support);
823 outArgs.setSupports(
824 MEB::OUT_ARG_DgDx, j,
825 updateDefaultLinearOpSupport(
826 DgDx_j_impl_support, DgDx_j_op_support
827 )
828 );
829 }
830
831 // DgDp
832 DgDp_default_op_support_.clear();
833 DgDp_default_mv_support_.clear();
834 for ( int j = 0; j < l_Ng; ++j ) {
835 DgDp_default_op_support_.push_back(Array<DefaultDerivLinearOpSupport>());
836 DgDp_default_mv_support_.push_back(Array<DefaultDerivMvAdjointSupport>());
837 for ( int l = 0; l < l_Np; ++l ) {
838 const MEB::DerivativeSupport DgDp_j_l_impl_support =
839 outArgsImpl.supports(MEB::OUT_ARG_DgDp,j,l);
840 // LinearOpBase support
841 const DefaultDerivLinearOpSupport DgDp_j_l_op_support =
842 determineDefaultDerivLinearOpSupport(DgDp_j_l_impl_support);
843 DgDp_default_op_support_[j].push_back(DgDp_j_l_op_support);
844 outArgs.setSupports(
845 MEB::OUT_ARG_DgDp, j, l,
846 updateDefaultLinearOpSupport(
847 DgDp_j_l_impl_support, DgDp_j_l_op_support
848 )
849 );
850 // MultiVectorBase
851 const DefaultDerivMvAdjointSupport DgDp_j_l_mv_support =
852 determineDefaultDerivMvAdjointSupport(
853 DgDp_j_l_impl_support, *this->get_g_space(j), *this->get_p_space(l)
854 );
855 DgDp_default_mv_support_[j].push_back(DgDp_j_l_mv_support);
856 outArgs.setSupports(
857 MEB::OUT_ARG_DgDp, j, l,
858 updateDefaultDerivMvAdjointSupport(
859 outArgs.supports(MEB::OUT_ARG_DgDp, j, l),
860 DgDp_j_l_mv_support
861 )
862 );
863 }
864 }
865 // 2007/09/09: rabart: ToDo: Move the above code into a private helper
866 // function!
867
868 // W (given W_op and W_factory)
869 default_W_support_ = false;
870 if ( outArgsImpl.supports(MEB::OUT_ARG_W_op) && !is_null(this->get_W_factory())
871 && !outArgsImpl.supports(MEB::OUT_ARG_W) )
872 {
873 default_W_support_ = true;
874 outArgs.setSupports(MEB::OUT_ARG_W);
875 outArgs.set_W_properties(outArgsImpl.get_W_properties());
876 }
877
878 //
879 // D) All done!
880 //
881
882 prototypeOutArgs_ = outArgs;
883 isInitialized_ = true;
884
885}
886
887template<class Scalar>
889{
890 isInitialized_ = false;
891}
892
893// Private functions with default implementaton to be overridden by subclasses
894
895
896template<class Scalar>
899{
900 typedef ModelEvaluatorBase MEB;
901 MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
903 outArgs.supports(MEB::OUT_ARG_DfDp,l).supports(MEB::DERIV_LINEAR_OP),
904 std::logic_error,
905 "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
906 " supports the LinearOpBase form of DfDp("<<l<<") (as determined from its"
907 " OutArgs object created by createOutArgsImpl())"
908 " but this function create_DfDp_op_impl(...) has not been overridden"
909 " to create such an object!"
910 );
911 return Teuchos::null;
912}
913
914
915template<class Scalar>
916RCP<LinearOpBase<Scalar> >
917ModelEvaluatorDefaultBase<Scalar>::create_DgDx_dot_op_impl(int j) const
918{
919 typedef ModelEvaluatorBase MEB;
920 MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
922 outArgs.supports(MEB::OUT_ARG_DgDx_dot,j).supports(MEB::DERIV_LINEAR_OP),
923 std::logic_error,
924 "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
925 " supports the LinearOpBase form of DgDx_dot("<<j<<") (as determined from"
926 " its OutArgs object created by createOutArgsImpl())"
927 " but this function create_DgDx_dot_op_impl(...) has not been overridden"
928 " to create such an object!"
929 );
930 return Teuchos::null;
931}
932
933
934template<class Scalar>
935RCP<LinearOpBase<Scalar> >
936ModelEvaluatorDefaultBase<Scalar>::create_DgDx_op_impl(int j) const
937{
938 typedef ModelEvaluatorBase MEB;
939 MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
941 outArgs.supports(MEB::OUT_ARG_DgDx,j).supports(MEB::DERIV_LINEAR_OP),
942 std::logic_error,
943 "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
944 " supports the LinearOpBase form of DgDx("<<j<<") (as determined from"
945 " its OutArgs object created by createOutArgsImpl())"
946 " but this function create_DgDx_op_impl(...) has not been overridden"
947 " to create such an object!"
948 );
949 return Teuchos::null;
950}
951
952
953template<class Scalar>
954RCP<LinearOpBase<Scalar> >
955ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op_impl(int j, int l) const
956{
957 typedef ModelEvaluatorBase MEB;
958 MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
960 outArgs.supports(MEB::OUT_ARG_DgDp,j,l).supports(MEB::DERIV_LINEAR_OP),
961 std::logic_error,
962 "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
963 " supports the LinearOpBase form of DgDp("<<j<<","<<l<<")"
964 " (as determined from its OutArgs object created by createOutArgsImpl())"
965 " but this function create_DgDp_op_impl(...) has not been overridden"
966 " to create such an object!"
967 );
968 return Teuchos::null;
969}
970
971
972template<class Scalar>
973RCP<const VectorSpaceBase<Scalar> >
975{
976 return this->get_f_space();
977}
978
979template<class Scalar>
982{
983 return this->get_g_space(j);
984}
985
986#ifdef Thyra_BUILD_HESSIAN_SUPPORT
987
988template<class Scalar>
991{
992 return Teuchos::null;
993}
994
995template<class Scalar>
996RCP<LinearOpBase<Scalar> >
997ModelEvaluatorDefaultBase<Scalar>::create_hess_f_xp(int l) const
998{
999 return Teuchos::null;
1000}
1001
1002template<class Scalar>
1003RCP<LinearOpBase<Scalar> >
1004ModelEvaluatorDefaultBase<Scalar>::create_hess_f_pp( int l1, int l2 ) const
1005{
1006 return Teuchos::null;
1007}
1008
1009template<class Scalar>
1010RCP<LinearOpBase<Scalar> >
1011ModelEvaluatorDefaultBase<Scalar>::create_hess_g_xx(int j) const
1012{
1013 return Teuchos::null;
1014}
1015
1016template<class Scalar>
1017RCP<LinearOpBase<Scalar> >
1018ModelEvaluatorDefaultBase<Scalar>::create_hess_g_xp( int j, int l ) const
1019{
1020 return Teuchos::null;
1021}
1022
1023template<class Scalar>
1024RCP<LinearOpBase<Scalar> >
1025ModelEvaluatorDefaultBase<Scalar>::create_hess_g_pp( int j, int l1, int l2 ) const
1026{
1027 return Teuchos::null;
1028}
1029
1030#endif // ifdef Thyra_BUILD_HESSIAN_SUPPORT
1031
1032// protected
1033
1034
1035template<class Scalar>
1037 :isInitialized_(false), default_W_support_(false)
1038{}
1039
1040
1041// private
1042
1043
1044template<class Scalar>
1046{
1047 if (!isInitialized_)
1048 const_cast<ModelEvaluatorDefaultBase<Scalar>*>(this)->initializeDefaultBase();
1049}
1050
1051
1052template<class Scalar>
1053void ModelEvaluatorDefaultBase<Scalar>::assert_l(const int l) const
1054{
1056}
1057
1058
1059template<class Scalar>
1060void ModelEvaluatorDefaultBase<Scalar>::assert_j(const int j) const
1061{
1063}
1064
1065
1066// Private static functions
1067
1068
1069template<class Scalar>
1070ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
1071ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivLinearOpSupport(
1072 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
1073 )
1074{
1075 typedef ModelEvaluatorBase MEB;
1076 if (
1077 (
1078 derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1079 ||
1080 derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW)
1081 )
1082 &&
1083 !derivSupportImpl.supports(MEB::DERIV_LINEAR_OP)
1084 )
1085 {
1086 return DefaultDerivLinearOpSupport(
1087 derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1088 ? MEB::DERIV_MV_BY_COL
1089 : MEB::DERIV_TRANS_MV_BY_ROW
1090 );
1091 }
1092 return DefaultDerivLinearOpSupport();
1093}
1094
1095
1096template<class Scalar>
1097RCP<LinearOpBase<Scalar> >
1098ModelEvaluatorDefaultBase<Scalar>::createDefaultLinearOp(
1099 const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
1100 const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
1101 const RCP<const VectorSpaceBase<Scalar> > &var_space
1102 )
1103{
1104 using Teuchos::rcp_implicit_cast;
1105 typedef LinearOpBase<Scalar> LOB;
1106 typedef ModelEvaluatorBase MEB;
1107 switch(defaultLinearOpSupport.mvImplOrientation()) {
1108 case MEB::DERIV_MV_BY_COL:
1109 // The MultiVector will do just fine as the LinearOpBase
1110 return createMembers(fnc_space, var_space->dim());
1111 case MEB::DERIV_TRANS_MV_BY_ROW:
1112 // We will have to implicitly transpose the underlying MultiVector
1113 return nonconstAdjoint<Scalar>(
1114 rcp_implicit_cast<LOB>(createMembers(var_space, fnc_space->dim()))
1115 );
1116#ifdef TEUCHOS_DEBUG
1117 default:
1119#endif
1120 }
1121 return Teuchos::null; // Will never be called!
1122}
1123
1124
1125template<class Scalar>
1126ModelEvaluatorBase::DerivativeSupport
1127ModelEvaluatorDefaultBase<Scalar>::updateDefaultLinearOpSupport(
1128 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1129 const DefaultDerivLinearOpSupport &defaultLinearOpSupport
1130 )
1131{
1132 typedef ModelEvaluatorBase MEB;
1133 MEB::DerivativeSupport derivSupport = derivSupportImpl;
1134 if (defaultLinearOpSupport.provideDefaultLinearOp())
1135 derivSupport.plus(MEB::DERIV_LINEAR_OP);
1136 return derivSupport;
1137}
1138
1139
1140template<class Scalar>
1141ModelEvaluatorBase::Derivative<Scalar>
1142ModelEvaluatorDefaultBase<Scalar>::getOutArgImplForDefaultLinearOpSupport(
1143 const ModelEvaluatorBase::Derivative<Scalar> &deriv,
1144 const DefaultDerivLinearOpSupport &defaultLinearOpSupport
1145 )
1146{
1147
1148 using Teuchos::rcp_dynamic_cast;
1149 typedef ModelEvaluatorBase MEB;
1150 typedef MultiVectorBase<Scalar> MVB;
1151 typedef ScaledAdjointLinearOpBase<Scalar> SALOB;
1152
1153 // If the derivative is not a LinearOpBase then it it either empty or is a
1154 // MultiVectorBase object, and in either case we just return it since the
1155 // underlying evalModelImpl(...) functions should be able to handle it!
1156 if (is_null(deriv.getLinearOp()))
1157 return deriv;
1158
1159 // The derivative is LinearOpBase so get out the underlying MultiVectorBase
1160 // object and return its derivative multi-vector form.
1161 switch(defaultLinearOpSupport.mvImplOrientation()) {
1162 case MEB::DERIV_MV_BY_COL: {
1163 return MEB::Derivative<Scalar>(
1164 rcp_dynamic_cast<MVB>(deriv.getLinearOp(),true),
1165 MEB::DERIV_MV_BY_COL
1166 );
1167 }
1168 case MEB::DERIV_TRANS_MV_BY_ROW: {
1169 return MEB::Derivative<Scalar>(
1170 rcp_dynamic_cast<MVB>(
1171 rcp_dynamic_cast<SALOB>(deriv.getLinearOp(),true)->getNonconstOrigOp()
1172 ),
1173 MEB::DERIV_TRANS_MV_BY_ROW
1174 );
1175 }
1176#ifdef TEUCHOS_DEBUG
1177 default:
1179#endif
1180 }
1181
1182 return ModelEvaluatorBase::Derivative<Scalar>(); // Should never get here!
1183
1184}
1185
1186
1187template<class Scalar>
1188ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
1189ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivMvAdjointSupport(
1190 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1191 const VectorSpaceBase<Scalar> &fnc_space,
1192 const VectorSpaceBase<Scalar> &var_space
1193 )
1194{
1195 typedef ModelEvaluatorBase MEB;
1196 // Here we will support the adjoint copy for of a multi-vector if both
1197 // spaces give rise to in-core vectors.
1198 const bool implSupportsMv =
1199 ( derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1200 || derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
1201 const bool implLacksMvOrientSupport =
1202 ( !derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1203 || !derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
1204 const bool bothSpacesHaveInCoreViews =
1205 ( fnc_space.hasInCoreView() && var_space.hasInCoreView() );
1206 if ( implSupportsMv && implLacksMvOrientSupport && bothSpacesHaveInCoreViews ) {
1207 return DefaultDerivMvAdjointSupport(
1208 derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1209 ? MEB::DERIV_TRANS_MV_BY_ROW
1210 : MEB::DERIV_MV_BY_COL
1211 );
1212 }
1213 // We can't provide an adjoint copy or such a copy is not needed!
1214 return DefaultDerivMvAdjointSupport();
1215}
1216
1217
1218template<class Scalar>
1219ModelEvaluatorBase::DerivativeSupport
1220ModelEvaluatorDefaultBase<Scalar>::updateDefaultDerivMvAdjointSupport(
1221 const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1222 const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
1223 )
1224{
1225 typedef ModelEvaluatorBase MEB;
1226 MEB::DerivativeSupport derivSupport = derivSupportImpl;
1227 if (defaultMvAdjointSupport.provideDefaultAdjoint())
1228 derivSupport.plus(defaultMvAdjointSupport.mvAdjointCopyOrientation());
1229 return derivSupport;
1230}
1231
1232
1233} // namespace Thyra
1234
1235
1236#endif // THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
size_type size() const
void push_back(const value_type &x)
RCP< const T > getConst() const
Ptr< T > ptr() const
Determines the forms of a general derivative that are supported.
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Derivative< Scalar > get_DfDp(int l) const
Precondition: supports(OUT_ARG_DfDp,l)==true.
RCP< LinearOpWithSolveBase< Scalar > > get_W() const
Precondition: supports(OUT_ARG_W)==true.
int Ng() const
Return the number of axillary response functions g(j)(...) supported (Ng >= 0).
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
Derivative< Scalar > get_DgDx_dot(int j) const
Precondition: supports(OUT_ARG_DgDx_dot,j)==true.
Derivative< Scalar > get_DgDx(int j) const
Precondition: supports(OUT_ARG_DgDx,j)==true.
bool supports(EOutArgsMembers arg) const
Determine if an input argument is supported or not.
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
Base subclass for ModelEvaluator that defines some basic types.
Default base class for concrete model evaluators.
RCP< LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
RCP< LinearOpBase< Scalar > > create_DgDx_op(int j) const
virtual RCP< const VectorSpaceBase< Scalar > > get_g_multiplier_space(int j) const
RCP< LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
RCP< LinearOpWithSolveBase< Scalar > > create_W() const
virtual RCP< const VectorSpaceBase< Scalar > > get_f_multiplier_space() const
void evalModel(const ModelEvaluatorBase::InArgs< Scalar > &inArgs, const ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
ModelEvaluatorBase::OutArgs< Scalar > createOutArgs() const
void resetDefaultBase()
Sets the the DefaultBase to an uninitialized state, forcing lazy initialization when needed.
void initializeDefaultBase()
Function called by subclasses to fully initialize this object on any important change.
RCP< LinearOpBase< Scalar > > create_DfDp_op(int l) const
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Abstract interface for objects that represent a space for vectors.
#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)