EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
46 namespace EpetraExt {
47 
48 // Constructors / initializers / accessors
49 
51 {}
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
122  ,Epetra_MultiVector *Y_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 
140 int ProductOperator::SetUseTranspose(bool useTranspose)
141 {
142  assertInitialized();
143  UseTranspose_ = useTranspose;
144  return 0;
145 }
146 
148 {
149  assertInitialized();
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 {
179  assertInitialized();
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 {
209  assertInitialized();
210  return -1.0;
211 }
212 
213 const char* ProductOperator::Label() const
214 {
215  assertInitialized();
216  return NULL;
217 }
218 
220 {
221  assertInitialized();
222  return UseTranspose_;
223 }
224 
226 {
227  assertInitialized();
228  return false;
229 }
230 
231 const Epetra_Comm&
233 {
234  assertInitialized();
235  return Op_.front()->OperatorRangeMap().Comm();
236 }
237 
238 const Epetra_Map&
240 {
241  assertInitialized();
242  return ( Op_trans_.back()==Teuchos::NO_TRANS
243  ? Op_.back()->OperatorDomainMap()
244  : Op_.back()->OperatorRangeMap()
245  );
246 }
247 
248 const Epetra_Map&
250 {
251  assertInitialized();
252  return ( Op_trans_.front()==Teuchos::NO_TRANS
253  ? Op_.front()->OperatorRangeMap()
254  : Op_.front()->OperatorDomainMap()
255  );
256 }
257 
258 // private
259 
260 void 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.
269  if (((! UseTranspose_ && ! applyInverse) || (UseTranspose_ && applyInverse))
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
virtual int SetUseTranspose(bool UseTranspose)=0
const Epetra_Map & OperatorDomainMap() 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.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
ProductOperator()
Construct to uninitialized.
const Epetra_Comm & Comm() const
int num_Op() const
Return the number of aggregate opeators.
const Epetra_Map & OperatorRangeMap() const
virtual bool UseTranspose() const =0
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
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.
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.