Thyra  Version of the Day
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 
73 namespace 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 
83 namespace 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.
89 class DefaultDerivLinearOpSupport {
90 public:
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_; }
105 private:
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.
114 class DefaultDerivMvAdjointSupport {
115 public:
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_; }
130 private:
131  bool provideDefaultAdjoint_;
133 };
134 
135 
136 // Type used to remember a pair of transposed multi-vectors to implement a
137 // adjoint copy.
138 template<class Scalar>
139 struct 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;
151 private:
152 };
153 
154 
155 
156 
157 } // namespace ModelEvaluatorDefaultBaseTypes
158 
159 
187 template<class Scalar>
188 class ModelEvaluatorDefaultBase : virtual public ModelEvaluator<Scalar>
189 {
190 public:
191 
194 
196  int Np() const;
198  int Ng() const;
200  RCP<LinearOpBase<Scalar> > create_DfDp_op(int l) const;
202  RCP<LinearOpBase<Scalar> > create_DgDx_dot_op(int j) const;
204  RCP<LinearOpBase<Scalar> > create_DgDx_op(int j) const;
206  RCP<LinearOpBase<Scalar> > create_DgDp_op(int j, int l) const;
208  RCP<LinearOpWithSolveBase<Scalar> > create_W() const;
212  void evalModel(
215  ) const;
216 
218 
219 protected:
220 
223 
232  void initializeDefaultBase();
233 
235 
236 private:
237 
240 
242  virtual RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
243 
245  virtual RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
246 
248  virtual RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
249 
251  virtual RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
252 
254 
257 
259  virtual ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const = 0;
260 
262  virtual void evalModelImpl(
265  ) const = 0;
266 
268 
269 protected:
270 
273 
274 private:
275 
276  // //////////////////////////////
277  // Private tpyes
278 
279  typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
280  DefaultDerivLinearOpSupport;
281 
282  typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
283  DefaultDerivMvAdjointSupport;
284 
285  typedef ModelEvaluatorDefaultBaseTypes::MultiVectorAdjointPair<Scalar>
286  MultiVectorAdjointPair;
287 
288  // //////////////////////////////
289  // Private data members
290 
291  bool isInitialized_;
292 
293  Array<DefaultDerivLinearOpSupport> DfDp_default_op_support_;
294 
295  Array<DefaultDerivLinearOpSupport> DgDx_dot_default_op_support_;
296 
297  Array<DefaultDerivLinearOpSupport> DgDx_default_op_support_;
298 
299  Array<Array<DefaultDerivLinearOpSupport> > DgDp_default_op_support_;
300  Array<Array<DefaultDerivMvAdjointSupport> > DgDp_default_mv_support_;
301 
302  bool default_W_support_;
303 
304  ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
305 
306  // //////////////////////////////
307  // Private member functions
308 
309  void lazyInitializeDefaultBase() const;
310 
311  void assert_l(const int l) const;
312 
313  void assert_j(const int j) const;
314 
315  // //////////////////////////////
316  // Private static functions
317 
318  static DefaultDerivLinearOpSupport
319  determineDefaultDerivLinearOpSupport(
320  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
321  );
322 
323  static RCP<LinearOpBase<Scalar> >
324  createDefaultLinearOp(
325  const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
326  const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
327  const RCP<const VectorSpaceBase<Scalar> > &var_space
328  );
329 
331  updateDefaultLinearOpSupport(
332  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
333  const DefaultDerivLinearOpSupport &defaultLinearOpSupport
334  );
335 
337  getOutArgImplForDefaultLinearOpSupport(
339  const DefaultDerivLinearOpSupport &defaultLinearOpSupport
340  );
341 
342  static DefaultDerivMvAdjointSupport
343  determineDefaultDerivMvAdjointSupport(
344  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
345  const VectorSpaceBase<Scalar> &fnc_space,
346  const VectorSpaceBase<Scalar> &var_space
347  );
348 
350  updateDefaultDerivMvAdjointSupport(
351  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
352  const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
353  );
354 
355 };
356 
357 
358 } // namespace Thyra
359 
360 
361 //
362 // Implementations
363 //
364 
365 
366 #include "Thyra_ModelEvaluatorHelpers.hpp"
367 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
368 #include "Thyra_DetachedMultiVectorView.hpp"
369 #include "Teuchos_Assert.hpp"
370 
371 
372 namespace Thyra {
373 
374 
375 // Overridden from ModelEvaluator
376 
377 
378 template<class Scalar>
380 {
381  lazyInitializeDefaultBase();
382  return prototypeOutArgs_.Np();
383 }
384 
385 
386 template<class Scalar>
388 {
389  lazyInitializeDefaultBase();
390  return prototypeOutArgs_.Ng();
391 }
392 
393 
394 template<class Scalar>
395 RCP<LinearOpWithSolveBase<Scalar> >
397 {
398  lazyInitializeDefaultBase();
399  if (default_W_support_)
400  return this->get_W_factory()->createOp();
401  return Teuchos::null;
402 }
403 
404 
405 template<class Scalar>
406 RCP<LinearOpBase<Scalar> >
408 {
409  lazyInitializeDefaultBase();
410 #ifdef TEUCHOS_DEBUG
411  assert_l(l);
412 #endif
413  const DefaultDerivLinearOpSupport
414  defaultLinearOpSupport = DfDp_default_op_support_[l];
415  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
416  return createDefaultLinearOp(
417  defaultLinearOpSupport,
418  this->get_f_space(),
419  this->get_p_space(l)
420  );
421  }
422  return this->create_DfDp_op_impl(l);
423 }
424 
425 
426 template<class Scalar>
427 RCP<LinearOpBase<Scalar> >
429 {
430  lazyInitializeDefaultBase();
431 #ifdef TEUCHOS_DEBUG
432  assert_j(j);
433 #endif
434  const DefaultDerivLinearOpSupport
435  defaultLinearOpSupport = DgDx_dot_default_op_support_[j];
436  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
437  return createDefaultLinearOp(
438  defaultLinearOpSupport,
439  this->get_g_space(j),
440  this->get_x_space()
441  );
442  }
443  return this->create_DgDx_dot_op_impl(j);
444 }
445 
446 
447 template<class Scalar>
448 RCP<LinearOpBase<Scalar> >
450 {
451  lazyInitializeDefaultBase();
452 #ifdef TEUCHOS_DEBUG
453  assert_j(j);
454 #endif
455  const DefaultDerivLinearOpSupport
456  defaultLinearOpSupport = DgDx_default_op_support_[j];
457  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
458  return createDefaultLinearOp(
459  defaultLinearOpSupport,
460  this->get_g_space(j),
461  this->get_x_space()
462  );
463  }
464  return this->create_DgDx_op_impl(j);
465 }
466 
467 
468 template<class Scalar>
469 RCP<LinearOpBase<Scalar> >
471 {
472  lazyInitializeDefaultBase();
473 #ifdef TEUCHOS_DEBUG
474  assert_j(j);
475  assert_l(l);
476 #endif
477  const DefaultDerivLinearOpSupport
478  defaultLinearOpSupport = DgDp_default_op_support_[j][l];
479  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
480  return createDefaultLinearOp(
481  defaultLinearOpSupport,
482  this->get_g_space(j),
483  this->get_p_space(l)
484  );
485  }
486  return this->create_DgDp_op_impl(j,l);
487 }
488 
489 
490 template<class Scalar>
493 {
494  lazyInitializeDefaultBase();
495  return prototypeOutArgs_;
496 }
497 
498 
499 template<class Scalar>
503  ) const
504 {
505 
506  using Teuchos::outArg;
507  typedef ModelEvaluatorBase MEB;
508 
509  lazyInitializeDefaultBase();
510 
511  const int l_Np = outArgs.Np();
512  const int l_Ng = outArgs.Ng();
513 
514  //
515  // A) Assert that the inArgs and outArgs object match this class!
516  //
517 
518 #ifdef TEUCHOS_DEBUG
519  assertInArgsEvalObjects(*this,inArgs);
520  assertOutArgsEvalObjects(*this,outArgs,&inArgs);
521 #endif
522 
523  //
524  // B) Setup the OutArgs object for the underlying implementation's
525  // evalModelImpl(...) function
526  //
527 
528  MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
529  Array<MultiVectorAdjointPair> DgDp_temp_adjoint_copies;
530 
531  {
532 
533  outArgsImpl.setArgs(outArgs,true);
534 
535  // DfDp(l)
536  if (outArgsImpl.supports(MEB::OUT_ARG_f)) {
537  for ( int l = 0; l < l_Np; ++l ) {
538  const DefaultDerivLinearOpSupport defaultLinearOpSupport =
539  DfDp_default_op_support_[l];
540  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
541  outArgsImpl.set_DfDp( l,
542  getOutArgImplForDefaultLinearOpSupport(
543  outArgs.get_DfDp(l), defaultLinearOpSupport
544  )
545  );
546  }
547  else {
548  // DfDp(l) already set by outArgsImpl.setArgs(...)!
549  }
550  }
551  }
552 
553  // DgDx_dot(j)
554  for ( int j = 0; j < l_Ng; ++j ) {
555  const DefaultDerivLinearOpSupport defaultLinearOpSupport =
556  DgDx_dot_default_op_support_[j];
557  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
558  outArgsImpl.set_DgDx_dot( j,
559  getOutArgImplForDefaultLinearOpSupport(
560  outArgs.get_DgDx_dot(j), defaultLinearOpSupport
561  )
562  );
563  }
564  else {
565  // DgDx_dot(j) already set by outArgsImpl.setArgs(...)!
566  }
567  }
568 
569  // DgDx(j)
570  for ( int j = 0; j < l_Ng; ++j ) {
571  const DefaultDerivLinearOpSupport defaultLinearOpSupport =
572  DgDx_default_op_support_[j];
573  if (defaultLinearOpSupport.provideDefaultLinearOp()) {
574  outArgsImpl.set_DgDx( j,
575  getOutArgImplForDefaultLinearOpSupport(
576  outArgs.get_DgDx(j), defaultLinearOpSupport
577  )
578  );
579  }
580  else {
581  // DgDx(j) already set by outArgsImpl.setArgs(...)!
582  }
583  }
584 
585  // DgDp(j,l)
586  for ( int j = 0; j < l_Ng; ++j ) {
587  const Array<DefaultDerivLinearOpSupport> &DgDp_default_op_support_j =
588  DgDp_default_op_support_[j];
589  const Array<DefaultDerivMvAdjointSupport> &DgDp_default_mv_support_j =
590  DgDp_default_mv_support_[j];
591  for ( int l = 0; l < l_Np; ++l ) {
592  const DefaultDerivLinearOpSupport defaultLinearOpSupport =
593  DgDp_default_op_support_j[l];
594  const DefaultDerivMvAdjointSupport defaultMvAdjointSupport =
595  DgDp_default_mv_support_j[l];
596  MEB::Derivative<Scalar> DgDp_j_l;
597  if (!outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none())
598  DgDp_j_l = outArgs.get_DgDp(j,l);
599  if (
600  defaultLinearOpSupport.provideDefaultLinearOp()
601  && !is_null(DgDp_j_l.getLinearOp())
602  )
603  {
604  outArgsImpl.set_DgDp( j, l,
605  getOutArgImplForDefaultLinearOpSupport(
606  DgDp_j_l, defaultLinearOpSupport
607  )
608  );
609  }
610  else if (
611  defaultMvAdjointSupport.provideDefaultAdjoint()
612  && !is_null(DgDp_j_l.getMultiVector())
613  )
614  {
615  const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv =
616  DgDp_j_l.getMultiVector();
617  if (
618  defaultMvAdjointSupport.mvAdjointCopyOrientation()
619  ==
620  DgDp_j_l.getMultiVectorOrientation()
621  )
622  {
623  // The orientation of the multi-vector is different so we need to
624  // create a temporary copy to pass to evalModelImpl(...) and then
625  // copy it back again!
626  const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv_adj =
627  createMembers(DgDp_j_l_mv->domain(), DgDp_j_l_mv->range()->dim());
628  outArgsImpl.set_DgDp( j, l,
629  MEB::Derivative<Scalar>(
630  DgDp_j_l_mv_adj,
631  getOtherDerivativeMultiVectorOrientation(
632  defaultMvAdjointSupport.mvAdjointCopyOrientation()
633  )
634  )
635  );
636  // Remember these multi-vectors so that we can do the transpose copy
637  // back after the evaluation!
638  DgDp_temp_adjoint_copies.push_back(
639  MultiVectorAdjointPair(DgDp_j_l_mv, DgDp_j_l_mv_adj)
640  );
641  }
642  else {
643  // The form of the multi-vector is supported by evalModelImpl(..)
644  // and is already set on the outArgsImpl object.
645  }
646  }
647  else {
648  // DgDp(j,l) already set by outArgsImpl.setArgs(...)!
649  }
650  }
651  }
652 
653  // W
654  {
655  RCP<LinearOpWithSolveBase<Scalar> > W;
656  if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
657  const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
658  W_factory = this->get_W_factory();
659  // Extract the underlying W_op object (if it exists)
660  RCP<const LinearOpBase<Scalar> > W_op_const;
661  uninitializeOp<Scalar>(*W_factory, W.ptr(), outArg(W_op_const));
662  RCP<LinearOpBase<Scalar> > W_op;
663  if (!is_null(W_op_const)) {
664  // Here we remove the const. This is perfectly safe since
665  // w.r.t. this class, we put a non-const object in there and we can
666  // expect to change that object after the fact. That is our
667  // prerogative.
668  W_op = Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(W_op_const);
669  }
670  else {
671  // The W_op object has not been initialized yet so create it. The
672  // next time through, we should not have to do this!
673  W_op = this->create_W_op();
674  }
675  outArgsImpl.set_W_op(W_op);
676  }
677  }
678 
679  }
680 
681  //
682  // C) Evaluate the underlying model implementation!
683  //
684 
685  this->evalModelImpl( inArgs, outArgsImpl );
686 
687  //
688  // D) Post-process the output arguments
689  //
690 
691  // Do explicit transposes for DgDp(j,l) if needed
692  const int numMvAdjointCopies = DgDp_temp_adjoint_copies.size();
693  for ( int adj_copy_i = 0; adj_copy_i < numMvAdjointCopies; ++adj_copy_i ) {
694  const MultiVectorAdjointPair adjPair =
695  DgDp_temp_adjoint_copies[adj_copy_i];
696  doExplicitMultiVectorAdjoint( *adjPair.mvImplAdjoint, &*adjPair.mvOuter );
697  }
698 
699  // Update W given W_op and W_factory
700  {
701  RCP<LinearOpWithSolveBase<Scalar> > W;
702  if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
703  const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
704  W_factory = this->get_W_factory();
705  W_factory->setOStream(this->getOStream());
706  W_factory->setVerbLevel(this->getVerbLevel());
707  initializeOp<Scalar>(*W_factory, outArgsImpl.get_W_op().getConst(), W.ptr());
708  }
709  }
710 
711 }
712 
713 
714 // protected
715 
716 
717 // Setup functions called by subclasses
718 
719 template<class Scalar>
721 {
722 
723  typedef ModelEvaluatorBase MEB;
724 
725  // In case we throw half way thorugh, set to uninitialized
726  isInitialized_ = false;
727  default_W_support_ = false;
728 
729  //
730  // A) Get the InArgs and OutArgs from the subclass
731  //
732 
733  const MEB::InArgs<Scalar> inArgs = this->createInArgs();
734  const MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
735 
736  //
737  // B) Validate the subclasses InArgs and OutArgs
738  //
739 
740 #ifdef TEUCHOS_DEBUG
741  assertInArgsOutArgsSetup( this->description(), inArgs, outArgsImpl );
742 #endif // TEUCHOS_DEBUG
743 
744  //
745  // C) Set up support for default derivative objects and prototype OutArgs
746  //
747 
748  const int l_Ng = outArgsImpl.Ng();
749  const int l_Np = outArgsImpl.Np();
750 
751  // Set support for all outputs supported in the underly implementation
752  MEB::OutArgsSetup<Scalar> outArgs;
753  outArgs.setModelEvalDescription(this->description());
754  outArgs.set_Np_Ng(l_Np,l_Ng);
755  outArgs.setSupports(outArgsImpl);
756 
757  // DfDp
758  DfDp_default_op_support_.clear();
759  if (outArgs.supports(MEB::OUT_ARG_f)) {
760  for ( int l = 0; l < l_Np; ++l ) {
761  const MEB::DerivativeSupport DfDp_l_impl_support =
762  outArgsImpl.supports(MEB::OUT_ARG_DfDp,l);
763  const DefaultDerivLinearOpSupport DfDp_l_op_support =
764  determineDefaultDerivLinearOpSupport(DfDp_l_impl_support);
765  DfDp_default_op_support_.push_back(DfDp_l_op_support);
766  outArgs.setSupports(
767  MEB::OUT_ARG_DfDp, l,
768  updateDefaultLinearOpSupport(
769  DfDp_l_impl_support, DfDp_l_op_support
770  )
771  );
772  }
773  }
774 
775  // DgDx_dot
776  DgDx_dot_default_op_support_.clear();
777  for ( int j = 0; j < l_Ng; ++j ) {
778  const MEB::DerivativeSupport DgDx_dot_j_impl_support =
779  outArgsImpl.supports(MEB::OUT_ARG_DgDx_dot,j);
780  const DefaultDerivLinearOpSupport DgDx_dot_j_op_support =
781  determineDefaultDerivLinearOpSupport(DgDx_dot_j_impl_support);
782  DgDx_dot_default_op_support_.push_back(DgDx_dot_j_op_support);
783  outArgs.setSupports(
784  MEB::OUT_ARG_DgDx_dot, j,
785  updateDefaultLinearOpSupport(
786  DgDx_dot_j_impl_support, DgDx_dot_j_op_support
787  )
788  );
789  }
790 
791  // DgDx
792  DgDx_default_op_support_.clear();
793  for ( int j = 0; j < l_Ng; ++j ) {
794  const MEB::DerivativeSupport DgDx_j_impl_support =
795  outArgsImpl.supports(MEB::OUT_ARG_DgDx,j);
796  const DefaultDerivLinearOpSupport DgDx_j_op_support =
797  determineDefaultDerivLinearOpSupport(DgDx_j_impl_support);
798  DgDx_default_op_support_.push_back(DgDx_j_op_support);
799  outArgs.setSupports(
800  MEB::OUT_ARG_DgDx, j,
801  updateDefaultLinearOpSupport(
802  DgDx_j_impl_support, DgDx_j_op_support
803  )
804  );
805  }
806 
807  // DgDp
808  DgDp_default_op_support_.clear();
809  DgDp_default_mv_support_.clear();
810  for ( int j = 0; j < l_Ng; ++j ) {
811  DgDp_default_op_support_.push_back(Array<DefaultDerivLinearOpSupport>());
812  DgDp_default_mv_support_.push_back(Array<DefaultDerivMvAdjointSupport>());
813  for ( int l = 0; l < l_Np; ++l ) {
814  const MEB::DerivativeSupport DgDp_j_l_impl_support =
815  outArgsImpl.supports(MEB::OUT_ARG_DgDp,j,l);
816  // LinearOpBase support
817  const DefaultDerivLinearOpSupport DgDp_j_l_op_support =
818  determineDefaultDerivLinearOpSupport(DgDp_j_l_impl_support);
819  DgDp_default_op_support_[j].push_back(DgDp_j_l_op_support);
820  outArgs.setSupports(
821  MEB::OUT_ARG_DgDp, j, l,
822  updateDefaultLinearOpSupport(
823  DgDp_j_l_impl_support, DgDp_j_l_op_support
824  )
825  );
826  // MultiVectorBase
827  const DefaultDerivMvAdjointSupport DgDp_j_l_mv_support =
828  determineDefaultDerivMvAdjointSupport(
829  DgDp_j_l_impl_support, *this->get_g_space(j), *this->get_p_space(l)
830  );
831  DgDp_default_mv_support_[j].push_back(DgDp_j_l_mv_support);
832  outArgs.setSupports(
833  MEB::OUT_ARG_DgDp, j, l,
834  updateDefaultDerivMvAdjointSupport(
835  outArgs.supports(MEB::OUT_ARG_DgDp, j, l),
836  DgDp_j_l_mv_support
837  )
838  );
839  }
840  }
841  // 2007/09/09: rabart: ToDo: Move the above code into a private helper
842  // function!
843 
844  // W (given W_op and W_factory)
845  default_W_support_ = false;
846  if ( outArgsImpl.supports(MEB::OUT_ARG_W_op) && !is_null(this->get_W_factory())
847  && !outArgsImpl.supports(MEB::OUT_ARG_W) )
848  {
849  default_W_support_ = true;
850  outArgs.setSupports(MEB::OUT_ARG_W);
851  outArgs.set_W_properties(outArgsImpl.get_W_properties());
852  }
853 
854  //
855  // D) All done!
856  //
857 
858  prototypeOutArgs_ = outArgs;
859  isInitialized_ = true;
860 
861 }
862 
863 
864 // Private functions with default implementaton to be overridden by subclasses
865 
866 
867 template<class Scalar>
868 RCP<LinearOpBase<Scalar> >
870 {
871  typedef ModelEvaluatorBase MEB;
872  MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
873  TEUCHOS_TEST_FOR_EXCEPTION(
874  outArgs.supports(MEB::OUT_ARG_DfDp,l).supports(MEB::DERIV_LINEAR_OP),
875  std::logic_error,
876  "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
877  " supports the LinearOpBase form of DfDp("<<l<<") (as determined from its"
878  " OutArgs object created by createOutArgsImpl())"
879  " but this function create_DfDp_op_impl(...) has not been overriden"
880  " to create such an object!"
881  );
882  return Teuchos::null;
883 }
884 
885 
886 template<class Scalar>
887 RCP<LinearOpBase<Scalar> >
888 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_dot_op_impl(int j) const
889 {
890  typedef ModelEvaluatorBase MEB;
891  MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
892  TEUCHOS_TEST_FOR_EXCEPTION(
893  outArgs.supports(MEB::OUT_ARG_DgDx_dot,j).supports(MEB::DERIV_LINEAR_OP),
894  std::logic_error,
895  "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
896  " supports the LinearOpBase form of DgDx_dot("<<j<<") (as determined from"
897  " its OutArgs object created by createOutArgsImpl())"
898  " but this function create_DgDx_dot_op_impl(...) has not been overriden"
899  " to create such an object!"
900  );
901  return Teuchos::null;
902 }
903 
904 
905 template<class Scalar>
906 RCP<LinearOpBase<Scalar> >
907 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_op_impl(int j) const
908 {
909  typedef ModelEvaluatorBase MEB;
910  MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
911  TEUCHOS_TEST_FOR_EXCEPTION(
912  outArgs.supports(MEB::OUT_ARG_DgDx,j).supports(MEB::DERIV_LINEAR_OP),
913  std::logic_error,
914  "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
915  " supports the LinearOpBase form of DgDx("<<j<<") (as determined from"
916  " its OutArgs object created by createOutArgsImpl())"
917  " but this function create_DgDx_op_impl(...) has not been overriden"
918  " to create such an object!"
919  );
920  return Teuchos::null;
921 }
922 
923 
924 template<class Scalar>
925 RCP<LinearOpBase<Scalar> >
926 ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op_impl(int j, int l) const
927 {
928  typedef ModelEvaluatorBase MEB;
929  MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
930  TEUCHOS_TEST_FOR_EXCEPTION(
931  outArgs.supports(MEB::OUT_ARG_DgDp,j,l).supports(MEB::DERIV_LINEAR_OP),
932  std::logic_error,
933  "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
934  " supports the LinearOpBase form of DgDp("<<j<<","<<l<<")"
935  " (as determined from its OutArgs object created by createOutArgsImpl())"
936  " but this function create_DgDp_op_impl(...) has not been overriden"
937  " to create such an object!"
938  );
939  return Teuchos::null;
940 }
941 
942 
943 // protected
944 
945 
946 template<class Scalar>
948  :isInitialized_(false), default_W_support_(false)
949 {}
950 
951 
952 // private
953 
954 
955 template<class Scalar>
957 {
958  if (!isInitialized_)
959  const_cast<ModelEvaluatorDefaultBase<Scalar>*>(this)->initializeDefaultBase();
960 }
961 
962 
963 template<class Scalar>
964 void ModelEvaluatorDefaultBase<Scalar>::assert_l(const int l) const
965 {
966  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(l,0,this->Np());
967 }
968 
969 
970 template<class Scalar>
971 void ModelEvaluatorDefaultBase<Scalar>::assert_j(const int j) const
972 {
973  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j,0,this->Ng());
974 }
975 
976 
977 // Private static functions
978 
979 
980 template<class Scalar>
981 ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
982 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivLinearOpSupport(
983  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
984  )
985 {
986  typedef ModelEvaluatorBase MEB;
987  if (
988  (
989  derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
990  ||
991  derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW)
992  )
993  &&
994  !derivSupportImpl.supports(MEB::DERIV_LINEAR_OP)
995  )
996  {
997  return DefaultDerivLinearOpSupport(
998  derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
999  ? MEB::DERIV_MV_BY_COL
1000  : MEB::DERIV_TRANS_MV_BY_ROW
1001  );
1002  }
1003  return DefaultDerivLinearOpSupport();
1004 }
1005 
1006 
1007 template<class Scalar>
1008 RCP<LinearOpBase<Scalar> >
1009 ModelEvaluatorDefaultBase<Scalar>::createDefaultLinearOp(
1010  const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
1011  const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
1012  const RCP<const VectorSpaceBase<Scalar> > &var_space
1013  )
1014 {
1015  using Teuchos::rcp_implicit_cast;
1016  typedef LinearOpBase<Scalar> LOB;
1017  typedef ModelEvaluatorBase MEB;
1018  switch(defaultLinearOpSupport.mvImplOrientation()) {
1019  case MEB::DERIV_MV_BY_COL:
1020  // The MultiVector will do just fine as the LinearOpBase
1021  return createMembers(fnc_space, var_space->dim());
1022  case MEB::DERIV_TRANS_MV_BY_ROW:
1023  // We will have to implicitly transpose the underlying MultiVector
1024  return nonconstAdjoint<Scalar>(
1025  rcp_implicit_cast<LOB>(createMembers(var_space, fnc_space->dim()))
1026  );
1027 #ifdef TEUCHOS_DEBUG
1028  default:
1029  TEUCHOS_TEST_FOR_EXCEPT(true);
1030 #endif
1031  }
1032  return Teuchos::null; // Will never be called!
1033 }
1034 
1035 
1036 template<class Scalar>
1037 ModelEvaluatorBase::DerivativeSupport
1038 ModelEvaluatorDefaultBase<Scalar>::updateDefaultLinearOpSupport(
1039  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1040  const DefaultDerivLinearOpSupport &defaultLinearOpSupport
1041  )
1042 {
1043  typedef ModelEvaluatorBase MEB;
1044  MEB::DerivativeSupport derivSupport = derivSupportImpl;
1045  if (defaultLinearOpSupport.provideDefaultLinearOp())
1046  derivSupport.plus(MEB::DERIV_LINEAR_OP);
1047  return derivSupport;
1048 }
1049 
1050 
1051 template<class Scalar>
1052 ModelEvaluatorBase::Derivative<Scalar>
1053 ModelEvaluatorDefaultBase<Scalar>::getOutArgImplForDefaultLinearOpSupport(
1054  const ModelEvaluatorBase::Derivative<Scalar> &deriv,
1055  const DefaultDerivLinearOpSupport &defaultLinearOpSupport
1056  )
1057 {
1058 
1059  using Teuchos::rcp_dynamic_cast;
1060  typedef ModelEvaluatorBase MEB;
1061  typedef MultiVectorBase<Scalar> MVB;
1062  typedef ScaledAdjointLinearOpBase<Scalar> SALOB;
1063 
1064  // If the derivative is not a LinearOpBase then it it either empty or is a
1065  // MultiVectorBase object, and in either case we just return it since the
1066  // underlying evalModelImpl(...) functions should be able to handle it!
1067  if (is_null(deriv.getLinearOp()))
1068  return deriv;
1069 
1070  // The derivative is LinearOpBase so get out the underlying MultiVectorBase
1071  // object and return its derivative multi-vector form.
1072  switch(defaultLinearOpSupport.mvImplOrientation()) {
1073  case MEB::DERIV_MV_BY_COL: {
1074  return MEB::Derivative<Scalar>(
1075  rcp_dynamic_cast<MVB>(deriv.getLinearOp(),true),
1076  MEB::DERIV_MV_BY_COL
1077  );
1078  }
1079  case MEB::DERIV_TRANS_MV_BY_ROW: {
1080  return MEB::Derivative<Scalar>(
1081  rcp_dynamic_cast<MVB>(
1082  rcp_dynamic_cast<SALOB>(deriv.getLinearOp(),true)->getNonconstOrigOp()
1083  ),
1084  MEB::DERIV_TRANS_MV_BY_ROW
1085  );
1086  }
1087 #ifdef TEUCHOS_DEBUG
1088  default:
1089  TEUCHOS_TEST_FOR_EXCEPT(true);
1090 #endif
1091  }
1092 
1093  return ModelEvaluatorBase::Derivative<Scalar>(); // Should never get here!
1094 
1095 }
1096 
1097 
1098 template<class Scalar>
1099 ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
1100 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivMvAdjointSupport(
1101  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1102  const VectorSpaceBase<Scalar> &fnc_space,
1103  const VectorSpaceBase<Scalar> &var_space
1104  )
1105 {
1106  typedef ModelEvaluatorBase MEB;
1107  // Here we will support the adjoint copy for of a multi-vector if both
1108  // spaces give rise to in-core vectors.
1109  const bool implSupportsMv =
1110  ( derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1111  || derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
1112  const bool implLacksMvOrientSupport =
1113  ( !derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1114  || !derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
1115  const bool bothSpacesHaveInCoreViews =
1116  ( fnc_space.hasInCoreView() && var_space.hasInCoreView() );
1117  if ( implSupportsMv && implLacksMvOrientSupport && bothSpacesHaveInCoreViews ) {
1118  return DefaultDerivMvAdjointSupport(
1119  derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
1120  ? MEB::DERIV_TRANS_MV_BY_ROW
1121  : MEB::DERIV_MV_BY_COL
1122  );
1123  }
1124  // We can't provide an adjoint copy or such a copy is not needed!
1125  return DefaultDerivMvAdjointSupport();
1126 }
1127 
1128 
1129 template<class Scalar>
1130 ModelEvaluatorBase::DerivativeSupport
1131 ModelEvaluatorDefaultBase<Scalar>::updateDefaultDerivMvAdjointSupport(
1132  const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
1133  const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
1134  )
1135 {
1136  typedef ModelEvaluatorBase MEB;
1137  MEB::DerivativeSupport derivSupport = derivSupportImpl;
1138  if (defaultMvAdjointSupport.provideDefaultAdjoint())
1139  derivSupport.plus(defaultMvAdjointSupport.mvAdjointCopyOrientation());
1140  return derivSupport;
1141 }
1142 
1143 
1144 } // namespace Thyra
1145 
1146 
1147 #endif // THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Default base class for concrete model evaluators.
RCP< LinearOpWithSolveBase< Scalar > > get_W() const
Precondition: supports(OUT_ARG_W)==true.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...
RCP< LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
Derivative< Scalar > get_DfDp(int l) const
Precondition: supports(OUT_ARG_DfDp,l)==true.
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
RCP< LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Derivative< Scalar > get_DgDx(int j) const
Precondition: supports(OUT_ARG_DgDx,j)==true.
Abstract interface for objects that represent a space for vectors.
void initializeDefaultBase()
Function called by subclasses to fully initialize this object on any important change.
Derivative< Scalar > get_DgDx_dot(int j) const
Precondition: supports(OUT_ARG_DgDx_dot,j)==true.
void evalModel(const ModelEvaluatorBase::InArgs< Scalar > &inArgs, const ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
ModelEvaluatorBase::OutArgs< Scalar > createOutArgs() const
RCP< LinearOpBase< Scalar > > create_DfDp_op(int l) const
bool supports(EOutArgsMembers arg) const
Determine if an input argument is supported or not.
RCP< LinearOpWithSolveBase< Scalar > > create_W() const
Base class for all linear operators.
int Ng() const
Return the number of axillary response functions g(j)(...) supported (Ng >= 0).
Base subclass for ModelEvaluator that defines some basic types.
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
Determines the forms of a general derivative that are supported.
RCP< MultiVectorBase< Scalar > > createMembers(const RCP< const VectorSpaceBase< Scalar > > &vs, int numMembers, const std::string &label="")
Create a set of vector members (a MultiVectorBase) from the vector space.
RCP< LinearOpBase< Scalar > > create_DgDx_op(int j) const
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...