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, double * array, 
00046                                  const int numvecs, const int stride)
00047     : Epetra_MultiVector(Copy, Map, array, stride, numvecs) 
00048   {
00049   }
00050   
00051   
00052   EpetraMultiVec::EpetraMultiVec(const Epetra_BlockMap& Map, const int numvecs)
00053     : Epetra_MultiVector(Map, 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::CloneView ( const std::vector<int>& index ) 
00107   {
00108     EpetraMultiVec * ptr_apv = new EpetraMultiVec(View, *this, index);
00109     return ptr_apv; // safe upcast.
00110   }
00111   
00112   
00113   void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index ) 
00114   {
00115     // this should be revisited to e
00116     EpetraMultiVec temp_vec(View, *this, index);
00117 
00118     int numvecs = index.size();
00119     if ( A.GetNumberVecs() != numvecs ) {
00120       std::vector<int> index2( numvecs );
00121       for(int i=0; i<numvecs; i++)
00122         index2[i] = i;
00123       EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00124       TEST_FOR_EXCEPTION( tmp_vec==NULL, std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
00125       EpetraMultiVec A_vec(View, *tmp_vec, index2);
00126       temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
00127     }
00128     else {
00129       temp_vec.MvAddMv( 1.0, A, 0.0, A );
00130     }
00131   }
00132 
00133   //-------------------------------------------------------------
00134   //
00135   // *this <- alpha * A * B + beta * (*this)
00136   //
00137   //-------------------------------------------------------------
00138   
00139   void EpetraMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00140       const Teuchos::SerialDenseMatrix<int,double>& B, double beta ) 
00141   {
00142     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00143     Epetra_MultiVector B_Pvec(Copy, LocalMap, B.values(), B.stride(), B.numCols());
00144     
00145     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00146     TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::SetBlocks() cast of MultiVec<double> to EpetraMultiVec failed.");
00147     
00148     TEST_FOR_EXCEPTION( 
00149         Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) != 0,
00150         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
00151   }
00152 
00153   //-------------------------------------------------------------
00154   //
00155   // *this <- alpha * A + beta * B
00156   //
00157   //-------------------------------------------------------------
00158   
00159   void EpetraMultiVec::MvAddMv ( double alpha , const MultiVec<double>& A, 
00160                                  double beta, const MultiVec<double>& B) 
00161   {
00162     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00163     TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
00164     EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B)); 
00165     TEST_FOR_EXCEPTION( B_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvAddMv() cast of MultiVec<double> to EpetraMultiVec failed.");
00166     
00167     TEST_FOR_EXCEPTION( 
00168         Update( alpha, *A_vec, beta, *B_vec, 0.0 ) != 0,
00169         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvAddMv() call to Epetra_MultiVec::Update() returned a nonzero value.");
00170   }
00171 
00172   //-------------------------------------------------------------
00173   //
00174   // dense B <- alpha * A^T * (*this)
00175   //
00176   //-------------------------------------------------------------
00177   
00178   void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
00179                                    Teuchos::SerialDenseMatrix<int,double>& B
00180 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00181                                    , ConjType conj
00182 #endif
00183                                   ) const
00184   {    
00185     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00186     
00187     if (A_vec) {
00188       Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00189       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00190       
00191     TEST_FOR_EXCEPTION( 
00192         B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
00193         EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
00194     }
00195   }
00196   
00197   //-------------------------------------------------------------
00198   //
00199   // b[i] = A[i]^T * this[i]
00200   // 
00201   //-------------------------------------------------------------
00202   
00203   void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double>* b
00204 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00205                                , ConjType conj
00206 #endif
00207                              ) const
00208   {
00209     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00210     TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraMultiVec::MvDot() cast of MultiVec<double> to EpetraMultiVec failed.");
00211 
00212     if ((b!=NULL) && ( (int)b->size() >= A_vec->NumVectors() ) ) {
00213       TEST_FOR_EXCEPTION( 
00214           this->Dot( *A_vec, &(*b)[0] ) != 0,
00215           EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvDot() call to Epetra_MultiVec::Dot() returned a nonzero value.");
00216     }
00217   }
00218 
00219   //-------------------------------------------------------------
00220   //
00221   // this[i] = alpha[i] * this[i]
00222   // 
00223   //-------------------------------------------------------------
00224   void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
00225   {
00226     // Check to make sure the vector is as long as the multivector has columns.
00227     int numvecs = this->NumVectors();
00228     TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument, 
00229         "Anasazi::EpetraMultiVec::MvScale() alpha argument size was inconsistent with number of vectors in mv.");
00230     
00231     std::vector<int> tmp_index( 1, 0 );
00232     for (int i=0; i<numvecs; i++) {
00233       Epetra_MultiVector temp_vec(View, *this, &tmp_index[0], 1);
00234       TEST_FOR_EXCEPTION( 
00235           temp_vec.Scale( alpha[i] ) != 0,
00236           EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvScale() call to Epetra_MultiVec::Scale() returned a nonzero value.");
00237       tmp_index[0]++;
00238     }
00239   }
00240 
00242   //
00243   //--------Anasazi::EpetraOp Implementation-------------------------------------
00244   //
00246 
00247   //
00248   // AnasaziOperator constructors
00249   //  
00250   EpetraOp::EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op) 
00251     : Epetra_Op(Op)
00252   {
00253   }
00254   
00255   EpetraOp::~EpetraOp() 
00256   {
00257   }
00258   //
00259   // AnasaziOperator applications
00260   //
00261   void EpetraOp::Apply ( const MultiVec<double>& X, 
00262                          MultiVec<double>& Y ) const 
00263   {
00264     //
00265     // This standard operator computes Y = A*X
00266     //
00267     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00268     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00269     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00270     
00271     TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00272     TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00273 
00274     int info = Epetra_Op->Apply( *vec_X, *vec_Y );
00275     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00276                         "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00277   }
00278 
00280   //
00281   //--------Anasazi::EpetraGenOp Implementation----------------------------------
00282   //
00284 
00285   //
00286   // AnasaziOperator constructors
00287   //
00288   
00289   EpetraGenOp::EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
00290                            const Teuchos::RCP<Epetra_Operator> &MOp,
00291                            bool isAInverse_) 
00292     : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp) 
00293   {
00294   }
00295     
00296   EpetraGenOp::~EpetraGenOp() 
00297   {
00298   }
00299   //
00300   // AnasaziOperator applications
00301   //
00302   void EpetraGenOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00303   {
00304     //
00305     // This generalized operator computes Y = A^{-1}*M*X
00306     //
00307     int info=0;
00308     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00309     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00310     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00311     Epetra_MultiVector temp_Y(*vec_Y); 
00312     
00313     TEST_FOR_EXCEPTION( vec_X==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00314     TEST_FOR_EXCEPTION( vec_Y==NULL, std::invalid_argument, "Anasazi::EpetraGenOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00315     //
00316     // Need to cast away constness because the member function Apply is not declared const.  
00317     // Change the transpose setting for the operator if necessary and change it back when done.
00318     //
00319     // Apply M
00320     info = Epetra_MOp->Apply( *vec_X, temp_Y );
00321     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00322                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00323     // Apply A or A^{-1}
00324     if (isAInverse) {
00325       info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
00326     }
00327     else {
00328       info = Epetra_AOp->Apply( temp_Y, *vec_Y );
00329     }
00330     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00331                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00332   }
00333   
00334   int EpetraGenOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00335   {
00336     //
00337     // This generalized operator computes Y = A^{-1}*M*X 
00338     //
00339     int info=0;
00340     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00341     
00342     // Apply M
00343     info = Epetra_MOp->Apply( X, temp_Y );
00344     if (info!=0) return info;
00345     
00346     // Apply A or A^{-1}
00347     if (isAInverse)
00348       info = Epetra_AOp->ApplyInverse( temp_Y, Y );
00349     else
00350       info = Epetra_AOp->Apply( temp_Y, Y );
00351 
00352     return info;
00353   }
00354   
00355   int EpetraGenOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00356   {
00357     //
00358     // This generalized operator computes Y = M^{-1}*A*X 
00359     //
00360     int info=0;
00361     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00362     
00363     // Apply A or A^{-1}
00364     if (isAInverse)
00365       info = Epetra_AOp->Apply( X, temp_Y );
00366     else
00367       info = Epetra_AOp->ApplyInverse( X, temp_Y );
00368     if (info!=0) return info;
00369     
00370     // Apply M^{-1}
00371     info = Epetra_MOp->ApplyInverse( temp_Y, Y );
00372     
00373     return info;
00374   }
00375 
00377   //
00378   //--------Anasazi::EpetraSymOp Implementation----------------------------------
00379   //
00381 
00382   //
00383   // AnasaziOperator constructors
00384   //
00385   EpetraSymOp::EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, 
00386                            bool isTrans) 
00387     : Epetra_Op(Op), isTrans_(isTrans)
00388   {
00389   }
00390   
00391   EpetraSymOp::~EpetraSymOp() 
00392   {
00393   }
00394   //
00395   // AnasaziOperator applications
00396   //
00397   void EpetraSymOp::Apply ( const MultiVec<double>& X, 
00398                             MultiVec<double>& Y ) const 
00399   {
00400     int info=0;
00401     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00402     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00403     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00404     Epetra_MultiVector* temp_vec = new Epetra_MultiVector( 
00405         (isTrans_) ? Epetra_Op->OperatorDomainMap() 
00406         : Epetra_Op->OperatorRangeMap(), 
00407         vec_X->NumVectors() );
00408     
00409     TEST_FOR_EXCEPTION( vec_X==NULL   , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00410     TEST_FOR_EXCEPTION( vec_Y==NULL   , std::invalid_argument, "Anasazi::EpetraSymOp::Apply() cast of MultiVec<double> to Epetra_MultiVector failed.");
00411     TEST_FOR_EXCEPTION( temp_vec==NULL, std::invalid_argument, "Anasazi::EpetraSymOp::Apply() allocation Epetra_MultiVector failed.");
00412     //
00413     // Need to cast away constness because the member function Apply
00414     // is not declared const.
00415     //
00416     // Transpose the operator (if isTrans_ = true)
00417     if (isTrans_) {
00418       info = Epetra_Op->SetUseTranspose( isTrans_ );
00419       if (info != 0) {
00420         delete temp_vec;
00421         TEST_FOR_EXCEPTION( true, OperatorError, 
00422                             "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00423       }
00424     }
00425     //
00426     // Compute A*X or A'*X 
00427     //
00428     info=Epetra_Op->Apply( *vec_X, *temp_vec );
00429     if (info!=0) { 
00430       delete temp_vec; 
00431       TEST_FOR_EXCEPTION( true, OperatorError, 
00432                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00433     }
00434     //
00435     // Transpose/Un-transpose the operator based on value of isTrans_
00436     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00437     if (info!=0) { 
00438       delete temp_vec; 
00439       TEST_FOR_EXCEPTION( true, OperatorError, 
00440                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00441     }
00442     
00443     // Compute A^T*(A*X) or A*A^T
00444     info=Epetra_Op->Apply( *temp_vec, *vec_Y );
00445     if (info!=0) { 
00446       delete temp_vec; 
00447       TEST_FOR_EXCEPTION( true, OperatorError, 
00448                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00449     }
00450     
00451     // Un-transpose the operator
00452     info=Epetra_Op->SetUseTranspose( false );
00453     delete temp_vec;
00454     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00455                         "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00456   }
00457   
00458   int EpetraSymOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00459   {
00460     int info=0;
00461     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00462     //
00463     // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
00464     //
00465     // Transpose the operator (if isTrans_ = true)
00466     if (isTrans_) {
00467       info=Epetra_Op->SetUseTranspose( isTrans_ );
00468       if (info!=0) { return info; }
00469     }
00470     //
00471     // Compute A*X or A^T*X 
00472     //
00473     info=Epetra_Op->Apply( X, temp_vec );
00474     if (info!=0) { return info; }
00475     //
00476     // Transpose/Un-transpose the operator based on value of isTrans_
00477     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00478     if (info!=0) { return info; }
00479     
00480     // Compute A^T*(A*X) or A*A^T
00481     info=Epetra_Op->Apply( temp_vec, Y );
00482     if (info!=0) { return info; }
00483     
00484     // Un-transpose the operator
00485     info=Epetra_Op->SetUseTranspose( false );
00486     return info;
00487   }
00488   
00489   int EpetraSymOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00490   {
00491     int info=0;
00492     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00493     //
00494     // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
00495     //
00496     // Transpose the operator (if isTrans_ = true)
00497     if (!isTrans_) {
00498       info=Epetra_Op->SetUseTranspose( !isTrans_ );
00499       if (info!=0) { return info; }
00500     }
00501     //
00502     // Compute A^{-1}*X or A^{-T}*X 
00503     //
00504     info=Epetra_Op->ApplyInverse( X, temp_vec );
00505     if (info!=0) { return info; }
00506     //
00507     // Transpose/Un-transpose the operator based on value of isTrans_
00508     info=Epetra_Op->SetUseTranspose( isTrans_ );
00509     if (info!=0) { return info; }
00510     
00511     // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
00512     info=Epetra_Op->Apply( temp_vec, Y );
00513     if (info!=0) { return info; }
00514     
00515     // Un-transpose the operator
00516     info=Epetra_Op->SetUseTranspose( false );
00517     return info;
00518   }
00519 
00521   //
00522   //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
00523   //
00525 
00526   //
00527   // Anasazi::Operator constructors
00528   //
00529   EpetraSymMVOp::EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00530                                bool isTrans) 
00531     : Epetra_MV(MV), isTrans_(isTrans)
00532   {
00533     if (isTrans)
00534       MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
00535     else
00536       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00537   }
00538   
00539   //
00540   // AnasaziOperator applications
00541   //
00542   void EpetraSymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00543   {
00544     int info=0;
00545     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00546     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00547     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00548     
00549     if (isTrans_) {
00550       
00551       Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
00552       
00553       /* A'*X */
00554       info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00555       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00556                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00557       
00558       /* A*(A'*X) */
00559       info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );      
00560       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00561                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00562     } 
00563     else {
00564       
00565       Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00566       
00567       /* A*X */
00568       info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00569       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00570                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00571       
00572       /* A'*(A*X) */
00573       info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
00574       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00575                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00576     }
00577   }
00578 
00580   //
00581   //--------Anasazi::EpetraWSymMVOp Implementation--------------------------------
00582   //
00584 
00585   //
00586   // Anasazi::Operator constructors
00587   //
00588   EpetraWSymMVOp::EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00589                                  const Teuchos::RCP<Epetra_Operator> &OP ) 
00590     : Epetra_MV(MV), Epetra_OP(OP)
00591   {
00592       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00593       Epetra_WMV = Teuchos::rcp( new Epetra_MultiVector( *MV_blockmap, Epetra_MV->NumVectors() ) ); 
00594       Epetra_OP->Apply( *Epetra_MV, *Epetra_WMV );
00595   }
00596   
00597   //
00598   // AnasaziOperator applications
00599   //
00600   void EpetraWSymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00601   {
00602     int info=0;
00603     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00604     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00605     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00606     
00607     Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00608       
00609     /* WA*X */
00610     info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_WMV, *vec_X, 0.0 );
00611     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00612                         "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00613       
00614     /* (WA)'*(WA*X) */
00615     info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_WMV, temp_vec, 0.0 );
00616     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00617                         "Anasazi::EpetraWSymMVOp::Apply(): Error returned from Epetra_MultiVector::Multiply()" );
00618     
00619   }
00620 } // end namespace Anasazi

Generated on Tue Oct 20 12:45:18 2009 for Anasazi by doxygen 1.4.7