BelosEpetraAdapter.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers 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 
00034 #include "BelosEpetraAdapter.hpp"
00035 
00036 using namespace Belos;
00037   
00038   //-------------------------------------------------------------
00039   
00041   // Construction/Destruction
00043   
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, bool zeroOut)
00053   : Epetra_MultiVector(Map_in, numvecs, zeroOut) 
00054 {
00055 }
00056 
00057 
00058 EpetraMultiVec::EpetraMultiVec(Epetra_DataAccess CV_in, const Epetra_MultiVector& P_vec,        
00059              const std::vector<int>& index )
00060   : Epetra_MultiVector(CV_in, 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 EpetraMultiVec::~EpetraMultiVec() 
00072 {
00073 }
00074 //
00075 //  member functions inherited from Belos::MultiVec
00076 //
00077 //
00078 //  Simulating a virtual copy constructor. If we could rely on the co-variance
00079 //  of virtual functions, we could return a pointer to EpetraMultiVec
00080 //  (the derived type) instead of a pointer to the pure virtual base class.
00081 //
00082 
00083 MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
00084 {
00085   EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs, false);
00086   return ptr_apv; // safe upcast.
00087 }
00088 //
00089 //  the following is a virtual copy constructor returning
00090 //  a pointer to the pure virtual class. std::vector values are
00091 //  copied.
00092 //
00093 
00094 MultiVec<double>* EpetraMultiVec::CloneCopy() const
00095 {
00096   EpetraMultiVec *ptr_apv = new EpetraMultiVec(*this);
00097   return ptr_apv; // safe upcast
00098 }
00099 
00100 
00101 MultiVec<double>* EpetraMultiVec::CloneCopy ( const std::vector<int>& index ) const
00102 {
00103   EpetraMultiVec * ptr_apv = new EpetraMultiVec(Copy, *this, index);
00104   return ptr_apv; // safe upcast.
00105 }
00106 
00107 
00108 MultiVec<double>* EpetraMultiVec::CloneView ( const std::vector<int>& index ) 
00109 {
00110   EpetraMultiVec * ptr_apv = new EpetraMultiVec(View, *this, index);
00111   return ptr_apv; // safe upcast.
00112 }
00113   
00114 
00115 void EpetraMultiVec::SetBlock( const MultiVec<double>& A, const std::vector<int>& index ) 
00116 { 
00117   EpetraMultiVec temp_vec(View, *this, index);
00118   
00119   int numvecs = index.size();
00120   if ( A.GetNumberVecs() != numvecs ) {
00121     std::vector<int> index2( numvecs );
00122     for(int i=0; i<numvecs; i++)
00123       index2[i] = i;
00124     EpetraMultiVec *tmp_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00125     TEST_FOR_EXCEPTION(tmp_vec==NULL, EpetraMultiVecFailure,
00126                        "Belos::EpetraMultiVec::SetBlock cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
00127     EpetraMultiVec A_vec(View, *tmp_vec, index2);
00128     temp_vec.MvAddMv( 1.0, A_vec, 0.0, A_vec );
00129   }
00130   else {
00131     temp_vec.MvAddMv( 1.0, A, 0.0, A );
00132   }
00133 }               
00134 
00135 //-------------------------------------------------------------
00136 //
00137 // *this <- alpha * A * B + beta * (*this)
00138 //
00139 //-------------------------------------------------------------
00140 
00141 void EpetraMultiVec::MvTimesMatAddMv ( const double alpha, const MultiVec<double>& A, 
00142                const Teuchos::SerialDenseMatrix<int,double>& B, const double beta ) 
00143 {
00144   Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00145   Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00146   
00147   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00148   TEST_FOR_EXCEPTION(A_vec==NULL, EpetraMultiVecFailure,
00149                      "Belos::EpetraMultiVec::MvTimesMatAddMv cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
00150   
00151   int info = Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta );
00152   TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00153          "Belos::EpetraMultiVec::MvTimesMatAddMv call to Multiply() returned a nonzero value.");
00154 
00155 }
00156 
00157 //-------------------------------------------------------------
00158 //
00159 // *this <- alpha * A + beta * B
00160 //
00161 //-------------------------------------------------------------
00162   
00163 void EpetraMultiVec::MvAddMv ( const double alpha , const MultiVec<double>& A, 
00164              const double beta, const MultiVec<double>& B) 
00165 {
00166   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00167   TEST_FOR_EXCEPTION( A_vec==NULL, EpetraMultiVecFailure,
00168                      "Belos::EpetraMultiVec::MvAddMv cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
00169   EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B));
00170   TEST_FOR_EXCEPTION( B_vec==NULL, EpetraMultiVecFailure,
00171                      "Belos::EpetraMultiVec::MvAddMv cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
00172   
00173   int info = Update( alpha, *A_vec, beta, *B_vec, 0.0 );
00174   TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00175          "Belos::EpetraMultiVec::MvAddMv call to Update() returned a nonzero value.");
00176 }
00177 
00178 //-------------------------------------------------------------
00179 //
00180 // this[i] = alpha[i] * this[i]
00181 //
00182 //-------------------------------------------------------------
00183 void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
00184 {
00185   // Check to make sure the vector is as long as the multivector has columns.
00186   int numvecs = this->NumVectors();
00187   TEST_FOR_EXCEPTION((int)alpha.size() != numvecs, EpetraMultiVecFailure, 
00188        "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvScale scaling vector (alpha) not same size as number of input vectors (mv).");
00189   int ret = 0;
00190   std::vector<int> tmp_index( 1, 0 );
00191   for (int i=0; i<numvecs; i++) {
00192     Epetra_MultiVector temp_vec(View, *this, &tmp_index[0], 1);
00193     ret = temp_vec.Scale( alpha[i] );
00194     TEST_FOR_EXCEPTION(ret!=0, EpetraMultiVecFailure, 
00195                       "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvScale call to Scale() returned a nonzero value.");
00196     tmp_index[0]++;
00197   }
00198 }
00199 
00200 //-------------------------------------------------------------
00201 //
00202 // dense B <- alpha * A^T * (*this)
00203 //
00204 //-------------------------------------------------------------
00205 
00206 void EpetraMultiVec::MvTransMv ( const double alpha, const MultiVec<double>& A,
00207          Teuchos::SerialDenseMatrix<int,double>& B) const
00208 {    
00209   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00210   
00211   if (A_vec) {
00212     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00213     Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00214     
00215     int info = B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 );
00216     TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00217            "Belos::EpetraMultiVec::MvTransMv call to Multiply() returned a nonzero value.");
00218   }
00219 }
00220 
00221 //-------------------------------------------------------------
00222 //
00223 // b[i] = A[i]^T * this[i]
00224 // 
00225 //-------------------------------------------------------------
00226 
00227 void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double>& b ) const
00228 {
00229   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
00230   TEST_FOR_EXCEPTION(A_vec==NULL, EpetraMultiVecFailure,
00231                      "Belos::EpetraMultiVec::MvDot cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
00232   if (A_vec && ( (int)b.size() >= A_vec->NumVectors() ) ) {
00233      int info = this->Dot( *A_vec, &b[0] );
00234      TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00235       "Belos::EpetraMultiVec::MvDot call to Dot() returned a nonzero value.");   
00236   }
00237 }
00238 
00239 //-------------------------------------------------------------
00240 //
00241 // alpha[i] = norm of i-th column of (*this)
00242 //
00243 //-------------------------------------------------------------
00244 
00245 void EpetraMultiVec::MvNorm ( std::vector<double>& normvec, NormType norm_type ) const {
00246   if ((int)normvec.size() >= GetNumberVecs()) {
00247     int info = 0;
00248     switch( norm_type ) {
00249     case ( OneNorm ) :
00250       info = Norm1(&normvec[0]);
00251       break;
00252     case ( TwoNorm ) :
00253       info = Norm2(&normvec[0]);
00254       break;
00255     case ( InfNorm ) :  
00256       info = NormInf(&normvec[0]);
00257       break;
00258     default:
00259       break;
00260     }
00261     TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00262            "Belos::EpetraMultiVec::MvNorm call to Norm() returned a nonzero value.");
00263   }
00264 }
00265 
00267 //
00268 // implementation of the BelosEpetraOp class.
00269 //
00271 //
00272 // BelosOperator constructors
00273 //
00274 
00275 EpetraOp::EpetraOp( const Teuchos::RCP<Epetra_Operator> &Op ) 
00276   : Epetra_Op(Op)
00277 {
00278 }
00279 
00280 //
00281 // BelosOperator applications
00282 //
00283 void EpetraOp::Apply ( const MultiVec<double>& x, 
00284                        MultiVec<double>& y, ETrans trans ) const {
00285   int info=0;
00286   MultiVec<double> & temp_x = const_cast<MultiVec<double> &>(x);
00287   Epetra_MultiVector* vec_x = dynamic_cast<Epetra_MultiVector* >(&temp_x);
00288   Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector* >(&y);
00289   //
00290   TEST_FOR_EXCEPTION( vec_x==NULL || vec_y==NULL, EpetraOpFailure, 
00291           "Belos::EpetraOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
00292   //
00293   // Set the operator to apply the transpose if that is requested.
00294   //
00295   if ( trans ) { 
00296     info = Epetra_Op->SetUseTranspose( true );
00297   }
00298   info = Epetra_Op->Apply( *vec_x, *vec_y );
00299   
00300   if ( trans ) { 
00301     info = Epetra_Op->SetUseTranspose( false );
00302   }
00303   
00304   TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00305          "Belos::EpetraOp::Apply call to Apply() returned a nonzero value."); 
00306 
00307 }
00308 
00310 //
00311 // implementation of the BelosEpetraPrecOp class.
00312 //
00314 //
00315 // BelosOperator constructors
00316 //
00317 
00318 EpetraPrecOp::EpetraPrecOp( const Teuchos::RCP<Epetra_Operator> &Op ) 
00319   : Epetra_Op(Op)
00320 {
00321 }
00322 
00323 //
00324 // BelosOperator applications
00325 //
00326 void EpetraPrecOp::Apply ( const MultiVec<double>& x, 
00327          MultiVec<double>& y, ETrans trans ) const {
00328   int info=0;
00329   MultiVec<double> & temp_x = const_cast<MultiVec<double> &>(x);
00330   Epetra_MultiVector* vec_x = dynamic_cast<Epetra_MultiVector* >(&temp_x);
00331   Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector* >(&y);
00332   //
00333   TEST_FOR_EXCEPTION( vec_x==NULL || vec_y==NULL, EpetraOpFailure, 
00334           "Belos::EpetraOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
00335   //
00336   // Set the operator to apply the transpose if that is requested.
00337   //
00338   if ( trans ) { 
00339     info = Epetra_Op->SetUseTranspose( true );
00340   }
00341   info = Epetra_Op->ApplyInverse( *vec_x, *vec_y );
00342   
00343   if ( trans ) { 
00344     info = Epetra_Op->SetUseTranspose( false );
00345   }
00346   
00347   TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00348          "Belos::EpetraOp::Apply call to Apply() returned a nonzero value."); 
00349 }
00350 
00351 int EpetraPrecOp::Apply ( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const 
00352 {
00353   //
00354   // This operator applies Y = A^{-1}*X
00355   //
00356   int info = Epetra_Op->ApplyInverse( X, Y );
00357   return info;
00358 }
00359 
00360 int EpetraPrecOp::ApplyInverse( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const
00361 {
00362   //
00363   // This operator applies Y = A*X
00364   //
00365   int info = Epetra_Op->Apply( X, Y );
00366   return info;
00367 }
00368 

Generated on Wed May 12 21:45:50 2010 for Belos by  doxygen 1.4.7