Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Thyra_MultiVectorLinearOp.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 Thyra_MultiVectorLinearOp_hpp
10#define Thyra_MultiVectorLinearOp_hpp
11
12#include "Thyra_RowStatLinearOpBase.hpp"
13#include "Thyra_ScaledLinearOpBase.hpp"
14#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
15#include "Thyra_DefaultMultiVectorProductVector.hpp"
16#include "Teuchos_ConstNonconstObjectContainer.hpp"
17#include "Thyra_VectorStdOps.hpp"
18#include "Thyra_MultiVectorStdOps.hpp"
19#include "Thyra_AssertOp.hpp"
20#include "Teuchos_dyn_cast.hpp"
21
22namespace Thyra {
23
28template<class Scalar>
30 virtual public RowStatLinearOpBase<Scalar>,
31 virtual public ScaledLinearOpBase<Scalar>
32{
33public:
34
37
40
42 const RCP<LinearOpBase<Scalar> > &op,
43 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
44 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
45 ) {
46 validateInitialize(op,multiVecRange,multiVecDomain);
47 op_ = op;
48 multiVecRange_ = multiVecRange;
49 multiVecDomain_ = multiVecDomain;
50 }
51
53 const RCP<const LinearOpBase<Scalar> > &op,
54 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
55 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
56 ) {
57 validateInitialize(op,multiVecRange,multiVecDomain);
58 op_ = op;
59 multiVecRange_ = multiVecRange;
60 multiVecDomain_ = multiVecDomain;
61 }
62
63 RCP<LinearOpBase<Scalar> >
64 getNonconstLinearOp() { return op_.getNonconstObj(); }
65
66 RCP<const LinearOpBase<Scalar> >
67 getLinearOp() const { return op_.getConstObj(); }
68
69 void uninitialize() {
70 op_.uninitialize();
71 multiVecRange_ = Teuchos::null;
72 multiVecDomain_ = Teuchos::null;
73 }
74
76
79
80 RCP< const VectorSpaceBase<Scalar> > range() const { return multiVecRange_; }
81
82 RCP< const VectorSpaceBase<Scalar> > domain() const { return multiVecDomain_; }
83
84 RCP<const LinearOpBase<Scalar> > clone() const { return Teuchos::null; // ToDo: Implement if needed ???
85 }
87
88protected:
89
92 bool opSupportedImpl(EOpTransp M_trans) const {
93 return Thyra::opSupported(*op_.getConstObj(),M_trans);
94 }
95
97 const EOpTransp M_trans,
98 const MultiVectorBase<Scalar> &XX,
99 const Ptr<MultiVectorBase<Scalar> > &YY,
100 const Scalar alpha,
101 const Scalar beta
102 ) const {
103
104 using Teuchos::dyn_cast;
105 typedef DefaultMultiVectorProductVector<Scalar> MVPV;
106
107 const Ordinal numCols = XX.domain()->dim();
108
109 for (Ordinal col_j = 0; col_j < numCols; ++col_j) {
110
111 const RCP<const VectorBase<Scalar> > x = XX.col(col_j);
112 const RCP<VectorBase<Scalar> > y = YY->col(col_j);
113
114 RCP<const MultiVectorBase<Scalar> >
115 X = dyn_cast<const MVPV>(*x).getMultiVector().assert_not_null();
116 RCP<MultiVectorBase<Scalar> >
117 Y = dyn_cast<MVPV>(*y).getNonconstMultiVector().assert_not_null();
118
119 Thyra::apply( *op_.getConstObj(), M_trans, *X, Y.ptr(), alpha, beta );
120 }
121
122 }
124
127
130 const RowStatLinearOpBaseUtils::ERowStat rowStat
131 ) const
132 { using Teuchos::rcp_dynamic_cast;
133 return rcp_dynamic_cast<const RowStatLinearOpBase<Scalar> >(op_.getConstObj())->rowStatIsSupported(rowStat); }
134
140 const RowStatLinearOpBaseUtils::ERowStat rowStat,
141 const Ptr<VectorBase<Scalar> > &rowStatVec
142 ) const
143 {
144 TEUCHOS_ASSERT(this->rowStatIsSupported(rowStat));
145
146 // Compute the scaling vector from the underlying operator and assign
147 // it to each column of the scaling multivector. We only use the first
148 // column in the scaling below, but this makes the scaling vector
149 // consistent with a Kronecker-product operator
150 typedef DefaultMultiVectorProductVector<Scalar> MVPV;
151 using Teuchos::dyn_cast;
152 using Teuchos::rcp_dynamic_cast;
153 RCP<MultiVectorBase<Scalar> > rowStatMultiVec =
154 dyn_cast<MVPV>(*rowStatVec).getNonconstMultiVector().assert_not_null();
155 const Ordinal numCols = rowStatMultiVec->domain()->dim();
156 if (numCols > 0) {
157 rcp_dynamic_cast<const RowStatLinearOpBase<Scalar> >(op_.getConstObj())->getRowStat(rowStat, rowStatMultiVec->col(0).ptr());
158 for (Ordinal col=1; col<numCols; ++col) {
159 Thyra::copy(*(rowStatMultiVec->col(0)),
160 rowStatMultiVec->col(col).ptr());
161 }
162 }
163 }
164
166
169
170 virtual bool supportsScaleLeftImpl() const {
171 using Teuchos::rcp_dynamic_cast;
172 return rcp_dynamic_cast<const ScaledLinearOpBase<Scalar> >(op_.getConstObj())->supportsScaleLeft();
173 }
174
175 virtual bool supportsScaleRightImpl() const {
176 using Teuchos::rcp_dynamic_cast;
177 return rcp_dynamic_cast<const ScaledLinearOpBase<Scalar> >(op_.getConstObj())->supportsScaleRight();
178 }
179
180 virtual void scaleLeftImpl(const VectorBase<Scalar> &row_scaling) {
181 TEUCHOS_ASSERT(this->supportsScaleLeft());
182
183 // Use the first column in the scaling vector to scale the operator
184 typedef DefaultMultiVectorProductVector<Scalar> MVPV;
185 using Teuchos::dyn_cast;
186 using Teuchos::rcp_dynamic_cast;
187 RCP<const MultiVectorBase<Scalar> > row_scaling_mv =
188 dyn_cast<const MVPV>(row_scaling).getMultiVector().assert_not_null();
189 const Ordinal numCols = row_scaling_mv->domain()->dim();
190 if (numCols > 0) {
191 rcp_dynamic_cast<ScaledLinearOpBase<Scalar> >(op_.getNonconstObj())->scaleLeft(*(row_scaling_mv->col(0)));
192 }
193
194 //row_scaling.describe(std::cout, Teuchos::VERB_EXTREME);
195 }
196
197 virtual void scaleRightImpl(const VectorBase<Scalar> &col_scaling) {
198 TEUCHOS_ASSERT(this->supportsScaleRight());
199
200 // Use the first column in the scaling vector to scale the operator
201 typedef DefaultMultiVectorProductVector<Scalar> MVPV;
202 using Teuchos::dyn_cast;
203 using Teuchos::rcp_dynamic_cast;
204 RCP<const MultiVectorBase<Scalar> > col_scaling_mv =
205 dyn_cast<const MVPV>(col_scaling).getMultiVector().assert_not_null();
206 const Ordinal numCols = col_scaling_mv->domain()->dim();
207 if (numCols > 0) {
208 rcp_dynamic_cast<ScaledLinearOpBase<Scalar> >(op_.getNonconstObj())->scaleRight(*(col_scaling_mv->col(0)));
209 }
210 }
211
213
214private:
215
216 // //////////////////////////////
217 // Private types
218
219 typedef Teuchos::ConstNonconstObjectContainer<LinearOpBase<Scalar> > CNOP;
220
221 // //////////////////////////////
222 // Private data members
223
225 RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > multiVecRange_;
226 RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > multiVecDomain_;
227
228 // //////////////////////////////
229 // Private member functions
230
232 const RCP<const LinearOpBase<Scalar> > &op,
233 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
234 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
235 ) {
236#ifdef TEUCHOS_DEBUG
237 TEUCHOS_TEST_FOR_EXCEPT(is_null(op));
238 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecRange));
239 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecDomain));
240 TEUCHOS_TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() );
241 if (op->range() != Teuchos::null)
242 THYRA_ASSERT_VEC_SPACES(
243 "MultiVectorLinearOp<Scalar>::initialize(op,multiVecRange,multiVecDomain)",
244 *op->range(), *multiVecRange->getBlock(0) );
245 if (op->domain() != Teuchos::null)
246 THYRA_ASSERT_VEC_SPACES(
247 "MultiVectorLinearOp<Scalar>::initialize(op,multiVecRange,multiVecDomain)",
248 *op->domain(), *multiVecDomain->getBlock(0) );
249#else
250 (void)op;
251 (void)multiVecRange;
252 (void)multiVecDomain;
253#endif
254}
255
256};
257
262template<class Scalar>
263RCP<MultiVectorLinearOp<Scalar> >
265{
266 return Teuchos::rcp(new MultiVectorLinearOp<Scalar>());
267}
268
273template<class Scalar>
274RCP<MultiVectorLinearOp<Scalar> >
276 const RCP<LinearOpBase<Scalar> > &op,
277 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
278 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
279 )
280{
281 RCP<MultiVectorLinearOp<Scalar> >
282 mvop = Teuchos::rcp(new MultiVectorLinearOp<Scalar>());
283 mvop->nonconstInitialize(op,multiVecRange,multiVecDomain);
284 return mvop;
285}
286
291template<class Scalar>
292RCP<MultiVectorLinearOp<Scalar> >
294 const RCP<LinearOpBase<Scalar> > &op,
295 const int num_blocks
296 )
297{
298 RCP<MultiVectorLinearOp<Scalar> >
299 mvop = Teuchos::rcp(new MultiVectorLinearOp<Scalar>());
300 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > mv_domain =
301 Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks);
302 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > mv_range =
303 Thyra::multiVectorProductVectorSpace(op->range(), num_blocks);
304 mvop->nonconstInitialize(op,mv_range,mv_domain);
305 return mvop;
306}
307
312template<class Scalar>
313RCP<MultiVectorLinearOp<Scalar> >
315 const RCP<const LinearOpBase<Scalar> > &op,
316 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
317 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
318 )
319{
320 RCP<MultiVectorLinearOp<Scalar> >
321 mvop = Teuchos::rcp(new MultiVectorLinearOp<Scalar>());
322 mvop->initialize(op,multiVecRange,multiVecDomain);
323 return mvop;
324}
325
330template<class Scalar>
331RCP<MultiVectorLinearOp<Scalar> >
333 const RCP<const LinearOpBase<Scalar> > &op,
334 const int num_blocks
335 )
336{
337 RCP<MultiVectorLinearOp<Scalar> >
338 mvop = Teuchos::rcp(new MultiVectorLinearOp<Scalar>());
339 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > mv_domain =
340 Thyra::multiVectorProductVectorSpace(op->domain(), num_blocks);
341 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > mv_range =
342 Thyra::multiVectorProductVectorSpace(op->range(), num_blocks);
343 mvop->initialize(op,mv_range,mv_domain);
344 return mvop;
345}
346
347} // end namespace Thyra
348
349
350#endif
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > multiVecDomain_
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
void initialize(const RCP< const LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
MultiVectorLinearOp()
Construct to uninitialized.
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &XX, const Ptr< MultiVectorBase< Scalar > > &YY, const Scalar alpha, const Scalar beta) const
bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
Determine if a given row stat is supported.
RCP< MultiVectorLinearOp< Scalar > > multiVectorLinearOp(const RCP< const LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
Nonmember constructor function.
void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< Scalar > > &rowStatVec) const
Get some statistics about a supported row.
static void validateInitialize(const RCP< const LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
bool opSupportedImpl(EOpTransp M_trans) const
RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > multiVecRange_
RCP< MultiVectorLinearOp< Scalar > > multiVectorLinearOp()
Nonmember constructor function.
RCP< LinearOpBase< Scalar > > getNonconstLinearOp()
RCP< const LinearOpBase< Scalar > > getLinearOp() const
Teuchos::ConstNonconstObjectContainer< LinearOpBase< Scalar > > CNOP
RCP< const VectorSpaceBase< Scalar > > domain() const
RCP< MultiVectorLinearOp< Scalar > > nonconstMultiVectorLinearOp(const RCP< LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
Nonmember constructor function.
void nonconstInitialize(const RCP< LinearOpBase< Scalar > > &op, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
RCP< MultiVectorLinearOp< Scalar > > nonconstMultiVectorLinearOp(const RCP< LinearOpBase< Scalar > > &op, const int num_blocks)
Nonmember constructor function.
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
RCP< const VectorSpaceBase< Scalar > > range() const
RCP< const LinearOpBase< Scalar > > clone() const
RCP< MultiVectorLinearOp< Scalar > > multiVectorLinearOp(const RCP< const LinearOpBase< Scalar > > &op, const int num_blocks)
Nonmember constructor function.