Teko Version of the Day
Loading...
Searching...
No Matches
Teko_InverseFactoryOperator.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_InverseFactoryOperator.hpp"
48
49// Teko includes
50#include "Teko_BasicMappingStrategy.hpp"
51
52// Thyra includes
53#include "Thyra_EpetraLinearOp.hpp"
54
55using Teuchos::RCP;
56using Teuchos::rcp;
57using Teuchos::rcpFromRef;
58using Teuchos::rcp_dynamic_cast;
59
60namespace Teko {
61namespace Epetra {
62
69InverseFactoryOperator::InverseFactoryOperator(const Teuchos::RCP<const InverseFactory> & ifp)
70 : inverseFactory_(ifp), firstBuildComplete_(false), setConstFwdOp_(true)
71{
72}
73
87{
88 if(not clearOld)
89 return;
90 invOperator_ = Teuchos::null;
91}
92
105void InverseFactoryOperator::buildInverseOperator(const Teuchos::RCP<const Epetra_Operator> & A,bool clear)
106{
107 Teko_DEBUG_SCOPE("InverseFactoryOperator::buildInverseOperator",10);
108
109 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
110 RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
111
112 // set the mapping strategy
113 SetMapStrategy(rcp(new InverseMappingStrategy(extractMappingStrategy(A))));
114
115 initInverse(clear);
116
117 // actually build the inverse operator
118 invOperator_ = Teko::buildInverse(*inverseFactory_,thyraA);
119
120 SetOperator(invOperator_,false);
121
122 firstBuildComplete_ = true;
123
124 if(setConstFwdOp_)
125 fwdOp_ = A;
126
127 setConstFwdOp_ = true;
128
129 TEUCHOS_ASSERT(invOperator_!=Teuchos::null);
130 TEUCHOS_ASSERT(getForwardOp()!=Teuchos::null);
131 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
132 TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
133 TEUCHOS_ASSERT(firstBuildComplete_==true);
134}
135
136void InverseFactoryOperator::buildInverseOperator(const Teuchos::RCP<Epetra_Operator> & A,bool /* clear */)
137{
138 setConstFwdOp_ = false;
139
140 fwdOp_.initialize(A);
141
142 buildInverseOperator(A.getConst());
143
144 TEUCHOS_ASSERT(setConstFwdOp_==true);
145}
146
159void InverseFactoryOperator::rebuildInverseOperator(const Teuchos::RCP<const Epetra_Operator> & A)
160{
161 Teko_DEBUG_SCOPE("InverseFactoryOperator::rebuildPreconditioner",10);
162
163 // if the inverse hasn't been built yet, rebuild from scratch
164 if(not firstBuildComplete_) {
165 buildInverseOperator(A,false);
166 return;
167 }
168
169 RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
170 Teko::rebuildInverse(*inverseFactory_,thyraA,invOperator_);
171
172 if(setConstFwdOp_)
173 fwdOp_.initialize(A);
174
175 SetOperator(invOperator_,false);
176
177 setConstFwdOp_ = true;
178
179 TEUCHOS_ASSERT(getForwardOp()!=Teuchos::null);
180 TEUCHOS_ASSERT(invOperator_!=Teuchos::null);
181 TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
182 TEUCHOS_ASSERT(firstBuildComplete_==true);
183}
184
185void InverseFactoryOperator::rebuildInverseOperator(const Teuchos::RCP<Epetra_Operator> & A)
186{
187 setConstFwdOp_ = false;
188
189 fwdOp_.initialize(A);
190
191 // build from constant epetra operator
192 rebuildInverseOperator(A.getConst());
193
194 TEUCHOS_ASSERT(setConstFwdOp_==true);
195}
196
197Teuchos::RCP<const Thyra::LinearOpBase<double> > InverseFactoryOperator::extractLinearOp(const Teuchos::RCP<const Epetra_Operator> & A) const
198{
199 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
200 const RCP<const EpetraOperatorWrapper> & eow = rcp_dynamic_cast<const EpetraOperatorWrapper>(A);
201
202 // if it is an EpetraOperatorWrapper, then get the Thyra operator
203 if(eow!=Teuchos::null)
204 return eow->getThyraOp();
205
206 // otherwise wrap it up as a thyra operator
207 return Thyra::epetraLinearOp(A);
208}
209
210Teuchos::RCP<const MappingStrategy> InverseFactoryOperator::extractMappingStrategy(const Teuchos::RCP<const Epetra_Operator> & A) const
211{
212 // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
213 const RCP<const EpetraOperatorWrapper> & eow = rcp_dynamic_cast<const EpetraOperatorWrapper>(A);
214
215 // if it is an EpetraOperatorWrapper, then get the Thyra operator
216 if(eow!=Teuchos::null)
217 return eow->getMapStrategy();
218
219 // otherwise wrap it up as a thyra operator
220 RCP<const Epetra_Map> range = rcpFromRef(A->OperatorRangeMap());
221 RCP<const Epetra_Map> domain = rcpFromRef(A->OperatorDomainMap());
222 return rcp(new BasicMappingStrategy(range,domain,A->Comm()));
223}
224
225} // end namespace Epetra
226} // end namespace Teko
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
const RCP< const MappingStrategy > getMapStrategy() const
Get the mapping strategy for this wrapper (translate between Thyra and Epetra)
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.
virtual void buildInverseOperator(const Teuchos::RCP< const Epetra_Operator > &A, bool clear=true)
Build this inverse operator from an Epetra_Operator passed in to this object.
virtual void initInverse(bool clearOld=false)
Build the underlying data structure for the inverse operator.
virtual void rebuildInverseOperator(const Teuchos::RCP< const Epetra_Operator > &A)
Rebuild this inverse from an Epetra_Operator passed in this to object.
Teuchos::RCP< const Epetra_Operator > getForwardOp() const
Flip a mapping strategy object around to give the "inverse" mapping strategy.