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, const Epetra_MultiVector& P_vec,         
00059          const std::vector<int>& index )
00060     : Epetra_MultiVector(CV, P_vec, &(const_cast<std::vector<int> &>(index))[0], index.size())
00061   {
00062   }
00063   
00064   
00065   EpetraMultiVec::EpetraMultiVec(const Epetra_MultiVector& P_vec)
00066     : Epetra_MultiVector(P_vec) 
00067   {
00068   }
00069   
00070   
00071   //
00072   //  member functions inherited from Anasazi::MultiVec
00073   //
00074   //
00075   //  Simulating a virtual copy constructor. If we could rely on the co-variance
00076   //  of virtual functions, we could return a pointer to EpetraMultiVec
00077   //  (the derived type) instead of a pointer to the pure virtual base class.
00078   //
00079   
00080   MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
00081   {
00082     EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
00083     return ptr_apv; // safe upcast.
00084   }
00085   //
00086   //  the following is a virtual copy constructor returning
00087   //  a pointer to the pure virtual class. vector values are
00088   //  copied.
00089   //
00090   
00091   MultiVec<double>* EpetraMultiVec::CloneCopy() const
00092   {
00093     EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
00094     return ptr_apv; // safe upcast
00095   }
00096   
00097   
00098   MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
00099   {
00100     EpetraMultiVec * ptr_apv = new EpetraMultiVec(Copy, *this, index);
00101     return ptr_apv; // safe upcast.
00102   }
00103   
00104   
00105   MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) 
00106   {
00107     EpetraMultiVec * ptr_apv = new EpetraMultiVec(View, *this, index);
00108     return ptr_apv; // safe upcast.
00109   }
00110   
00111   
00112   void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index ) 
00113   { 
00114     EpetraMultiVec temp_vec(View, *this, index);
00115 
00116     int numvecs = index.size();
00117     if ( A.GetNumberVecs() != numvecs ) {
00118       std::vector<int> index2( numvecs );
00119       for(int i=0; i<numvecs; i++)
00120   index2[i] = i;
00121       EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00122       assert(tmp_vec!=NULL);
00123       EpetraMultiVec A_vec(View, *tmp_vec, index2);
00124       temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
00125     }
00126     else {
00127       temp_vec.MvAddMv( 1.0, A, 0.0, A );
00128     }
00129   }               
00130 
00131   //-------------------------------------------------------------
00132   //
00133   // *this <- alpha * A * B + beta * (*this)
00134   //
00135   //-------------------------------------------------------------
00136   
00137   void EpetraMultiVec::MvTimesMatAddMv ( const double alpha, const MultiVec<double>& A, 
00138            const Teuchos::SerialDenseMatrix<int,double>& B, const double beta ) 
00139   {
00140     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00141     Epetra_MultiVector B_Pvec(Copy, LocalMap, B.values(), B.stride(), B.numCols());
00142     
00143     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00144     assert(A_vec!=NULL);
00145     
00146     int ret = Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta );
00147     assert( ret == 0 );
00148   }
00149 
00150   //-------------------------------------------------------------
00151   //
00152   // *this <- alpha * A + beta * B
00153   //
00154   //-------------------------------------------------------------
00155   
00156   void EpetraMultiVec::MvAddMv ( const double alpha , const MultiVec<double>& A, 
00157          const double beta, const MultiVec<double>& B) 
00158   {
00159     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00160     assert(A_vec!=NULL);
00161     EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B)); 
00162     assert(B_vec!=NULL);
00163     
00164     int ret = Update( alpha, *A_vec, beta, *B_vec, 0.0 );
00165     assert( ret == 0 ); 
00166   }
00167 
00168   //-------------------------------------------------------------
00169   //
00170   // dense B <- alpha * A^T * (*this)
00171   //
00172   //-------------------------------------------------------------
00173   
00174   void EpetraMultiVec::MvTransMv ( const double alpha, const MultiVec<double>& A,
00175            Teuchos::SerialDenseMatrix<int,double>& B
00176 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00177            , ConjType conj
00178 #endif
00179            ) const
00180   {    
00181     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00182     
00183     if (A_vec) {
00184       Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00185       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00186       
00187       int ret = B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 );
00188       assert( ret == 0 ); 
00189     }
00190   }
00191   
00192   //-------------------------------------------------------------
00193   //
00194   // b[i] = A[i]^T * this[i]
00195   // 
00196   //-------------------------------------------------------------
00197   
00198   void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double>* b
00199 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00200              , ConjType conj
00201 #endif
00202              ) const
00203   {
00204     EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00205     assert(A_vec!=NULL);
00206     if ((A_vec!=NULL) && (b!=NULL) && ( (int)b->size() >= A_vec->NumVectors() ) ) {
00207       int ret = this->Dot( *A_vec, &(*b)[0] );
00208       assert( ret == 0 );
00209     }
00210   }
00211 
00212   //-------------------------------------------------------------
00213   //
00214   // this[i] = alpha[i] * this[i]
00215   // 
00216   //-------------------------------------------------------------
00217   void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
00218   {
00219     // Check to make sure the vector is as long as the multivector has columns.
00220     int numvecs = this->NumVectors();
00221     assert( (int)alpha.size() == numvecs );
00222     
00223     int ret = 0;
00224     std::vector<int> tmp_index( 1, 0 );
00225     for (int i=0; i<numvecs; i++) {
00226       Epetra_MultiVector temp_vec(View, *this, &tmp_index[0], 1);
00227       ret = temp_vec.Scale( alpha[i] );
00228       assert (ret == 0);
00229       tmp_index[0]++;
00230     }
00231   }
00232 
00234   //
00235   //--------Anasazi::EpetraOp Implementation-------------------------------------
00236   //
00238 
00239   //
00240   // AnasaziOperator constructors
00241   //  
00242   EpetraOp::EpetraOp(const Teuchos::RefCountPtr<Epetra_Operator> &Op) 
00243     : Epetra_Op(Op)
00244   {
00245   }
00246   
00247   EpetraOp::~EpetraOp() 
00248   {
00249   }
00250   //
00251   // AnasaziOperator applications
00252   //
00253   void EpetraOp::Apply ( const MultiVec<double>& X, 
00254              MultiVec<double>& Y ) const 
00255   {
00256     //
00257     // This standard operator computes Y = A*X
00258     //
00259     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00260     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00261     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00262     
00263     assert( vec_X!=NULL && vec_Y!=NULL );
00264 
00265     int info = Epetra_Op->Apply( *vec_X, *vec_Y );
00266     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00267                         "Anasazi::EpetraOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00268   }
00269 
00271   //
00272   //--------Anasazi::EpetraGenOp Implementation----------------------------------
00273   //
00275 
00276   //
00277   // AnasaziOperator constructors
00278   //
00279   
00280   EpetraGenOp::EpetraGenOp(const Teuchos::RefCountPtr<Epetra_Operator> &AOp,
00281          const Teuchos::RefCountPtr<Epetra_Operator> &MOp,
00282          bool isAInverse_) 
00283     : isAInverse( isAInverse_ ), Epetra_AOp(AOp), Epetra_MOp(MOp) 
00284   {
00285   }
00286     
00287   EpetraGenOp::~EpetraGenOp() 
00288   {
00289   }
00290   //
00291   // AnasaziOperator applications
00292   //
00293   void EpetraGenOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00294   {
00295     //
00296     // This generalized operator computes Y = A^{-1}*M*X
00297     //
00298     int info=0;
00299     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00300     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00301     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00302     Epetra_MultiVector temp_Y(*vec_Y); 
00303     
00304     assert( vec_X!=NULL && vec_Y!=NULL );
00305     //
00306     // Need to cast away constness because the member function Apply is not declared const.  
00307     // Change the transpose setting for the operator if necessary and change it back when done.
00308     //
00309     // Apply M
00310     info = Epetra_MOp->Apply( *vec_X, temp_Y );
00311     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00312                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00313     // Apply A or A^{-1}
00314     if (isAInverse) {
00315       info = Epetra_AOp->ApplyInverse( temp_Y, *vec_Y );
00316     }
00317     else {
00318       info = Epetra_AOp->Apply( temp_Y, *vec_Y );
00319     }
00320     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00321                         "Anasazi::EpetraGenOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00322   }
00323   
00324   int EpetraGenOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00325   {
00326     //
00327     // This generalized operator computes Y = A^{-1}*M*X 
00328     //
00329     int info=0;
00330     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00331     
00332     // Apply M
00333     info = Epetra_MOp->Apply( X, temp_Y );
00334     if (info!=0) return info;
00335     
00336     // Apply A or A^{-1}
00337     if (isAInverse)
00338       info = Epetra_AOp->ApplyInverse( temp_Y, Y );
00339     else
00340       info = Epetra_AOp->Apply( temp_Y, Y );
00341 
00342     return info;
00343   }
00344   
00345   int EpetraGenOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00346   {
00347     //
00348     // This generalized operator computes Y = M^{-1}*A*X 
00349     //
00350     int info=0;
00351     Epetra_MultiVector temp_Y(OperatorDomainMap(), Y.NumVectors()); 
00352     
00353     // Apply A or A^{-1}
00354     if (isAInverse)
00355       info = Epetra_AOp->Apply( X, temp_Y );
00356     else
00357       info = Epetra_AOp->ApplyInverse( X, temp_Y );
00358     if (info!=0) return info;
00359     
00360     // Apply M^{-1}
00361     info = Epetra_MOp->ApplyInverse( temp_Y, Y );
00362     
00363     return info;
00364   }
00365 
00367   //
00368   //--------Anasazi::EpetraSymOp Implementation----------------------------------
00369   //
00371 
00372   //
00373   // AnasaziOperator constructors
00374   //
00375   EpetraSymOp::EpetraSymOp(const Teuchos::RefCountPtr<Epetra_Operator> &Op, 
00376          const bool isTrans) 
00377     : Epetra_Op(Op), isTrans_(isTrans)
00378   {
00379   }
00380   
00381   EpetraSymOp::~EpetraSymOp() 
00382   {
00383   }
00384   //
00385   // AnasaziOperator applications
00386   //
00387   void EpetraSymOp::Apply ( const MultiVec<double>& X, 
00388           MultiVec<double>& Y ) const 
00389   {
00390     int info=0;
00391     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00392     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00393     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00394     Epetra_MultiVector* temp_vec = new Epetra_MultiVector( 
00395                 (isTrans_) ? Epetra_Op->OperatorDomainMap() 
00396                 : Epetra_Op->OperatorRangeMap(), 
00397                 vec_X->NumVectors() );
00398     
00399     assert( vec_X!=NULL && vec_Y!=NULL && temp_vec!=NULL );
00400     //
00401     // Need to cast away constness because the member function Apply
00402     // is not declared const.
00403     //
00404     // Transpose the operator (if isTrans_ = true)
00405     if (isTrans_) {
00406       info = Epetra_Op->SetUseTranspose( isTrans_ );
00407       if (info != 0) {
00408         delete temp_vec;
00409         TEST_FOR_EXCEPTION( true, OperatorError, 
00410                             "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00411       }
00412     }
00413     //
00414     // Compute A*X or A'*X 
00415     //
00416     info=Epetra_Op->Apply( *vec_X, *temp_vec );
00417     if (info!=0) { 
00418       delete temp_vec; 
00419       TEST_FOR_EXCEPTION( true, OperatorError, 
00420                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00421     }
00422     //
00423     // Transpose/Un-transpose the operator based on value of isTrans_
00424     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00425     if (info!=0) { 
00426       delete temp_vec; 
00427       TEST_FOR_EXCEPTION( true, OperatorError, 
00428                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00429     }
00430     
00431     // Compute A^T*(A*X) or A*A^T
00432     info=Epetra_Op->Apply( *temp_vec, *vec_Y );
00433     if (info!=0) { 
00434       delete temp_vec; 
00435       TEST_FOR_EXCEPTION( true, OperatorError, 
00436                           "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00437     }
00438     
00439     // Un-transpose the operator
00440     info=Epetra_Op->SetUseTranspose( false );
00441     delete temp_vec;
00442     TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00443                         "Anasazi::EpetraSymOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00444   }
00445   
00446   int EpetraSymOp::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00447   {
00448     int info=0;
00449     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00450     //
00451     // This generalized operator computes Y = A^T*A*X or Y = A*A^T*X
00452     //
00453     // Transpose the operator (if isTrans_ = true)
00454     if (isTrans_) {
00455       info=Epetra_Op->SetUseTranspose( isTrans_ );
00456       if (info!=0) { return info; }
00457     }
00458     //
00459     // Compute A*X or A^T*X 
00460     //
00461     info=Epetra_Op->Apply( X, temp_vec );
00462     if (info!=0) { return info; }
00463     //
00464     // Transpose/Un-transpose the operator based on value of isTrans_
00465     info=Epetra_Op->SetUseTranspose( !isTrans_ );
00466     if (info!=0) { return info; }
00467     
00468     // Compute A^T*(A*X) or A*A^T
00469     info=Epetra_Op->Apply( temp_vec, Y );
00470     if (info!=0) { return info; }
00471     
00472     // Un-transpose the operator
00473     info=Epetra_Op->SetUseTranspose( false );
00474     return info;
00475   }
00476   
00477   int EpetraSymOp::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
00478   {
00479     int info=0;
00480     Epetra_MultiVector temp_vec(OperatorDomainMap(), Y.NumVectors()); 
00481     //
00482     // This generalized operator computes Y = (A^T*A)^{-1}*X or Y = (A*A^T)^{-1}*X
00483     //
00484     // Transpose the operator (if isTrans_ = true)
00485     if (!isTrans_) {
00486       info=Epetra_Op->SetUseTranspose( !isTrans_ );
00487       if (info!=0) { return info; }
00488     }
00489     //
00490     // Compute A^{-1}*X or A^{-T}*X 
00491     //
00492     info=Epetra_Op->ApplyInverse( X, temp_vec );
00493     if (info!=0) { return info; }
00494     //
00495     // Transpose/Un-transpose the operator based on value of isTrans_
00496     info=Epetra_Op->SetUseTranspose( isTrans_ );
00497     if (info!=0) { return info; }
00498     
00499     // Compute A^{-T}*(A^{-1}*X) or A^{-1}*A^{-T}
00500     info=Epetra_Op->Apply( temp_vec, Y );
00501     if (info!=0) { return info; }
00502     
00503     // Un-transpose the operator
00504     info=Epetra_Op->SetUseTranspose( false );
00505     return info;
00506   }
00507 
00509   //
00510   //--------Anasazi::EpetraSymMVOp Implementation--------------------------------
00511   //
00513 
00514   //
00515   // Anasazi::Operator constructors
00516   //
00517   EpetraSymMVOp::EpetraSymMVOp(const Teuchos::RefCountPtr<Epetra_MultiVector> &MV, const bool isTrans) 
00518     : Epetra_MV(MV), isTrans_(isTrans)
00519   {
00520     if (isTrans)
00521       MV_localmap = Teuchos::rcp( new Epetra_LocalMap( Epetra_MV->NumVectors(), 0, Epetra_MV->Map().Comm() ) );
00522     else
00523       MV_blockmap = Teuchos::rcp( &Epetra_MV->Map(), false );
00524   }
00525   
00526   //
00527   // AnasaziOperator applications
00528   //
00529   void EpetraSymMVOp::Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const 
00530   {
00531     int info=0;
00532     MultiVec<double> & temp_X = const_cast<MultiVec<double> &>(X);
00533     Epetra_MultiVector* vec_X = dynamic_cast<Epetra_MultiVector* >(&temp_X);
00534     Epetra_MultiVector* vec_Y = dynamic_cast<Epetra_MultiVector* >(&Y);
00535     
00536     if (isTrans_) {
00537       
00538       Epetra_MultiVector temp_vec( *MV_localmap, temp_X.GetNumberVecs() );
00539       
00540       /* A'*X */
00541       info = temp_vec.Multiply( 'T', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00542       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00543                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00544       
00545       /* A*(A'*X) */
00546       info = vec_Y->Multiply( 'N', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );      
00547       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00548                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00549     } 
00550     else {
00551       
00552       Epetra_MultiVector temp_vec( *MV_blockmap, temp_X.GetNumberVecs() );
00553       
00554       /* A*X */
00555       info = temp_vec.Multiply( 'N', 'N', 1.0, *Epetra_MV, *vec_X, 0.0 );
00556       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00557                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00558       
00559       /* A'*(A*X) */
00560       info = vec_Y->Multiply( 'T', 'N', 1.0, *Epetra_MV, temp_vec, 0.0 );
00561       TEST_FOR_EXCEPTION( info != 0, OperatorError, 
00562                           "Anasazi::EpetraSymMVOp::Apply(): Error returned from Epetra_Operator::Apply()" );
00563     }
00564   }
00565 
00566 } // end namespace Anasazi

Generated on Thu Sep 18 12:31:36 2008 for Anasazi by doxygen 1.3.9.1