EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_ProductOperator.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
43#include "Epetra_Vector.h"
44#include "Epetra_Map.h"
45
46namespace EpetraExt {
47
48// Constructors / initializers / accessors
49
52
54 const int numOp
55 ,const Teuchos::RCP<const Epetra_Operator> op[]
56 ,const Teuchos::ETransp opTrans[]
57 ,const EApplyMode opInverse[]
58 )
59{
60 initialize(numOp,op,opTrans,opInverse);
61}
62
64 const int numOp
65 ,const Teuchos::RCP<const Epetra_Operator> op[]
66 ,const Teuchos::ETransp opTrans[]
67 ,const EApplyMode opInverse[]
68 )
69{
70#ifdef _DEBUG
71 TEUCHOS_TEST_FOR_EXCEPTION(
72 numOp < 1, std::invalid_argument
73 ,"ProductOperator::initialize(...): Error!"
74 );
75 // ToDo: Validate maps for operators!
76#endif // _DEBUG
77 Op_.resize(numOp);
78 Op_trans_.resize(numOp);
79 Op_inverse_.resize(numOp);
80 std::copy( op, op + numOp, Op_.begin() );
81 std::copy( opTrans, opTrans + numOp, Op_trans_.begin() );
82 std::copy( opInverse, opInverse + numOp, Op_inverse_.begin() );
83 UseTranspose_ = false;
84 // Wipe cache vectors so that they will be reset just to be safe
85 range_vecs_.resize(0);
86 domain_vecs_.resize(0);
87}
88
90 int *numOp
91 ,Teuchos::RCP<const Epetra_Operator> op[]
92 ,Teuchos::ETransp opTrans[]
93 ,EApplyMode opInverse[]
94 )
95{
96#ifdef _DEBUG
97 TEUCHOS_TEST_FOR_EXCEPTION(
98 (op != NULL || opTrans != NULL || opInverse!=NULL) && numOp==NULL
99 ,std::invalid_argument
100 ,"ProductOperator::uninitialize(...): Error!"
101 );
102#endif // _DEBUG
103 if(numOp) {
104 *numOp = Op_.size();
105 if(op) std::copy( Op_.begin(), Op_.end(), op );
106 if(opTrans) std::copy( Op_trans_.begin(), Op_trans_.end(), opTrans );
107 if(opInverse) std::copy( Op_inverse_.begin(), Op_inverse_.end(), opInverse );
108 }
109 UseTranspose_ = false;
110 Op_.resize(0);
111 Op_trans_.resize(0);
112 Op_inverse_.resize(0);
113 range_vecs_.resize(0);
114 domain_vecs_.resize(0);
115}
116
118 const int k
119 ,Teuchos::ETransp opTrans
120 ,EApplyMode opInverse
121 ,const Epetra_MultiVector &X_k
123 ) const
124{
125 Epetra_Operator &Op_k = const_cast<Epetra_Operator&>(*Op_[k]); // Okay since we put back UseTranspose!
126 bool oldUseTranspose = Op_k.UseTranspose();
127 Op_k.SetUseTranspose((opTrans==Teuchos::NO_TRANS) != (Op_trans_[k]==Teuchos::NO_TRANS));
128 const bool applyInverse_k = (opInverse==APPLY_MODE_APPLY) != (Op_inverse_[k]==APPLY_MODE_APPLY);
129 const int err = ! applyInverse_k ? Op_[k]->Apply(X_k,*Y_k) : Op_[k]->ApplyInverse(X_k,*Y_k);
130 Op_k.SetUseTranspose(oldUseTranspose);
131 TEUCHOS_TEST_FOR_EXCEPTION(
132 err!=0, std::runtime_error, "ProductOperator::applyConstituent(...): Error,"
133 " Op["<<k<<"]." << (!applyInverse_k?"Apply":"ApplyInverse") << "(...) "
134 "returned err = " << err << " with Op["<<k<<"].UseTranspose() = "<<
135 Op_[k]->UseTranspose() << "!");
136}
137
138// Overridden from Epetra_Operator
139
141{
143 UseTranspose_ = useTranspose;
144 return 0;
145}
146
148{
150 const int numOp = this->num_Op();
151 // Setup the temporary vectors
152 initializeTempVecs(false);
153 // Apply the constituent operators one at a time!
154 if( !UseTranspose_ ) {
155 //
156 // Forward Mat-vec: Y = M * X (See main documenation)
157 //
158 for( int k = numOp-1; k >= 0; --k ) {
159 const Epetra_MultiVector &X_k = ( k==numOp-1 ? X : *range_vecs_[k] );
160 Epetra_MultiVector &Y_k = ( k==0 ? Y : *range_vecs_[k-1] );
161 applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
162 }
163 }
164 else if( UseTranspose_ ) {
165 //
166 // Adjoint Mat-vec: Y = M' * X (See main documentation)
167 //
168 for( int k = 0; k <= numOp-1; ++k ) {
169 const Epetra_MultiVector &X_k = ( k==0 ? X : *domain_vecs_[k-1] );
170 Epetra_MultiVector &Y_k = ( k==numOp-1 ? Y : *domain_vecs_[k] );
171 applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
172 }
173 }
174 return 0;
175}
176
178{
180 const int numOp = this->num_Op();
181 // Setup the temporary vectors
182 initializeTempVecs(true);
183 // Apply the constituent operators one at a time!
184 if( !UseTranspose_ ) {
185 //
186 // Forward Inverse Mat-vec: Y = inv(M) * X (See main documenation)
187 //
188 for( int k = 0; k <= numOp-1; ++k ) {
189 const Epetra_MultiVector &X_k = ( k==0 ? X : *domain_vecs_[k-1] );
190 Epetra_MultiVector &Y_k = ( k==numOp-1 ? Y : *domain_vecs_[k] );
191 applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
192 }
193 }
194 else if( UseTranspose_ ) {
195 //
196 // Adjoint Invese Mat-vec: Y = inv(M') * X (See main documentation)
197 //
198 for( int k = numOp-1; k >= 0; --k ) {
199 const Epetra_MultiVector &X_k = ( k==numOp-1 ? X : *range_vecs_[k] );
200 Epetra_MultiVector &Y_k = ( k==0 ? Y : *range_vecs_[k-1] );
201 applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
202 }
203 }
204 return 0;
205}
206
208{
210 return -1.0;
211}
212
213const char* ProductOperator::Label() const
214{
216 return NULL;
217}
218
220{
222 return UseTranspose_;
223}
224
226{
228 return false;
229}
230
231const Epetra_Comm&
233{
235 return Op_.front()->OperatorRangeMap().Comm();
236}
237
238const Epetra_Map&
240{
242 return ( Op_trans_.back()==Teuchos::NO_TRANS
243 ? Op_.back()->OperatorDomainMap()
244 : Op_.back()->OperatorRangeMap()
245 );
246}
247
248const Epetra_Map&
250{
252 return ( Op_trans_.front()==Teuchos::NO_TRANS
253 ? Op_.front()->OperatorRangeMap()
254 : Op_.front()->OperatorDomainMap()
255 );
256}
257
258// private
259
260void ProductOperator::initializeTempVecs(bool applyInverse) const
261{
262 const int numOp = this->num_Op ();
263 if (numOp > 0) {
264 // FIXME (mfh 24 Mar 2014): I added the parentheses around the ||
265 // below to silence a compiler warning. I'm concerned that the
266 // original author of that code didn't understand that && takes
267 // precedence over ||, but I didn't want to change the meaning of
268 // the original code.
270 && range_vecs_.size () == 0) {
271 //
272 // Forward Mat-vec
273 //
274 // We need to create storage to hold:
275 //
276 // T[k-1] = M[k]*T[k]
277 //
278 // for k = numOp-1...1
279 //
280 // where: T[numOp-1] = X (input vector)
281 //
282 range_vecs_.resize (numOp - 1);
283 for (int k = numOp-1; k >= 1; --k) {
284 range_vecs_[k-1] = Teuchos::rcp (new Epetra_Vector ((Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
285 ? Op_[k]->OperatorRangeMap ()
286 : Op_[k]->OperatorDomainMap ()));
287 }
288 }
289 // FIXME (mfh 24 Mar 2014): I added the parentheses around the ||
290 // below to silence a compiler warning. I'm concerned that the
291 // original author of that code didn't understand that && takes
292 // precedence over ||, but I didn't want to change the meaning of
293 // the original code.
294 else if (((UseTranspose_ && ! applyInverse) || (! UseTranspose_ && applyInverse))
295 && domain_vecs_.size () == 0) {
296 //
297 // Adjoint Mat-vec
298 //
299 // We need to create storage to hold:
300 //
301 // T[k] = M[k]'*T[k-1]
302 //
303 // for k = 0...numOp-2
304 //
305 // where: T[-1] = X (input vector)
306 //
307 domain_vecs_.resize (numOp - 1);
308 for (int k = 0; k <= numOp - 2; ++k) {
309 domain_vecs_[k] = Teuchos::rcp (new Epetra_Vector ((Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
310 ? Op_[k]->OperatorDomainMap ()
311 : Op_[k]->OperatorRangeMap ()));
312 }
313 }
314 }
315}
316
317} // namespace EpetraExt
void uninitialize(int *num_Op, Teuchos::RCP< const Epetra_Operator > Op[], Teuchos::ETransp Op_trans[], EApplyMode p_inverse[])
Set to an uninitialized state and wipe out memory.
int num_Op() const
Return the number of aggregate opeators.
void initializeTempVecs(bool applyInverse) const
const Epetra_Map & OperatorRangeMap() const
const Epetra_Map & OperatorDomainMap() const
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
void initialize(const int num_Op, const Teuchos::RCP< const Epetra_Operator > Op[], const Teuchos::ETransp Op_trans[], const EApplyMode Op_inverse[])
Setup with constituent operators.
void applyConstituent(const int k, Teuchos::ETransp Op_trans, EApplyMode Op_inverse, const Epetra_MultiVector &X_k, Epetra_MultiVector *Y_k) const
Apply the kth aggregate operator M[k] correctly.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
ProductOperator()
Construct to uninitialized.
virtual bool UseTranspose() const=0
int applyInverse(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Ifpack_CrsRiluk *M, bool verbose)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.