Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultMultiVectorLinearOpWithSolve_def.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_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
43#define THYRA_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
44
45#include "Thyra_DefaultMultiVectorLinearOpWithSolve_decl.hpp"
46#include "Thyra_DefaultDiagonalLinearOp.hpp"
47#include "Thyra_LinearOpWithSolveBase.hpp"
48#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
49#include "Thyra_DefaultMultiVectorProductVector.hpp"
50#include "Thyra_AssertOp.hpp"
51#include "Teuchos_dyn_cast.hpp"
52
53
54namespace Thyra {
55
56
57// Constructors/initializers/accessors
58
59
60template<class Scalar>
62{}
63
64
65template<class Scalar>
68 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
69 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
70 )
71{
72 validateInitialize(lows,multiVecRange,multiVecDomain);
73 lows_ = lows;
74 multiVecRange_ = multiVecRange;
75 multiVecDomain_ = multiVecDomain;
76}
77
78
79template<class Scalar>
81 const RCP<const LinearOpWithSolveBase<Scalar> > &lows,
82 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
83 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
84 )
85{
86 validateInitialize(lows,multiVecRange,multiVecDomain);
87 lows_ = lows;
88 multiVecRange_ = multiVecRange;
89 multiVecDomain_ = multiVecDomain;
90}
91
92
93template<class Scalar>
96{
97 return lows_.getNonconstObj();
98}
99
100
101template<class Scalar>
104{
105 return lows_.getConstObj();
106}
107
108
109template<class Scalar>
111{
112 lows_.uninitialize();
113 multiVecRange_ = Teuchos::null;
114 multiVecDomain_ = Teuchos::null;
115}
116
117
118// Overridden from LinearOpBase
119
120
121template<class Scalar>
124{
125 return multiVecRange_;
126}
127
128
129template<class Scalar>
132{
133 return multiVecDomain_;
134}
135
136
137template<class Scalar>
140{
141 return Teuchos::null; // ToDo: Implement if needed ???
142}
143
144
145// protected
146
147
148// Overridden from LinearOpBase
149
150
151template<class Scalar>
153 EOpTransp M_trans
154 ) const
155{
156 return Thyra::opSupported(*lows_.getConstObj(),M_trans);
157}
158
159
160template<class Scalar>
162 const EOpTransp M_trans,
163 const MultiVectorBase<Scalar> &XX,
164 const Ptr<MultiVectorBase<Scalar> > &YY,
165 const Scalar alpha,
166 const Scalar beta
167 ) const
168{
169
170 using Teuchos::dyn_cast;
172
173 const Ordinal numCols = XX.domain()->dim();
174
175 for (Ordinal col_j = 0; col_j < numCols; ++col_j) {
176
177 const RCP<const VectorBase<Scalar> > x = XX.col(col_j);
178 const RCP<VectorBase<Scalar> > y = YY->col(col_j);
179
181 X = dyn_cast<const MVPV>(*x).getMultiVector().assert_not_null();
183 Y = dyn_cast<MVPV>(*y).getNonconstMultiVector().assert_not_null();
184
185 Thyra::apply( *lows_.getConstObj(), M_trans, *X, Y.ptr(), alpha, beta );
186
187 }
188
189}
190
191
192// Overridden from LinearOpWithSolveBase
193
194
195template<class Scalar>
196bool
198 EOpTransp M_trans
199 ) const
200{
201 return Thyra::solveSupports(*lows_.getConstObj(),M_trans);
202}
203
204
205template<class Scalar>
206bool
208 EOpTransp M_trans, const SolveMeasureType& solveMeasureType
209 ) const
210{
211 return Thyra::solveSupportsSolveMeasureType(
212 *lows_.getConstObj(),M_trans,solveMeasureType);
213}
214
215
216template<class Scalar>
219 const EOpTransp transp,
220 const MultiVectorBase<Scalar> &BB,
221 const Ptr<MultiVectorBase<Scalar> > &XX,
222 const Ptr<const SolveCriteria<Scalar> > solveCriteria
223 ) const
224{
225
226 using Teuchos::dyn_cast;
227 using Teuchos::outArg;
228 using Teuchos::inOutArg;
230
231 const Ordinal numCols = BB.domain()->dim();
232
233 SolveStatus<Scalar> overallSolveStatus;
234 accumulateSolveStatusInit(outArg(overallSolveStatus));
235
236 for (Ordinal col_j = 0; col_j < numCols; ++col_j) {
237
238 const RCP<const VectorBase<Scalar> > b = BB.col(col_j);
239 const RCP<VectorBase<Scalar> > x = XX->col(col_j);
240
242 B = dyn_cast<const MVPV>(*b).getMultiVector().assert_not_null();
244 X = dyn_cast<MVPV>(*x).getNonconstMultiVector().assert_not_null();
245
246 const SolveStatus<Scalar> solveStatus =
247 Thyra::solve(*lows_.getConstObj(), transp, *B, X.ptr(), solveCriteria);
248
249 accumulateSolveStatus(
250 SolveCriteria<Scalar>(), // Never used
251 solveStatus, inOutArg(overallSolveStatus) );
252
253 }
254
255 return overallSolveStatus;
256
257}
258
259
260// private
261
262
263template<class Scalar>
265 const RCP<const LinearOpWithSolveBase<Scalar> > &lows,
266 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
267 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
268 )
269{
270#ifdef TEUCHOS_DEBUG
271 TEUCHOS_TEST_FOR_EXCEPT(is_null(lows));
272 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecRange));
273 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVecDomain));
274 TEUCHOS_TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() );
275 if (lows->range() != Teuchos::null)
277 "DefaultMultiVectorLinearOpWithSolve<Scalar>::initialize(lows,multiVecRange,multiVecDomain)",
278 *lows->range(), *multiVecRange->getBlock(0) );
279 if (lows->domain() != Teuchos::null)
281 "DefaultMultiVectorLinearOpWithSolve<Scalar>::initialize(lows,multiVecRange,multiVecDomain)",
282 *lows->domain(), *multiVecDomain->getBlock(0) );
283#else
284 (void)lows;
285 (void)multiVecRange;
286 (void)multiVecDomain;
287#endif
288}
289
290
291} // end namespace Thyra
292
293
294#endif // THYRA_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
Ptr< T > ptr() const
const RCP< T > & assert_not_null() const
Implicit concrete LinearOpWithSolveBase subclass that takes a flattended out multi-vector and perform...
RCP< const LinearOpWithSolveBase< Scalar > > getLinearOpWithSolve() const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void nonconstInitialize(const RCP< LinearOpWithSolveBase< Scalar > > &lows, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void initialize(const RCP< const LinearOpWithSolveBase< Scalar > > &lows, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecRange, const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &multiVecDomain)
Standard concrete implementation of a product vector space that creates product vectors fromed implic...
Concrete implementation of a product vector which is really composed out of the columns of a multi-ve...
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Base class for all linear operators that can support a high-level solve operation.
Interface for a collection of column vectors called a multi-vector.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
T_To & dyn_cast(T_From &from)
Simple struct that defines the requested solution criteria for a solve.
Simple struct for the return status from a solve.