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, 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, bool zeroOut)
00053   : Epetra_MultiVector(Map, numvecs, zeroOut) 
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 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)); assert(tmp_vec!=NULL);
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 ( const double alpha, const MultiVec<double>& A, 
00140                const Teuchos::SerialDenseMatrix<int,double>& B, const 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)); assert(A_vec!=NULL);
00146   
00147   int info = Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta );
00148   TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00149          "Belos::EpetraMultiVec::MvTimesMatAddMv call to Multiply() returned a nonzero value.");
00150 
00151 }
00152 
00153 //-------------------------------------------------------------
00154 //
00155 // *this <- alpha * A + beta * B
00156 //
00157 //-------------------------------------------------------------
00158   
00159 void EpetraMultiVec::MvAddMv ( const double alpha , const MultiVec<double>& A, 
00160              const double beta, const MultiVec<double>& B) 
00161 {
00162   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); assert(A_vec!=NULL);
00163   EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B)); assert(B_vec!=NULL);
00164   
00165   int info = Update( alpha, *A_vec, beta, *B_vec, 0.0 );
00166   TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00167          "Belos::EpetraMultiVec::MvAddMv call to Update() returned a nonzero value.");
00168 }
00169 
00170 //-------------------------------------------------------------
00171 //
00172 // this[i] = alpha[i] * this[i]
00173 //
00174 //-------------------------------------------------------------
00175 void EpetraMultiVec::MvScale ( const std::vector<double>& alpha )
00176 {
00177   // Check to make sure the vector is as long as the multivector has columns.
00178   int numvecs = this->NumVectors();
00179   assert( (int)alpha.size() == numvecs );
00180   int ret = 0;
00181   std::vector<int> tmp_index( 1, 0 );
00182   for (int i=0; i<numvecs; i++) {
00183     Epetra_MultiVector temp_vec(View, *this, &tmp_index[0], 1);
00184     ret = temp_vec.Scale( alpha[i] );
00185     assert (ret == 0);
00186     tmp_index[0]++;
00187   }
00188 }
00189 
00190 //-------------------------------------------------------------
00191 //
00192 // dense B <- alpha * A^T * (*this)
00193 //
00194 //-------------------------------------------------------------
00195 
00196 void EpetraMultiVec::MvTransMv ( const double alpha, const MultiVec<double>& A,
00197          Teuchos::SerialDenseMatrix<int,double>& B) const
00198 {    
00199   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00200   
00201   if (A_vec) {
00202     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00203     Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00204     
00205     int info = B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 );
00206     TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00207            "Belos::EpetraMultiVec::MvTransMv call to Multiply() returned a nonzero value.");
00208   }
00209 }
00210 
00211 //-------------------------------------------------------------
00212 //
00213 // b[i] = A[i]^T * this[i]
00214 // 
00215 //-------------------------------------------------------------
00216 
00217 void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double>& b ) const
00218 {
00219   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); assert(A_vec!=NULL);
00220   if (A_vec && ( (int)b.size() >= A_vec->NumVectors() ) ) {
00221      int info = this->Dot( *A_vec, &b[0] );
00222      TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00223       "Belos::EpetraMultiVec::MvDot call to Dot() returned a nonzero value.");   
00224   }
00225 }
00226 
00227 //-------------------------------------------------------------
00228 //
00229 // alpha[i] = norm of i-th column of (*this)
00230 //
00231 //-------------------------------------------------------------
00232 
00233 void EpetraMultiVec::MvNorm ( std::vector<double>& normvec, NormType norm_type ) const {
00234   if ((int)normvec.size() >= GetNumberVecs()) {
00235     int info = 0;
00236     switch( norm_type ) {
00237     case ( OneNorm ) :
00238       info = Norm1(&normvec[0]);
00239       break;
00240     case ( TwoNorm ) :
00241       info = Norm2(&normvec[0]);
00242       break;
00243     case ( InfNorm ) :  
00244       info = NormInf(&normvec[0]);
00245       break;
00246     default:
00247       break;
00248     }
00249     TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00250            "Belos::EpetraMultiVec::MvNorm call to Norm() returned a nonzero value.");
00251   }
00252 }
00253 
00255 //
00256 // implementation of the BelosEpetraOp class.
00257 //
00259 //
00260 // BelosOperator constructors
00261 //
00262 
00263 EpetraOp::EpetraOp( const Teuchos::RCP<Epetra_Operator> &Op ) 
00264   : Epetra_Op(Op)
00265 {
00266 }
00267 
00268 //
00269 // BelosOperator applications
00270 //
00271 void EpetraOp::Apply ( const MultiVec<double>& x, 
00272                        MultiVec<double>& y, ETrans trans ) const {
00273   int info=0;
00274   MultiVec<double> & temp_x = const_cast<MultiVec<double> &>(x);
00275   Epetra_MultiVector* vec_x = dynamic_cast<Epetra_MultiVector* >(&temp_x);
00276   Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector* >(&y);
00277   //
00278   TEST_FOR_EXCEPTION( vec_x==NULL || vec_y==NULL, EpetraOpFailure, 
00279           "Belos::EpetraOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
00280   //
00281   // Set the operator to apply the transpose if that is requested.
00282   //
00283   if ( trans ) { 
00284     info = Epetra_Op->SetUseTranspose( true );
00285   }
00286   info = Epetra_Op->Apply( *vec_x, *vec_y );
00287   
00288   if ( trans ) { 
00289     info = Epetra_Op->SetUseTranspose( false );
00290   }
00291   
00292   TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00293          "Belos::EpetraOp::Apply call to Apply() returned a nonzero value."); 
00294 
00295 }
00296 
00298 //
00299 // implementation of the BelosEpetraPrecOp class.
00300 //
00302 //
00303 // BelosOperator constructors
00304 //
00305 
00306 EpetraPrecOp::EpetraPrecOp( const Teuchos::RCP<Epetra_Operator> &Op ) 
00307   : Epetra_Op(Op)
00308 {
00309 }
00310 
00311 //
00312 // BelosOperator applications
00313 //
00314 void EpetraPrecOp::Apply ( const MultiVec<double>& x, 
00315          MultiVec<double>& y, ETrans trans ) const {
00316   int info=0;
00317   MultiVec<double> & temp_x = const_cast<MultiVec<double> &>(x);
00318   Epetra_MultiVector* vec_x = dynamic_cast<Epetra_MultiVector* >(&temp_x);
00319   Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector* >(&y);
00320   //
00321   TEST_FOR_EXCEPTION( vec_x==NULL || vec_y==NULL, EpetraOpFailure, 
00322           "Belos::EpetraOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
00323   //
00324   // Set the operator to apply the transpose if that is requested.
00325   //
00326   if ( trans ) { 
00327     info = Epetra_Op->SetUseTranspose( true );
00328   }
00329   info = Epetra_Op->ApplyInverse( *vec_x, *vec_y );
00330   
00331   if ( trans ) { 
00332     info = Epetra_Op->SetUseTranspose( false );
00333   }
00334   
00335   TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00336          "Belos::EpetraOp::Apply call to Apply() returned a nonzero value."); 
00337 }
00338 
00339 int EpetraPrecOp::Apply ( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const 
00340 {
00341   //
00342   // This operator applies Y = A^{-1}*X
00343   //
00344   int info = Epetra_Op->ApplyInverse( X, Y );
00345   return info;
00346 }
00347 
00348 int EpetraPrecOp::ApplyInverse( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const
00349 {
00350   //
00351   // This operator applies Y = A*X
00352   //
00353   int info = Epetra_Op->Apply( X, Y );
00354   return info;
00355 }
00356 

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7