Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_AdjointSensitivityModelEvaluator_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_AdjointSensitivityModelEvaluator_impl_hpp
10 #define Tempus_AdjointSensitivityModelEvaluator_impl_hpp
11 
12 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
13 #include "Thyra_DefaultMultiVectorProductVector.hpp"
16 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
17 #include "Thyra_DefaultAdjointLinearOpWithSolve.hpp"
19 #include "Thyra_VectorStdOps.hpp"
20 #include "Thyra_MultiVectorStdOps.hpp"
21 
22 namespace Tempus {
23 
24 template <typename Scalar>
27  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
28  const Scalar& t_final,
29  const bool is_pseudotransient,
31  model_(model),
32  t_final_(t_final),
33  is_pseudotransient_(is_pseudotransient),
34  mass_matrix_is_computed_(false),
35  t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
36 {
37  typedef Thyra::ModelEvaluatorBase MEB;
38 
39  // Set parameters
42  if (pList != Teuchos::null)
43  *pl = *pList;
45  mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
46  mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
47  p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
48  g_index_ = pl->get<int>("Response Function Index", 0);
49  num_adjoint_ = model_->get_g_space(g_index_)->dim();
50 
51  // We currently do not support a non-constant mass matrix
53  mass_matrix_is_constant_ == false, std::logic_error,
54  "AdjointSensitivityModelEvaluator currently does not support " <<
55  "non-constant mass matrix df/dx_dot!");
56 
58  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
60  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
62  Thyra::multiVectorProductVectorSpace(model_->get_p_space(p_index_),
63  num_adjoint_);
64 
65  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
66  MEB::InArgsSetup<Scalar> inArgs;
67  inArgs.setModelEvalDescription(this->description());
68  inArgs.setSupports(MEB::IN_ARG_x);
69  inArgs.setSupports(MEB::IN_ARG_t);
70  if (me_inArgs.supports(MEB::IN_ARG_x_dot))
71  inArgs.setSupports(MEB::IN_ARG_x_dot);
72  if (me_inArgs.supports(MEB::IN_ARG_alpha))
73  inArgs.setSupports(MEB::IN_ARG_alpha);
74  if (me_inArgs.supports(MEB::IN_ARG_beta))
75  inArgs.setSupports(MEB::IN_ARG_beta);
76 
77  // Support additional parameters for x and xdot
78  inArgs.set_Np(me_inArgs.Np());
79  prototypeInArgs_ = inArgs;
80 
81  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
82  MEB::OutArgsSetup<Scalar> outArgs;
83  outArgs.setModelEvalDescription(this->description());
84  outArgs.set_Np_Ng(me_inArgs.Np(),1);
85  outArgs.setSupports(MEB::OUT_ARG_f);
86  outArgs.setSupports(MEB::OUT_ARG_W_op);
87  prototypeOutArgs_ = outArgs;
88 
89  // ME must support W_op to define adjoint ODE/DAE.
90  // Must support alpha, beta if it suports x_dot
91  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
92  TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_W_op));
93  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
94  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
95  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
96  }
97  MEB::DerivativeSupport dgdx_support =
98  me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
99  MEB::DerivativeSupport dgdp_support =
100  me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
101  TEUCHOS_ASSERT(dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
102  TEUCHOS_ASSERT(dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
103 }
104 
105 template <typename Scalar>
106 void
110 {
111  sh_ = sh;
112  if (is_pseudotransient_)
113  forward_state_ = sh_->getCurrentState();
114  else {
116  forward_state_ = Teuchos::null;
117  }
118 }
119 
120 template <typename Scalar>
123 get_p_space(int p) const
124 {
125  TEUCHOS_ASSERT(p < model_->Np());
126  return model_->get_p_space(p);
127 }
128 
129 template <typename Scalar>
132 get_p_names(int p) const
133 {
134  TEUCHOS_ASSERT(p < model_->Np());
135  return model_->get_p_names(p);
136 }
137 
138 template <typename Scalar>
141 get_x_space() const
142 {
143  return adjoint_space_;
144 }
145 
146 template <typename Scalar>
149 get_f_space() const
150 {
151  return residual_space_;
152 }
153 
154 template <typename Scalar>
157 get_g_space(int j) const
158 {
159  TEUCHOS_ASSERT(j == 0);
160  return response_space_;
161 }
162 
163 template <typename Scalar>
166 create_W_op() const
167 {
168  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = model_->create_W_op();
169  return Thyra::nonconstMultiVectorLinearOp( Thyra::nonconstAdjoint(op),
170  num_adjoint_);
171 }
172 
173 template <typename Scalar>
177 {
178  using Teuchos::RCP;
179  using Teuchos::rcp_dynamic_cast;
181 
182  RCP<const LOWSFB> factory = model_->get_W_factory();
183  if (factory == Teuchos::null)
184  return Teuchos::null; // model_ doesn't support W_factory
185 
186  RCP<const LOWSFB> alowsfb = Thyra::adjointLinearOpWithSolveFactory(factory);
187  return Thyra::multiVectorLinearOpWithSolveFactory(
188  alowsfb, residual_space_, adjoint_space_);
189 }
190 
191 template <typename Scalar>
195 {
196  return prototypeInArgs_;
197 }
198 
199 template <typename Scalar>
203 {
204  typedef Thyra::ModelEvaluatorBase MEB;
205  using Teuchos::RCP;
206  using Teuchos::rcp_dynamic_cast;
207 
208  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
209  MEB::InArgs<Scalar> nominal = this->createInArgs();
210 
211  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
212 
213  // Set initial x, x_dot
214  RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
215  Thyra::assign(x.ptr(), zero);
216  nominal.set_x(x);
217 
218  if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
219  RCP< Thyra::VectorBase<Scalar> > x_dot =
220  Thyra::createMember(*adjoint_space_);
221  Thyra::assign(x_dot.ptr(), zero);
222  nominal.set_x_dot(x_dot);
223  }
224 
225  const int np = model_->Np();
226  for (int i=0; i<np; ++i)
227  nominal.set_p(i, me_nominal.get_p(i));
228 
229  return nominal;
230 }
231 
232 template <typename Scalar>
236 {
237  return prototypeOutArgs_;
238 }
239 
240 template <typename Scalar>
241 void
244  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
245 {
246  typedef Thyra::ModelEvaluatorBase MEB;
247  using Teuchos::RCP;
248  using Teuchos::rcp_dynamic_cast;
249 
250  // Interpolate forward solution at supplied time, reusing previous
251  // interpolation if possible
252  TEUCHOS_ASSERT(sh_ != Teuchos::null);
253  const Scalar t = inArgs.get_t();
254  Scalar forward_t;
255  if (is_pseudotransient_)
256  forward_t = forward_state_->getTime();
257  else {
258  forward_t = t_final_ - t;
259  if (forward_state_ == Teuchos::null || t_interp_ != t) {
260  if (forward_state_ == Teuchos::null)
261  forward_state_ = sh_->interpolateState(forward_t);
262  else
263  sh_->interpolateState(forward_t, forward_state_.get());
264  t_interp_ = t;
265  }
266  }
267 
268  // setup input arguments for model
269  MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
270  me_inArgs.set_x(forward_state_->getX());
271  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
272  if (inArgs.get_x_dot() != Teuchos::null)
273  me_inArgs.set_x_dot(forward_state_->getXDot());
274  else {
275  if (is_pseudotransient_) {
276  // For pseudo-transient, we have to always use the same form of the
277  // residual in order to reuse df/dx, df/dx_dot,..., so we force
278  // the underlying ME to always compute the implicit form with x_dot == 0
279  // if it wasn't provided.
280  if (my_x_dot_ == Teuchos::null) {
281  my_x_dot_ = Thyra::createMember(model_->get_x_space());
282  Thyra::assign(my_x_dot_.ptr(), Scalar(0.0));
283  }
284  me_inArgs.set_x_dot(my_x_dot_);
285  }
286  else {
287  // clear out xdot if it was set in nominalValues to get to ensure we
288  // get the explicit form
289  me_inArgs.set_x_dot(Teuchos::null);
290  }
291  }
292  }
293  if (me_inArgs.supports(MEB::IN_ARG_t))
294  me_inArgs.set_t(forward_t);
295  const int np = me_inArgs.Np();
296  for (int i=0; i<np; ++i)
297  me_inArgs.set_p(i, inArgs.get_p(i));
298 
299  // compute adjoint W == model W
300  // It would be nice to not reevaluate W in the psuedo-transient case, but
301  // it isn't clear how to do this in a clean way. Probably just need to
302  // control that with the nonlinear solver.
303  RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
304  if (op != Teuchos::null) {
305  if (me_inArgs.supports(MEB::IN_ARG_alpha))
306  me_inArgs.set_alpha(inArgs.get_alpha());
307  if (me_inArgs.supports(MEB::IN_ARG_beta))
308  me_inArgs.set_beta(inArgs.get_beta());
309 
310  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
311  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
312  RCP<Thyra::DefaultScaledAdjointLinearOp<Scalar> > adjoint_op =
314  mv_adjoint_op->getNonconstLinearOp(),true);
315  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
316  me_outArgs.set_W_op(adjoint_op->getNonconstOp());
317  model_->evalModel(me_inArgs, me_outArgs);
318  }
319 
320  RCP<Thyra::VectorBase<Scalar> > adjoint_f = outArgs.get_f();
321  RCP<Thyra::VectorBase<Scalar> > adjoint_g = outArgs.get_g(0);
322  RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv;
323  if (adjoint_f != Teuchos::null || adjoint_g != Teuchos::null) {
324  RCP<const Thyra::VectorBase<Scalar> > adjoint_x =
325  inArgs.get_x().assert_not_null();
326  adjoint_x_mv =
327  rcp_dynamic_cast<const DMVPV>(adjoint_x,true)->getMultiVector();
328  }
329 
330  // Compute adjoint residual F(y):
331  // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y - dg/dx^T
332  // * For explict form, F(y) = -df/dx^T*y + dg/dx^T
333  // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
334  // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
335  if (adjoint_f != Teuchos::null) {
336  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
337  rcp_dynamic_cast<DMVPV>(adjoint_f,true)->getNonconstMultiVector();
338 
339  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
340 
341  // dg/dx^T
342  // Don't re-evaluate dg/dx for pseudotransient
343  if (!is_pseudotransient_ || my_dgdx_mv_ == Teuchos::null) {
344  if (my_dgdx_mv_ == Teuchos::null)
345  my_dgdx_mv_ =
346  Thyra::createMembers(model_->get_x_space(),
347  model_->get_g_space(g_index_)->dim());
348  me_outArgs.set_DgDx(g_index_,
349  MEB::Derivative<Scalar>(my_dgdx_mv_,
350  MEB::DERIV_MV_GRADIENT_FORM));
351  model_->evalModel(me_inArgs, me_outArgs);
352  me_outArgs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
353  }
354  Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
355 
356  // Explicit form of the residual F(y) = -df/dx^T*y + dg/dx^T
357  // Don't re-evaluate df/dx for pseudotransient
358  if (!is_pseudotransient_ || my_dfdx_ == Teuchos::null) {
359  if (my_dfdx_ == Teuchos::null)
360  my_dfdx_ = model_->create_W_op();
361  me_outArgs.set_W_op(my_dfdx_);
362  if (me_inArgs.supports(MEB::IN_ARG_alpha))
363  me_inArgs.set_alpha(0.0);
364  if (me_inArgs.supports(MEB::IN_ARG_beta))
365  me_inArgs.set_beta(1.0);
366  model_->evalModel(me_inArgs, me_outArgs);
367  }
368  my_dfdx_->apply(Thyra::CONJTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
369  Scalar(-1.0), Scalar(1.0));
370 
371  // Implicit form residual F(y) df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
372  // using the second scalar argument to apply() to change the explicit term
373  // above.
374  // Don't re-evaluate df/dx_dot for pseudotransient
375  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
376  RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot = inArgs.get_x_dot();
377  if (adjoint_x_dot != Teuchos::null) {
378  RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
379  rcp_dynamic_cast<const DMVPV>(adjoint_x_dot,true)->getMultiVector();
380  if (mass_matrix_is_identity_) {
381  // F = -F + y_dot
382  Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
383  *adjoint_x_dot_mv);
384  }
385  else {
386  if (!is_pseudotransient_ || my_dfdxdot_ == Teuchos::null) {
387  if (my_dfdxdot_ == Teuchos::null)
388  my_dfdxdot_ = model_->create_W_op();
389  if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
390  me_outArgs.set_W_op(my_dfdxdot_);
391  me_inArgs.set_alpha(1.0);
392  me_inArgs.set_beta(0.0);
393  model_->evalModel(me_inArgs, me_outArgs);
394  mass_matrix_is_computed_ = true;
395  }
396  }
397  my_dfdxdot_->apply(Thyra::CONJTRANS, *adjoint_x_dot_mv,
398  adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
399  }
400  }
401  }
402  }
403 
404  // Compute g = dg/dp^T - df/dp^T*y for computing the model parameter term in
405  // the adjoint sensitivity formula.
406  // We don't add pseudotransient logic here because this part is only
407  // evaluated once in that case anyway.
408  if (adjoint_g != Teuchos::null) {
409  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
410  rcp_dynamic_cast<DMVPV>(adjoint_g,true)->getNonconstMultiVector();
411 
412  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
413 
414  // dg/dp
415  MEB::DerivativeSupport dgdp_support =
416  me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
417  if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
418  me_outArgs.set_DgDp(g_index_, p_index_,
419  MEB::Derivative<Scalar>(adjoint_g_mv,
420  MEB::DERIV_MV_GRADIENT_FORM));
421  model_->evalModel(me_inArgs, me_outArgs);
422  }
423  else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
424  const int num_g = model_->get_g_space(g_index_)->dim();
425  const int num_p = model_->get_p_space(p_index_)->dim();
426  RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans =
427  createMembers(model_->get_g_space(g_index_), num_p);
428  me_outArgs.set_DgDp(g_index_, p_index_,
429  MEB::Derivative<Scalar>(dgdp_trans,
430  MEB::DERIV_MV_JACOBIAN_FORM));
431  model_->evalModel(me_inArgs, me_outArgs);
432  Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*adjoint_g_mv);
433  Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
434  for (int i=0; i<num_p; ++i)
435  for (int j=0; j<num_g; ++j)
436  dgdp_view(i,j) = dgdp_trans_view(j,i);
437  }
438  else
439  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
440  "Invalid dg/dp support");
441  me_outArgs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
442 
443  // dg/dp - df/dp^T*y
444  MEB::DerivativeSupport dfdp_support =
445  me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_);
447  if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
448  if (my_dfdp_op_ == Teuchos::null)
449  my_dfdp_op_ = model_->create_DfDp_op(p_index_);
450  me_outArgs.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
451  trans = Thyra::CONJTRANS;
452  }
453  else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
454  if (my_dfdp_mv_ == Teuchos::null)
455  my_dfdp_mv_ = createMembers(model_->get_f_space(),
456  model_->get_p_space(p_index_)->dim());
457  me_outArgs.set_DfDp(p_index_,
458  MEB::Derivative<Scalar>(my_dfdp_mv_,
459  MEB::DERIV_MV_JACOBIAN_FORM));
460  my_dfdp_op_ = my_dfdp_mv_;
461  trans = Thyra::CONJTRANS;
462  }
463  else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
464  if (my_dfdp_mv_ == Teuchos::null)
465  my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
466  model_->get_f_space()->dim());
467  me_outArgs.set_DfDp(p_index_,
468  MEB::Derivative<Scalar>(my_dfdp_mv_,
469  MEB::DERIV_MV_GRADIENT_FORM));
470  my_dfdp_op_ = my_dfdp_mv_;
471  trans = Thyra::CONJ;
472  }
473  else
475  true, std::logic_error, "Invalid df/dp support");
476  model_->evalModel(me_inArgs, me_outArgs);
477  my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
478  Scalar(-1.0), Scalar(1.0));
479  }
480 }
481 
482 template<class Scalar>
486 {
487  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
488  pl->set<int>("Sensitivity Parameter Index", 0);
489  pl->set<int>("Response Function Index", 0);
490  pl->set<bool>("Mass Matrix Is Constant", true);
491  pl->set<bool>("Mass Matrix Is Identity", false);
492  return pl;
493 }
494 
495 } // namespace Tempus
496 
497 #endif
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
AdjointSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Scalar &t_final, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
EOpTransp
RCP< const VectorBase< Scalar > > get_x_dot() const
T & get(const std::string &name, T def_value)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
static magnitudeType rmax()
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra::ModelEvaluatorBase::OutArgs< Scalar > prototypeOutArgs_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
virtual std::string description() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
RCP< const VectorBase< Scalar > > get_p(int l) const