Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 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 
42 #include "AnasaziEpetraAdapter.hpp"
43 
48 namespace Anasazi {
49 
51  //
52  //--------Anasazi::EpetraMultiVec Implementation-------------------------------
53  //
55 
56  // Construction/Destruction
57 
58  EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array,
59  const int numvecs, const int stride)
60  : Epetra_MultiVector(Epetra_DataAccess::Copy, Map_in, array, stride, numvecs)
61  {
62  }
63 
64 
65  EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs)
66  : Epetra_MultiVector(Map_in, numvecs)
67  {
68  }
69 
70 
72  const Epetra_MultiVector& P_vec,
73  const std::vector<int>& index )
74  : Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size())
75  {
76  }
77 
78 
80  : Epetra_MultiVector(P_vec)
81  {
82  }
83 
84 
85  //
86  // member functions inherited from Anasazi::MultiVec
87  //
88  //
89  // Simulating a virtual copy constructor. If we could rely on the co-variance
90  // of virtual functions, we could return a pointer to EpetraMultiVec
91  // (the derived type) instead of a pointer to the pure virtual base class.
92  //
93 
94  MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
95  {
96  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
97  return ptr_apv; // safe upcast.
98  }
99  //
100  // the following is a virtual copy constructor returning
101  // a pointer to the pure virtual class. vector values are
102  // copied.
103  //
104 
106  {
107  EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
108  return ptr_apv; // safe upcast
109  }
110 
111 
112  MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
113  {
114  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::Copy, *this, index);
115  return ptr_apv; // safe upcast.
116  }
117 
118 
119  MultiVec<double>* EpetraMultiVec::CloneViewNonConst ( const std::vector<int>& index )
120  {
121  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
122  return ptr_apv; // safe upcast.
123  }
124 
125  const MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) const
126  {
127  EpetraMultiVec * ptr_apv = new EpetraMultiVec(Epetra_DataAccess::View, *this, index);
128  return ptr_apv; // safe upcast.
129  }
130 
131 
132  void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index )
133  {
134  // this should be revisited to e
135  EpetraMultiVec temp_vec(Epetra_DataAccess::View, *this, index);
136 
137  int numvecs = index.size();
138  if ( A.GetNumberVecs() != numvecs ) {
139  std::vector<int> index2( numvecs );
140  for(int i=0; i<numvecs; i++)
141  index2[i] = i;
142  EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
143  TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
144  EpetraMultiVec A_vec(Epetra_DataAccess::View, *tmp_vec, index2);
145  temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
146  }
147  else {
148  temp_vec.MvAddMv( 1.0, A, 0.0, A );
149  }
150  }
151 
152  //-------------------------------------------------------------
153  //
154  // *this <- alpha * A * B + beta * (*this)
155  //
156  //-------------------------------------------------------------
157 
158  void EpetraMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
159  const Teuchos::SerialDenseMatrix<int,double>& B, double beta )
160  {
161  Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
162  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
163 
164  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
165  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
166 
168  Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) != 0,
169  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
170  }
171 
172  //-------------------------------------------------------------
173  //
174  // *this <- alpha * A + beta * B
175  //
176  //-------------------------------------------------------------
177 
178  void EpetraMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A,
179  double beta, const MultiVec<double>& B)
180  {
181  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
182  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
183  EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B));
184  TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
185 
187  Update( alpha, *A_vec, beta, *B_vec, 0.0 ) != 0,
188  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
189  }
190 
191  //-------------------------------------------------------------
192  //
193  // dense B <- alpha * A^T * (*this)
194  //
195  //-------------------------------------------------------------
196 
197  void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
199 #ifdef HAVE_ANASAZI_EXPERIMENTAL
200  , ConjType conj
201 #endif
202  ) const
203  {
204  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
205 
206  if (A_vec) {
207  Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
208  Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
209 
211  B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
212  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
213  }
214  }
215 
216  //-------------------------------------------------------------
217  //
218  // b[i] = A[i]^T * this[i]
219  //
220  //-------------------------------------------------------------
221 
222  void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
223 #ifdef HAVE_ANASAZI_EXPERIMENTAL
224  , ConjType conj
225 #endif
226  ) const
227  {
228  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
229  TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::MvDot() cast of MultiVec<double> to EpetraMultiVec failed.");
230 
231  if (( (int)b.size() >= A_vec->NumVectors() ) ) {
233  this->Dot( *A_vec, &b[0] ) != 0,
234  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvDot() call to Epetra_MultiVec::Dot() returned a nonzero value.");
235  }
236  }
237 
238  //-------------------------------------------------------------
239  //
240  // this[i] = alpha[i] * this[i]
241  //
242  //-------------------------------------------------------------
243  void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
244  {
245  // Check to make sure the vector is as long as the multivector has columns.
246  int numvecs = this->NumVectors();
247  TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
248  "Anasazi::EpetraMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
249 
250  std::vector<int> tmp_index( 1, 0 );
251  for (int i=0; i<numvecs; i++) {
252  Epetra_MultiVector temp_vec(Epetra_DataAccess::View, *this, &tmp_index[0], 1);
254  temp_vec.Scale( alpha[i] ) != 0,
255  EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvScale() call to Epetra_MultiVec::Scale() returned a nonzero value.");
256  tmp_index[0]++;
257  }
258  }
259 
261  //
262  //--------Anasazi::EpetraOp Implementation-------------------------------------
263  //
265 
266  //
267  // AnasaziOperator constructors
268  //
270  : Epetra_Op(Op)
271  {
272  }
273 
275  {
276  }
277  //
278  // AnasaziOperator applications
279  //
281  MultiVec<double>& Y ) const
282  {
283  //
284  // This standard operator computes Y = A*X
285  //
286  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
287  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
288  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
289 
290  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
291  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
292 
293  int info = Epetra_Op->Apply( *vec_X, *vec_Y );
295  "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
296  }
297 
299  //
300  //--------Anasazi::EpetraGenOp Implementation----------------------------------
301  //
303 
304  //
305  // AnasaziOperator constructors
306  //
307 
310  bool isAInverse_)
311  : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp)
312  {
313  }
314 
316  {
317  }
318  //
319  // AnasaziOperator applications
320  //
322  {
323  //
324  // This generalized operator computes Y = A^{-1}*M*X
325  //
326  int info=0;
327  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
328  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
329  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
330  Epetra_MultiVector temp_Y(*vec_Y);
331 
332  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
333  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
334  //
335  // Need to cast away constness because the member function Apply is not declared const.
336  // Change the transpose setting for the operator if necessary and change it back when done.
337  //
338  // Apply M
339  info = Epetra_MOp->Apply( *vec_X, temp_Y );
341  "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
342  // Apply A or A^{-1}
343  if (isAInverse) {
344  info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
345  }
346  else {
347  info = Epetra_AOp->Apply( temp_Y, *vec_Y );
348  }
350  "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
351  }
352 
354  {
355  //
356  // This generalized operator computes Y = A^{-1}*M*X
357  //
358  int info=0;
359  Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
360 
361  // Apply M
362  info = Epetra_MOp->Apply( X, temp_Y );
363  if (info!=0) return info;
364 
365  // Apply A or A^{-1}
366  if (isAInverse)
367  info = Epetra_AOp->ApplyInverse( temp_Y, Y );
368  else
369  info = Epetra_AOp->Apply( temp_Y, Y );
370 
371  return info;
372  }
373 
375  {
376  //
377  // This generalized operator computes Y = M^{-1}*A*X
378  //
379  int info=0;
380  Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors());
381 
382  // Apply A or A^{-1}
383  if (isAInverse)
384  info = Epetra_AOp->Apply( X, temp_Y );
385  else
386  info = Epetra_AOp->ApplyInverse( X, temp_Y );
387  if (info!=0) return info;
388 
389  // Apply M^{-1}
390  info = Epetra_MOp->ApplyInverse( temp_Y, Y );
391 
392  return info;
393  }
394 
396  //
397  //--------Anasazi::EpetraSymOp Implementation----------------------------------
398  //
400 
401  //
402  // AnasaziOperator constructors
403  //
405  bool isTrans)
406  : Epetra_Op(Op), isTrans_(isTrans)
407  {
408  }
409 
411  {
412  }
413  //
414  // AnasaziOperator applications
415  //
417  MultiVec<double>& Y ) const
418  {
419  int info=0;
420  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
421  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
422  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
423  Epetra_MultiVector* temp_vec = new Epetra_MultiVector(
424  (isTrans_) ? Epetra_Op->OperatorDomainMap()
425  : Epetra_Op->OperatorRangeMap(),
426  vec_X->NumVectors() );
427 
428  TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
429  TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
430  TEUCHOS_TEST_FOR_EXCEPTION( temp_vec==NULL, std::invalid_argument, "Anasazi::EpetraSymOp::Apply() allocation Epetra_MultiVector failed.");
431  //
432  // Need to cast away constness because the member function Apply
433  // is not declared const.
434  //
435  // Transpose the operator (if isTrans_ = true)
436  if (isTrans_) {
437  info = Epetra_Op->SetUseTranspose( isTrans_ );
438  if (info != 0) {
439  delete temp_vec;
440  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
441  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
442  }
443  }
444  //
445  // Compute A*X or A'*X
446  //
447  info=Epetra_Op->Apply( *vec_X, *temp_vec );
448  if (info!=0) {
449  delete temp_vec;
450  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
451  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
452  }
453  //
454  // Transpose/Un-transpose the operator based on value of isTrans_
455  info=Epetra_Op->SetUseTranspose( !isTrans_ );
456  if (info!=0) {
457  delete temp_vec;
458  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
459  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
460  }
461 
462  // Compute A^T*(A*X) or A*A^T
463  info=Epetra_Op->Apply( *temp_vec, *vec_Y );
464  if (info!=0) {
465  delete temp_vec;
466  TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError,
467  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
468  }
469 
470  // Un-transpose the operator
471  info=Epetra_Op->SetUseTranspose( false );
472  delete temp_vec;
473  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError,
474  "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
475  }
476 
478  {
479  int info=0;
480  Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
481  //
482  // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
483  //
484  // Transpose the operator (if isTrans_ = true)
485  if (isTrans_) {
486  info=Epetra_Op->SetUseTranspose( isTrans_ );
487  if (info!=0) { return info; }
488  }
489  //
490  // Compute A*X or A^T*X
491  //
492  info=Epetra_Op->Apply( X, temp_vec );
493  if (info!=0) { return info; }
494  //
495  // Transpose/Un-transpose the operator based on value of isTrans_
496  info=Epetra_Op->SetUseTranspose( !isTrans_ );
497  if (info!=0) { return info; }
498 
499  // Compute A^T*(A*X) or A*A^T
500  info=Epetra_Op->Apply( temp_vec, Y );
501  if (info!=0) { return info; }
502 
503  // Un-transpose the operator
504  info=Epetra_Op->SetUseTranspose( false );
505  return info;
506  }
507 
509  {
510  int info=0;
511  Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors());
512  //
513  // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
514  //
515  // Transpose the operator (if isTrans_ = true)
516  if (!isTrans_) {
517  info=Epetra_Op->SetUseTranspose( !isTrans_ );
518  if (info!=0) { return info; }
519  }
520  //
521  // Compute A^{-1}*X or A^{-T}*X
522  //
523  info=Epetra_Op->ApplyInverse( X, temp_vec );
524  if (info!=0) { return info; }
525  //
526  // Transpose/Un-transpose the operator based on value of isTrans_
527  info=Epetra_Op->SetUseTranspose( isTrans_ );
528  if (info!=0) { return info; }
529 
530  // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
531  info=Epetra_Op->Apply( temp_vec, Y );
532  if (info!=0) { return info; }
533 
534  // Un-transpose the operator
535  info=Epetra_Op->SetUseTranspose( false );
536  return info;
537  }
538 
540  //
541  //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
542  //
544 
545  //
546  // Anasazi::Operator constructors
547  //
549  bool isTrans)
550  : Epetra_MV(MV), isTrans_(isTrans)
551  {
552  if (isTrans)
553  MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
554  else
555  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
556  }
557 
558  //
559  // AnasaziOperator applications
560  //
562  {
563  int info=0;
564  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
565  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
566  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
567 
568  if (isTrans_) {
569 
570  Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
571 
572  /* A'*X */
573  info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
575  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
576 
577  /* A*(A'*X) */
578  info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
580  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
581  }
582  else {
583 
584  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
585 
586  /* A*X */
587  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
589  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
590 
591  /* A'*(A*X) */
592  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
594  "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
595  }
596  }
597 
599  //
600  //--------Anasazi::EpetraWSymMVOp Implementation--------------------------------
601  //
603 
604  //
605  // Anasazi::Operator constructors
606  //
608  const Teuchos::RCP<Epetra_Operator> &OP )
609  : Epetra_MV(MV), Epetra_OP(OP)
610  {
611  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
612  Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
613  Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
614  }
615 
616  //
617  // AnasaziOperator applications
618  //
620  {
621  int info=0;
622  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
623  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
624  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
625 
626  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
627 
628  /* WA*X */
629  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
631  "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
632 
633  /* A'*(WA*X) */
634  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
636  "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
637  }
638 
640  //
641  //--------Anasazi::EpetraW2SymMVOp Implementation--------------------------------
642  //
644 
645  //
646  // Anasazi::Operator constructors
647  //
649  const Teuchos::RCP<Epetra_Operator> &OP )
650  : Epetra_MV(MV), Epetra_OP(OP)
651  {
652  MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
653  Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) );
654  Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
655  }
656 
657  //
658  // AnasaziOperator applications
659  //
661  {
662  int info=0;
663  MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
664  Epetra_MultiVector* vec_X = dynamic_cast<EpetraMultiVecAccessor*>(&temp_X)->GetEpetraMultiVec();
665  Epetra_MultiVector* vec_Y = dynamic_cast<EpetraMultiVecAccessor*>(&Y)->GetEpetraMultiVec();
666 
667  Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
668 
669  /* WA*X */
670  info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
672  "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
673 
674  /* (WA)'*(WA*X) */
675  info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_WMV, temp_vec, 0.0 );
677  "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
678 
679  }
680 } // end namespace Anasazi
ScalarType * values() const
EpetraMultiVecAccessor is an interfaceto allow any Anasazi::MultiVec implementation that is based on ...
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.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int SetUseTranspose(bool UseTranspose)=0
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].
void MvTransMv(double alpha, const MultiVec< double > &A, Teuchos::SerialDenseMatrix< int, double > &B) const
Compute a dense matrix B through the matrix-matrix multiply .
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual const Epetra_Map & OperatorDomainMap() const =0
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
EpetraSymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, bool isTrans=false)
Basic constructor for applying operator [default] or .
virtual const Epetra_Map & OperatorRangeMap() const =0
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.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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...
EpetraSymOp(const Teuchos::RCP< Epetra_Operator > &Op, bool isTrans=false)
Basic constructor for applying operator [default] or .
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.
MultiVec< double > * Clone(const int numvecs) const
Creates a new empty EpetraMultiVec containing numvecs columns.
OrdinalType numCols() const
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
EpetraW2SymMVOp(const Teuchos::RCP< const Epetra_MultiVector > &MV, const Teuchos::RCP< Epetra_Operator > &OP)
Basic constructor for applying operator .
MultiVec< double > * CloneCopy() const
Creates a new EpetraMultiVec and copies contents of *this into the new vector (deep copy)...
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 .
Epetra_DataAccess
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...
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Apply inverse method [inherited from Epetra_Operator class].
OrdinalType stride() const
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...
OrdinalType numRows() const
const MultiVec< double > * CloneView(const std::vector< int > &index) const
Creates a new EpetraMultiVec that shares the selected contents of *this.
Exceptions thrown to signal error in operator application.
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.
void Apply(const MultiVec< double > &X, MultiVec< double > &Y) const
Apply method [inherited from Anasazi::Operator class].