Anasazi Version of the Day
AnasaziEpetraAdapter.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "AnasaziEpetraAdapter.hpp"
00030 
00035 namespace Anasazi {
00036 
00038   //
00039   //--------Anasazi::EpetraMultiVec Implementation-------------------------------
00040   //
00042 
00043   // Construction/Destruction
00044   
00045   EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, 
00046                                  const int numvecs, const int stride)
00047     : Epetra_MultiVector(Copy, Map_in, array, stride, numvecs) 
00048   {
00049   }
00050   
00051   
00052   EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs)
00053     : Epetra_MultiVector(Map_in, numvecs) 
00054   {
00055   }
00056   
00057   
00058   EpetraMultiVec::EpetraMultiVec(Epetra_DataAccess CV, 
00059                                  const Epetra_MultiVector& P_vec, 
00060                                  const std::vector<int>& index )
00061     : Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size())
00062   {
00063   }
00064   
00065   
00066   EpetraMultiVec::EpetraMultiVec(const Epetra_MultiVector& P_vec)
00067     : Epetra_MultiVector(P_vec) 
00068   {
00069   }
00070   
00071   
00072   //
00073   //  member functions inherited from Anasazi::MultiVec
00074   //
00075   //
00076   //  Simulating a virtual copy constructor. If we could rely on the co-variance
00077   //  of virtual functions, we could return a pointer to EpetraMultiVec
00078   //  (the derived type) instead of a pointer to the pure virtual base class.
00079   //
00080   
00081   MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
00082   {
00083     EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
00084     return ptr_apv; // safe upcast.
00085   }
00086   //
00087   //  the following is a virtual copy constructor returning
00088   //  a pointer to the pure virtual class. vector values are
00089   //  copied.
00090   //
00091   
00092   MultiVec<double>* EpetraMultiVec::CloneCopy() const
00093   {
00094     EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
00095     return ptr_apv; // safe upcast
00096   }
00097   
00098   
00099   MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
00100   {
00101     EpetraMultiVec * ptr_apv = new EpetraMultiVec(Copy, *this, index);
00102     return ptr_apv; // safe upcast.
00103   }
00104   
00105   
00106   MultiVec<double>* EpetraMultiVec::CloneViewNonConst ( const std::vector<int>& index ) 
00107   {
00108     EpetraMultiVec * ptr_apv = new EpetraMultiVec(View, *this, index);
00109     return ptr_apv; // safe upcast.
00110   }
00111   
00112   const MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) const
00113   {
00114     EpetraMultiVec * ptr_apv = new EpetraMultiVec(View, *this, index);
00115     return ptr_apv; // safe upcast.
00116   }
00117   
00118   
00119   void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index ) 
00120   {
00121     // this should be revisited to e
00122     EpetraMultiVec temp_vec(View, *this, index);
00123 
00124     int numvecs = index.size();
00125     if ( A.GetNumberVecs() != numvecs ) {
00126       std::vector<int> index2( numvecs );
00127       for(int i=0; i<numvecs; i++)
00128         index2[i] = i;
00129       EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00130       TEUCHOS_TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
00131       EpetraMultiVec A_vec(View, *tmp_vec, index2);
00132       temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
00133     }
00134     else {
00135       temp_vec.MvAddMv( 1.0, A, 0.0, A );
00136     }
00137   }
00138 
00139   //-------------------------------------------------------------
00140   //
00141   // *this <- alpha * A * B + beta * (*this)
00142   //
00143   //-------------------------------------------------------------
00144   
00145   void EpetraMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00146       const Teuchos::SerialDenseMatrix<int,double>& B, double beta ) 
00147   {
00148     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00149     Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00150     
00151     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00152     TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
00153     
00154     TEUCHOS_TEST_FOR_EXCEPTION( 
00155         Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) != 0,
00156         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
00157   }
00158 
00159   //-------------------------------------------------------------
00160   //
00161   // *this <- alpha * A + beta * B
00162   //
00163   //-------------------------------------------------------------
00164   
00165   void EpetraMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A, 
00166                                  double beta, const MultiVec<double>& B) 
00167   {
00168     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00169     TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
00170     EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B)); 
00171     TEUCHOS_TEST_FOR_EXCEPTION( B_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
00172     
00173     TEUCHOS_TEST_FOR_EXCEPTION( 
00174         Update( alpha, *A_vec, beta, *B_vec, 0.0 ) != 0,
00175         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
00176   }
00177 
00178   //-------------------------------------------------------------
00179   //
00180   // dense B <- alpha * A^T * (*this)
00181   //
00182   //-------------------------------------------------------------
00183   
00184   void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
00185                                    Teuchos::SerialDenseMatrix<int,double>& B
00186 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00187                                    , ConjType conj
00188 #endif
00189                                   ) const
00190   {    
00191     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00192     
00193     if (A_vec) {
00194       Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00195       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00196       
00197     TEUCHOS_TEST_FOR_EXCEPTION( 
00198         B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
00199         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
00200     }
00201   }
00202   
00203   //-------------------------------------------------------------
00204   //
00205   // b[i] = A[i]^T * this[i]
00206   // 
00207   //-------------------------------------------------------------
00208   
00209   void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double> & b
00210 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00211                                , ConjType conj
00212 #endif
00213                              ) const
00214   {
00215     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00216     TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvDot() cast of MultiVec<double> to EpetraMultiVec failed.");
00217 
00218     if (( (int)b.size() >= A_vec->NumVectors() ) ) {
00219       TEUCHOS_TEST_FOR_EXCEPTION( 
00220           this->Dot( *A_vec, &b[0] ) != 0,
00221           EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvDot() call to Epetra_MultiVec::Dot() returned a nonzero value.");
00222     }
00223   }
00224 
00225   //-------------------------------------------------------------
00226   //
00227   // this[i] = alpha[i] * this[i]
00228   // 
00229   //-------------------------------------------------------------
00230   void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
00231   {
00232     // Check to make sure the vector is as long as the multivector has columns.
00233     int numvecs = this->NumVectors();
00234     TEUCHOS_TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument, 
00235         "Anasazi::EpetraMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
00236     
00237     std::vector<int> tmp_index( 1, 0 );
00238     for (int i=0; i<numvecs; i++) {
00239       Epetra_MultiVector temp_vec(View, *this, &tmp_index[0], 1);
00240       TEUCHOS_TEST_FOR_EXCEPTION( 
00241           temp_vec.Scale( alpha[i] ) != 0,
00242           EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvScale() call to Epetra_MultiVec::Scale() returned a nonzero value.");
00243       tmp_index[0]++;
00244     }
00245   }
00246 
00248   //
00249   //--------Anasazi::EpetraOp Implementation-------------------------------------
00250   //
00252 
00253   //
00254   // AnasaziOperator constructors
00255   //  
00256   EpetraOp::EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op) 
00257     : Epetra_Op(Op)
00258   {
00259   }
00260   
00261   EpetraOp::~EpetraOp() 
00262   {
00263   }
00264   //
00265   // AnasaziOperator applications
00266   //
00267   void EpetraOp::Apply ( const MultiVec<double>& X, 
00268                          MultiVec<double>& Y ) const 
00269   {
00270     //
00271     // This standard operator computes Y = A*X
00272     //
00273     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00274     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00275     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00276     
00277     TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00278     TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00279 
00280     int info = Epetra_Op->Apply( *vec_X, *vec_Y );
00281     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00282                         "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00283   }
00284 
00286   //
00287   //--------Anasazi::EpetraGenOp Implementation----------------------------------
00288   //
00290 
00291   //
00292   // AnasaziOperator constructors
00293   //
00294   
00295   EpetraGenOp::EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
00296                            const Teuchos::RCP<Epetra_Operator> &MOp,
00297                            bool isAInverse_) 
00298     : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp) 
00299   {
00300   }
00301     
00302   EpetraGenOp::~EpetraGenOp() 
00303   {
00304   }
00305   //
00306   // AnasaziOperator applications
00307   //
00308   void EpetraGenOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00309   {
00310     //
00311     // This generalized operator computes Y = A^{-1}*M*X
00312     //
00313     int info=0;
00314     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00315     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00316     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00317     Epetra_MultiVector temp_Y(*vec_Y); 
00318     
00319     TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00320     TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00321     //
00322     // Need to cast away constness because the member function Apply is not declared const.  
00323     // Change the transpose setting for the operator if necessary and change it back when done.
00324     //
00325     // Apply M
00326     info = Epetra_MOp->Apply( *vec_X, temp_Y );
00327     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00328                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00329     // Apply A or A^{-1}
00330     if (isAInverse) {
00331       info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
00332     }
00333     else {
00334       info = Epetra_AOp->Apply( temp_Y, *vec_Y );
00335     }
00336     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00337                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00338   }
00339   
00340   int EpetraGenOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00341   {
00342     //
00343     // This generalized operator computes Y = A^{-1}*M*X 
00344     //
00345     int info=0;
00346     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00347     
00348     // Apply M
00349     info = Epetra_MOp->Apply( X, temp_Y );
00350     if (info!=0) return info;
00351     
00352     // Apply A or A^{-1}
00353     if (isAInverse)
00354       info = Epetra_AOp->ApplyInverse( temp_Y, Y );
00355     else
00356       info = Epetra_AOp->Apply( temp_Y, Y );
00357 
00358     return info;
00359   }
00360   
00361   int EpetraGenOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00362   {
00363     //
00364     // This generalized operator computes Y = M^{-1}*A*X 
00365     //
00366     int info=0;
00367     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00368     
00369     // Apply A or A^{-1}
00370     if (isAInverse)
00371       info = Epetra_AOp->Apply( X, temp_Y );
00372     else
00373       info = Epetra_AOp->ApplyInverse( X, temp_Y );
00374     if (info!=0) return info;
00375     
00376     // Apply M^{-1}
00377     info = Epetra_MOp->ApplyInverse( temp_Y, Y );
00378     
00379     return info;
00380   }
00381 
00383   //
00384   //--------Anasazi::EpetraSymOp Implementation----------------------------------
00385   //
00387 
00388   //
00389   // AnasaziOperator constructors
00390   //
00391   EpetraSymOp::EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, 
00392                            bool isTrans) 
00393     : Epetra_Op(Op), isTrans_(isTrans)
00394   {
00395   }
00396   
00397   EpetraSymOp::~EpetraSymOp() 
00398   {
00399   }
00400   //
00401   // AnasaziOperator applications
00402   //
00403   void EpetraSymOp::Apply ( const MultiVec<double>& X, 
00404                             MultiVec<double>& Y ) const 
00405   {
00406     int info=0;
00407     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00408     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00409     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00410     Epetra_MultiVector* temp_vec = new Epetra_MultiVector( 
00411         (isTrans_) ? Epetra_Op->OperatorDomainMap() 
00412         : Epetra_Op->OperatorRangeMap(), 
00413         vec_X->NumVectors() );
00414     
00415     TEUCHOS_TEST_FOR_EXCEPTION( vec_X==NULL   , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00416     TEUCHOS_TEST_FOR_EXCEPTION( vec_Y==NULL   , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00417     TEUCHOS_TEST_FOR_EXCEPTION( temp_vec==NULL, std::invalid_argument, "Anasazi::EpetraSymOp::Apply() allocation Epetra_MultiVector failed.");
00418     //
00419     // Need to cast away constness because the member function Apply
00420     // is not declared const.
00421     //
00422     // Transpose the operator (if isTrans_ = true)
00423     if (isTrans_) {
00424       info = Epetra_Op->SetUseTranspose( isTrans_ );
00425       if (info != 0) {
00426         delete temp_vec;
00427         TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError, 
00428                             "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00429       }
00430     }
00431     //
00432     // Compute A*X or A'*X 
00433     //
00434     info=Epetra_Op->Apply( *vec_X, *temp_vec );
00435     if (info!=0) { 
00436       delete temp_vec; 
00437       TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError, 
00438                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00439     }
00440     //
00441     // Transpose/Un-transpose the operator based on value of isTrans_
00442     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00443     if (info!=0) { 
00444       delete temp_vec; 
00445       TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError, 
00446                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00447     }
00448     
00449     // Compute A^T*(A*X) or A*A^T
00450     info=Epetra_Op->Apply( *temp_vec, *vec_Y );
00451     if (info!=0) { 
00452       delete temp_vec; 
00453       TEUCHOS_TEST_FOR_EXCEPTION( true, OperatorError, 
00454                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00455     }
00456     
00457     // Un-transpose the operator
00458     info=Epetra_Op->SetUseTranspose( false );
00459     delete temp_vec;
00460     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00461                         "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00462   }
00463   
00464   int EpetraSymOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00465   {
00466     int info=0;
00467     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00468     //
00469     // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
00470     //
00471     // Transpose the operator (if isTrans_ = true)
00472     if (isTrans_) {
00473       info=Epetra_Op->SetUseTranspose( isTrans_ );
00474       if (info!=0) { return info; }
00475     }
00476     //
00477     // Compute A*X or A^T*X 
00478     //
00479     info=Epetra_Op->Apply( X, temp_vec );
00480     if (info!=0) { return info; }
00481     //
00482     // Transpose/Un-transpose the operator based on value of isTrans_
00483     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00484     if (info!=0) { return info; }
00485     
00486     // Compute A^T*(A*X) or A*A^T
00487     info=Epetra_Op->Apply( temp_vec, Y );
00488     if (info!=0) { return info; }
00489     
00490     // Un-transpose the operator
00491     info=Epetra_Op->SetUseTranspose( false );
00492     return info;
00493   }
00494   
00495   int EpetraSymOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00496   {
00497     int info=0;
00498     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00499     //
00500     // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
00501     //
00502     // Transpose the operator (if isTrans_ = true)
00503     if (!isTrans_) {
00504       info=Epetra_Op->SetUseTranspose( !isTrans_ );
00505       if (info!=0) { return info; }
00506     }
00507     //
00508     // Compute A^{-1}*X or A^{-T}*X 
00509     //
00510     info=Epetra_Op->ApplyInverse( X, temp_vec );
00511     if (info!=0) { return info; }
00512     //
00513     // Transpose/Un-transpose the operator based on value of isTrans_
00514     info=Epetra_Op->SetUseTranspose( isTrans_ );
00515     if (info!=0) { return info; }
00516     
00517     // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
00518     info=Epetra_Op->Apply( temp_vec, Y );
00519     if (info!=0) { return info; }
00520     
00521     // Un-transpose the operator
00522     info=Epetra_Op->SetUseTranspose( false );
00523     return info;
00524   }
00525 
00527   //
00528   //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
00529   //
00531 
00532   //
00533   // Anasazi::Operator constructors
00534   //
00535   EpetraSymMVOp::EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00536                                bool isTrans) 
00537     : Epetra_MV(MV), isTrans_(isTrans)
00538   {
00539     if (isTrans)
00540       MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
00541     else
00542       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00543   }
00544   
00545   //
00546   // AnasaziOperator applications
00547   //
00548   void EpetraSymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00549   {
00550     int info=0;
00551     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00552     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00553     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00554     
00555     if (isTrans_) {
00556       
00557       Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
00558       
00559       /* A'*X */
00560       info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00561       TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00562                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00563       
00564       /* A*(A'*X) */
00565       info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );      
00566       TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00567                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00568     } 
00569     else {
00570       
00571       Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00572       
00573       /* A*X */
00574       info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00575       TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00576                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00577       
00578       /* A'*(A*X) */
00579       info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
00580       TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00581                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00582     }
00583   }
00584 
00586   //
00587   //--------Anasazi::EpetraWSymMVOp Implementation--------------------------------
00588   //
00590 
00591   //
00592   // Anasazi::Operator constructors
00593   //
00594   EpetraWSymMVOp::EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00595                                  const Teuchos::RCP<Epetra_Operator> &OP ) 
00596     : Epetra_MV(MV), Epetra_OP(OP)
00597   {
00598       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00599       Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) ); 
00600       Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
00601   }
00602   
00603   //
00604   // AnasaziOperator applications
00605   //
00606   void EpetraWSymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00607   {
00608     int info=0;
00609     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00610     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00611     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00612     
00613     Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00614       
00615     /* WA*X */
00616     info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
00617     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00618                         "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00619       
00620     /* A'*(WA*X) */
00621     info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
00622     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00623                         "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00624   }
00625   
00627   //
00628   //--------Anasazi::EpetraW2SymMVOp Implementation--------------------------------
00629   //
00631 
00632   //
00633   // Anasazi::Operator constructors
00634   //
00635   EpetraW2SymMVOp::EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00636                                  const Teuchos::RCP<Epetra_Operator> &OP ) 
00637     : Epetra_MV(MV), Epetra_OP(OP)
00638   {
00639       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00640       Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) ); 
00641       Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
00642   }
00643   
00644   //
00645   // AnasaziOperator applications
00646   //
00647   void EpetraW2SymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00648   {
00649     int info=0;
00650     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00651     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00652     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00653     
00654     Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00655       
00656     /* WA*X */
00657     info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
00658     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00659                         "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00660       
00661     /* (WA)'*(WA*X) */
00662     info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_WMV, temp_vec, 0.0 );
00663     TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00664                         "Anasazi::EpetraW2SymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00665     
00666   }
00667 } // end namespace Anasazi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends