Anasazi  Version of the Day
AnasaziEpetraAdapter.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
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 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #include "AnasaziEpetraAdapter.hpp"
30 
35 namespace Anasazi {
36 
38  //
39  //--------Anasazi::EpetraMultiVec Implementation-------------------------------
40  //
42 
43  // Construction/Destruction
44 
45  EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array,
46  const int numvecs, const int stride)
47  : Epetra_MultiVector(Epetra_DataAccess::Copy, Map_in, array, stride, numvecs)
48  {
49  }
50 
51 
52  EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs)
53  : Epetra_MultiVector(Map_in, numvecs)
54  {
55  }
56 
57 
58  EpetraMultiVec::EpetraMultiVec(Epetra_DataAccess CV,
59  const Epetra_MultiVector& P_vec,
60  const std::vector<int>& index )
61  : Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size())
62  {
63  }
64 
65 
66  EpetraMultiVec::EpetraMultiVec(const Epetra_MultiVector& P_vec)
67  : Epetra_MultiVector(P_vec)
68  {
69  }
70 
71 
72  //
73  // member functions inherited from Anasazi::MultiVec
74  //
75  //
76  // Simulating a virtual copy constructor. If we could rely on the co-variance
77  // of virtual functions, we could return a pointer to EpetraMultiVec
78  // (the derived type) instead of a pointer to the pure virtual base class.
79  //
80 
81  MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
82  {
83  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
84  return ptr_apv; // safe upcast.
85  }
86  //
87  // the following is a virtual copy constructor returning
88  // a pointer to the pure virtual class. vector values are
89  // copied.
90  //
91 
93  {
94  EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
95  return ptr_apv; // safe upcast
96  }
97 
98 
99  MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
100  {
101  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::Copy, *this, index);
102  return ptr_apv; // safe upcast.
103  }
104 
105 
106  MultiVec<double>* EpetraMultiVec::CloneViewNonConst ( const std::vector<int>& index )
107  {
108  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
109  return ptr_apv; // safe upcast.
110  }
111 
112  const MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) const
113  {
114  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
115  return ptr_apv; // safe upcast.
116  }
117 
118 
119  void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index )
120  {
121  // this should be revisited to e
122  EpetraMultiVec temp_vec(Epetra_DataAccess::View, *this, index);
123 
124  int numvecs = index.size();
125  if ( A.GetNumberVecs() != numvecs ) {
126  std::vector<int> index2( numvecs );
127  for(int i=0; i<numvecs; i++)
128  index2[i] = i;
129  EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
130  TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
131  EpetraMultiVec A_vec(Epetra_DataAccess::View, *tmp_vec, index2);
132  temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
133  }
134  else {
135  temp_vec.MvAddMv( 1.0, A, 0.0, A );
136  }
137  }
138 
139  //-------------------------------------------------------------
140  //
141  // *this <- alpha * A * B + beta * (*this)
142  //
143  //-------------------------------------------------------------
144 
145  void EpetraMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
146  const Teuchos::SerialDenseMatrix<int,double>& B, double beta )
147  {
148  Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
149  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
150 
151  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
152  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
153 
154  TEUCHOS_TEST_FOR_EXCEPTION(
155  Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) != 0,
156  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
157  }
158 
159  //-------------------------------------------------------------
160  //
161  // *this <- alpha * A + beta * B
162  //
163  //-------------------------------------------------------------
164 
165  void EpetraMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A,
166  double beta, const MultiVec<double>& B)
167  {
168  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
169  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
170  EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B));
171  TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
172 
173  TEUCHOS_TEST_FOR_EXCEPTION(
174  Update( alpha, *A_vec, beta, *B_vec, 0.0 ) != 0,
175  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
176  }
177 
178  //-------------------------------------------------------------
179  //
180  // dense B <- alpha * A^T * (*this)
181  //
182  //-------------------------------------------------------------
183 
184  void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
185  Teuchos::SerialDenseMatrix<int,double>& B
186 #ifdef HAVE_ANASAZI_EXPERIMENTAL
187  , ConjType conj
188 #endif
189  ) const
190  {
191  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
192 
193  if (A_vec) {
194  Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
195  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
196 
197  TEUCHOS_TEST_FOR_EXCEPTION(
198  B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
199  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
200  }
201  }
202 
203  //-------------------------------------------------------------
204  //
205  // b[i] = A[i]^T * this[i]
206  //
207  //-------------------------------------------------------------
208 
209  void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
210 #ifdef HAVE_ANASAZI_EXPERIMENTAL
211  , ConjType conj
212 #endif
213  ) const
214  {
215  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
216  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvDot() cast of MultiVec<double> to EpetraMultiVec failed.");
217 
218  if (( (int)b.size() >= A_vec->NumVectors() ) ) {
219  TEUCHOS_TEST_FOR_EXCEPTION(
220  this->Dot( *A_vec, &b[0] ) != 0,
221  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvDot() call to Epetra_MultiVec::Dot() returned a nonzero value.");
222  }
223  }
224 
225  //-------------------------------------------------------------
226  //
227  // this[i] = alpha[i] * this[i]
228  //
229  //-------------------------------------------------------------
230  void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
231  {
232  // Check to make sure the vector is as long as the multivector has columns.
233  int numvecs = this->NumVectors();
234  TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
235  "Anasazi::EpetraMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
236 
237  std::vector<int> tmp_index( 1, 0 );
238  for (int i=0; i<numvecs; i++) {
239  Epetra_MultiVector temp_vec(Epetra_DataAccess::View, *this, &tmp_index[0], 1);
240  TEUCHOS_TEST_FOR_EXCEPTION(
241  temp_vec.Scale( alpha[i] ) != 0,
242  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvScale() call to Epetra_MultiVec::Scale() returned a nonzero value.");
243  tmp_index[0]++;
244  }
245  }
246 
248  //
249  //--------Anasazi::EpetraOp Implementation-------------------------------------
250  //
252 
253  //
254  // AnasaziOperator constructors
255  //
256  EpetraOp::EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op)
257  : Epetra_Op(Op)
258  {
259  }
260 
262  {
263  }
264  //
265  // AnasaziOperator applications
266  //
268  MultiVec<double>& Y ) const
269  {
270  //
271  // This standard operator computes Y = A*X
272  //
273  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
274  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
275  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
276 
277  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
278  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
279 
280  int info = Epetra_Op->Apply( *vec_X, *vec_Y );
281  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
282  "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
283  }
284 
286  //
287  //--------Anasazi::EpetraGenOp Implementation----------------------------------
288  //
290 
291  //
292  // AnasaziOperator constructors
293  //
294 
295  EpetraGenOp::EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
296  const Teuchos::RCP<Epetra_Operator> &MOp,
297  bool isAInverse_)
298  : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp)
299  {
300  }
301 
303  {
304  }
305  //
306  // AnasaziOperator applications
307  //
309  {
310  //
311  // This generalized operator computes Y = A^{-1}*M*X
312  //
313  int info=0;
314  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
315  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
316  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
317  Epetra_MultiVector temp_Y(*vec_Y);
318 
319  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
320  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
321  //
322  // Need to cast away constness because the member function Apply is not declared const.
323  // Change the transpose setting for the operator if necessary and change it back when done.
324  //
325  // Apply M
326  info = Epetra_MOp->Apply( *vec_X, temp_Y );
327  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
328  "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
329  // Apply A or A^{-1}
330  if (isAInverse) {
331  info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
332  }
333  else {
334  info = Epetra_AOp->Apply( temp_Y, *vec_Y );
335  }
336  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
337  "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
338  }
339 
340  int EpetraGenOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
341  {
342  //
343  // This generalized operator computes Y = A^{-1}*M*X
344  //
345  int info=0;
346  Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
347 
348  // Apply M
349  info = Epetra_MOp->Apply( X, temp_Y );
350  if (info!=0) return info;
351 
352  // Apply A or A^{-1}
353  if (isAInverse)
354  info = Epetra_AOp->ApplyInverse( temp_Y, Y );
355  else
356  info = Epetra_AOp->Apply( temp_Y, Y );
357 
358  return info;
359  }
360 
361  int EpetraGenOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
362  {
363  //
364  // This generalized operator computes Y = M^{-1}*A*X
365  //
366  int info=0;
367  Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
368 
369  // Apply A or A^{-1}
370  if (isAInverse)
371  info = Epetra_AOp->Apply( X, temp_Y );
372  else
373  info = Epetra_AOp->ApplyInverse( X, temp_Y );
374  if (info!=0) return info;
375 
376  // Apply M^{-1}
377  info = Epetra_MOp->ApplyInverse( temp_Y, Y );
378 
379  return info;
380  }
381 
383  //
384  //--------Anasazi::EpetraSymOp Implementation----------------------------------
385  //
387 
388  //
389  // AnasaziOperator constructors
390  //
391  EpetraSymOp::EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op,
392  bool isTrans)
393  : Epetra_Op(Op), isTrans_(isTrans)
394  {
395  }
396 
398  {
399  }
400  //
401  // AnasaziOperator applications
402  //
404  MultiVec<double>& Y ) const
405  {
406  int info=0;
407  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
408  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
409  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
410  Epetra_MultiVector* temp_vec = new Epetra_MultiVector(
411  (isTrans_) ? Epetra_Op->OperatorDomainMap()
412  : Epetra_Op->OperatorRangeMap(),
413  vec_X->NumVectors() );
414 
415  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
416  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
417  TEUCHOS_TEST_FOR_EXCEPTION( temp_vec==NULL, std::invalid_argument, "Anasazi::EpetraSymOp::Apply() allocation Epetra_MultiVector failed.");
418  //
419  // Need to cast away constness because the member function Apply
420  // is not declared const.
421  //
422  // Transpose the operator (if isTrans_ = true)
423  if (isTrans_) {
424  info = Epetra_Op->SetUseTranspose( isTrans_ );
425  if (info != 0) {
426  delete temp_vec;
427  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
428  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
429  }
430  }
431  //
432  // Compute A*X or A'*X
433  //
434  info=Epetra_Op->Apply( *vec_X, *temp_vec );
435  if (info!=0) {
436  delete temp_vec;
437  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
438  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
439  }
440  //
441  // Transpose/Un-transpose the operator based on value of isTrans_
442  info=Epetra_Op->SetUseTranspose( !isTrans_ );
443  if (info!=0) {
444  delete temp_vec;
445  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
446  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
447  }
448 
449  // Compute A^T*(A*X) or A*A^T
450  info=Epetra_Op->Apply( *temp_vec, *vec_Y );
451  if (info!=0) {
452  delete temp_vec;
453  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
454  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
455  }
456 
457  // Un-transpose the operator
458  info=Epetra_Op->SetUseTranspose( false );
459  delete temp_vec;
460  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
461  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
462  }
463 
464  int EpetraSymOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
465  {
466  int info=0;
467  Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
468  //
469  // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
470  //
471  // Transpose the operator (if isTrans_ = true)
472  if (isTrans_) {
473  info=Epetra_Op->SetUseTranspose( isTrans_ );
474  if (info!=0) { return info; }
475  }
476  //
477  // Compute A*X or A^T*X
478  //
479  info=Epetra_Op->Apply( X, temp_vec );
480  if (info!=0) { return info; }
481  //
482  // Transpose/Un-transpose the operator based on value of isTrans_
483  info=Epetra_Op->SetUseTranspose( !isTrans_ );
484  if (info!=0) { return info; }
485 
486  // Compute A^T*(A*X) or A*A^T
487  info=Epetra_Op->Apply( temp_vec, Y );
488  if (info!=0) { return info; }
489 
490  // Un-transpose the operator
491  info=Epetra_Op->SetUseTranspose( false );
492  return info;
493  }
494 
495  int EpetraSymOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
496  {
497  int info=0;
498  Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
499  //
500  // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
501  //
502  // Transpose the operator (if isTrans_ = true)
503  if (!isTrans_) {
504  info=Epetra_Op->SetUseTranspose( !isTrans_ );
505  if (info!=0) { return info; }
506  }
507  //
508  // Compute A^{-1}*X or A^{-T}*X
509  //
510  info=Epetra_Op->ApplyInverse( X, temp_vec );
511  if (info!=0) { return info; }
512  //
513  // Transpose/Un-transpose the operator based on value of isTrans_
514  info=Epetra_Op->SetUseTranspose( isTrans_ );
515  if (info!=0) { return info; }
516 
517  // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
518  info=Epetra_Op->Apply( temp_vec, Y );
519  if (info!=0) { return info; }
520 
521  // Un-transpose the operator
522  info=Epetra_Op->SetUseTranspose( false );
523  return info;
524  }
525 
527  //
528  //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
529  //
531 
532  //
533  // Anasazi::Operator constructors
534  //
535  EpetraSymMVOp::EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
536  bool isTrans)
537  : Epetra_MV(MV), isTrans_(isTrans)
538  {
539  if (isTrans)
540  MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
541  else
542  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
543  }
544 
545  //
546  // AnasaziOperator applications
547  //
549  {
550  int info=0;
551  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
552  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
553  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
554 
555  if (isTrans_) {
556 
557  Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
558 
559  /* A'*X */
560  info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
561  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
562  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
563 
564  /* A*(A'*X) */
565  info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
566  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
567  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
568  }
569  else {
570 
571  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
572 
573  /* A*X */
574  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
575  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
576  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
577 
578  /* A'*(A*X) */
579  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
580  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
581  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
582  }
583  }
584 
586  //
587  //--------Anasazi::EpetraWSymMVOp Implementation--------------------------------
588  //
590 
591  //
592  // Anasazi::Operator constructors
593  //
594  EpetraWSymMVOp::EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
595  const Teuchos::RCP<Epetra_Operator> &OP )
596  : Epetra_MV(MV), Epetra_OP(OP)
597  {
598  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
599  Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
600  Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
601  }
602 
603  //
604  // AnasaziOperator applications
605  //
607  {
608  int info=0;
609  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
610  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
611  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
612 
613  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
614 
615  /* WA*X */
616  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
617  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
618  "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
619 
620  /* A'*(WA*X) */
621  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
622  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
623  "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
624  }
625 
627  //
628  //--------Anasazi::EpetraW2SymMVOp Implementation--------------------------------
629  //
631 
632  //
633  // Anasazi::Operator constructors
634  //
635  EpetraW2SymMVOp::EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
636  const Teuchos::RCP<Epetra_Operator> &OP )
637  : Epetra_MV(MV), Epetra_OP(OP)
638  {
639  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
640  Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
641  Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
642  }
643 
644  //
645  // AnasaziOperator applications
646  //
648  {
649  int info=0;
650  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
651  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
652  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
653 
654  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
655 
656  /* WA*X */
657  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
658  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
659  "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
660 
661  /* (WA)'*(WA*X) */
662  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_WMV, temp_vec, 0.0 );
663  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
664  "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
665 
666  }
667 } // end namespace Anasazi
EpetraMultiVecAccessor is an interfaceto allow any Anasazi::MultiVec implementation that is based on ...
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
void MvTransMv(double alpha, const MultiVec< double > &A, Teuchos::SerialDenseMatrix< int, double > &B) const
Compute a dense matrix B through the matrix-matrix multiply .
MultiVec< double > * CloneViewNonConst(const std::vector< int > &index)
Creates a new EpetraMultiVec that shares the selected contents of *this.
void SetBlock(const MultiVec< double > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
EpetraSymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, bool isTrans=false)
Basic constructor for applying operator [default] or .
void MvDot(const MultiVec< double > &A, std::vector< double > &b) const
Compute a vector b where the components are the individual dot-products, i.e. where A[i] is the i-th...
void MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
ConjType
Enumerated types used to specify conjugation arguments.
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraMultiVec that shares the selected contents of *this.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
MultiVec< double > * CloneCopy() const
Creates a new EpetraMultiVec and copies contents of *this into the new vector (deep copy)...
EpetraSymOp(const Teuchos::RCP< Epetra_Operator > &Op, bool isTrans=false)
Basic constructor for applying operator [default] or .
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
MultiVec< double > * Clone(const int numvecs) const
Creates a new empty EpetraMultiVec containing numvecs columns.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
void MvAddMv(double alpha, const MultiVec< double > &A, double beta, const MultiVec< double > &B)
Replace *this with .
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
EpetraW2SymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, const Teuchos::RCP< Epetra_Operator > &OP)
Basic constructor for applying operator .
EpetraWSymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, const Teuchos::RCP< Epetra_Operator > &OP)
Basic constructor for applying operator .
EpetraMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_MultiVector is n...
EpetraGenOp(const Teuchos::RCP< Epetra_Operator > &AOp, const Teuchos::RCP< Epetra_Operator > &MOp, bool isAInverse=true)
Basic constructor for applying operator [default] or .
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
This method takes the Anasazi::MultiVec X and applies the operator to it resulting in the Anasazi::Mu...
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
EpetraMultiVec(const Epetra_BlockMap &Map_in, const int numvecs)
Basic EpetraMultiVec constructor.
EpetraOp(const Teuchos::RCP< Epetra_Operator > &Op)
Basic constructor. Accepts reference-counted pointer to an Epetra_Operator.
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
Exceptions thrown to signal error in operator application.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.