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)
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 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);
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. 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   assert( Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta ) == 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)); assert(A_vec!=NULL);
00160   EpetraMultiVec *B_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(B)); assert(B_vec!=NULL);
00161   
00162   assert ( Update( alpha, *A_vec, beta, *B_vec, 0.0 ) == 0 ); 
00163 }
00164 
00165 //-------------------------------------------------------------
00166 //
00167 // dense B <- alpha * A^T * (*this)
00168 //
00169 //-------------------------------------------------------------
00170 
00171 void EpetraMultiVec::MvTransMv ( const double alpha, const MultiVec<double>& A,
00172          Teuchos::SerialDenseMatrix<int,double>& B) const
00173 {    
00174   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
00175   
00176   if (A_vec) {
00177     Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
00178     Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00179     
00180     assert ( B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) == 0 ); 
00181   }
00182 }
00183 
00184 //-------------------------------------------------------------
00185 //
00186 // b[i] = A[i]^T * this[i]
00187 // 
00188 //-------------------------------------------------------------
00189 
00190 void EpetraMultiVec::MvDot ( const MultiVec<double>& A, std::vector<double>* b ) const
00191 {
00192   EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); assert(A_vec!=NULL);
00193   if (A_vec && b && ( (int)b->size() >= A_vec->NumVectors() ) ) {
00194     assert( this->Dot( *A_vec, &(*b)[0] ) == 0 );
00195   }
00196 }
00197 
00198 //-------------------------------------------------------------
00199 //
00200 // alpha[i] = norm of i-th column of (*this)
00201 //
00202 //-------------------------------------------------------------
00203 
00204 void EpetraMultiVec::MvNorm ( std::vector<double>* normvec, NormType norm_type ) const {
00205   if (normvec && ((int)normvec->size() >= GetNumberVecs())) {
00206     switch( norm_type ) {
00207     case ( OneNorm ) :
00208       assert( Norm1(&(*normvec)[0]) == 0 );
00209       break;
00210     case ( TwoNorm ) :
00211       assert( Norm2(&(*normvec)[0]) == 0 );
00212       break;
00213     case ( InfNorm ) :  
00214       assert( NormInf(&(*normvec)[0]) == 0 );
00215       break;
00216     default:
00217       break;
00218     }
00219   }
00220 }
00221 
00223 //
00224 // implementation of the BelosEpetraOp class.
00225 //
00227 //
00228 // BelosOperator constructors
00229 //
00230 
00231 EpetraOp::EpetraOp( const Teuchos::RefCountPtr<Epetra_Operator> &Op ) 
00232   : Epetra_Op(Op)
00233 {
00234 }
00235 
00236 //
00237 // BelosOperator applications
00238 //
00239 ReturnType EpetraOp::Apply ( const MultiVec<double>& x, 
00240            MultiVec<double>& y, ETrans trans ) const {
00241   int info=0;
00242   MultiVec<double> & temp_x = const_cast<MultiVec<double> &>(x);
00243   Epetra_MultiVector* vec_x = dynamic_cast<Epetra_MultiVector* >(&temp_x);
00244   Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector* >(&y);
00245   //
00246   assert( vec_x!=NULL && vec_y!=NULL );
00247   //
00248   // Set the operator to apply the transpose if that is requested.
00249   // (TO DO:  insert a throw if the application returns a nonzero )
00250   //
00251   if ( trans ) { 
00252     info = Epetra_Op->SetUseTranspose( true );
00253     if (info != 0) { return Undefined; }
00254   }
00255   info = Epetra_Op->Apply( *vec_x, *vec_y );
00256   
00257   if ( trans ) { 
00258     info = Epetra_Op->SetUseTranspose( false );
00259     if (info != 0) { return Undefined; }
00260   }
00261   
00262   if (info != 0) { return Error; }
00263   return Ok;  
00264 }
00265 
00267 //
00268 // implementation of the BelosEpetraPrecOp class.
00269 //
00271 //
00272 // BelosOperator constructors
00273 //
00274 
00275 EpetraPrecOp::EpetraPrecOp( const Teuchos::RefCountPtr<Epetra_Operator> &Op ) 
00276   : Epetra_Op(Op)
00277 {
00278 }
00279 
00280 //
00281 // BelosOperator applications
00282 //
00283 ReturnType EpetraPrecOp::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   assert( vec_x!=NULL && vec_y!=NULL );
00291   //
00292   // Set the operator to apply the transpose if that is requested.
00293   // (TO DO:  insert a throw if the application returns a nonzero )
00294   //
00295   if ( trans ) { 
00296     info = Epetra_Op->SetUseTranspose( true );
00297     if (info != 0) { return Undefined; }
00298   }
00299   info = Epetra_Op->ApplyInverse( *vec_x, *vec_y );
00300   
00301   if ( trans ) { 
00302     info = Epetra_Op->SetUseTranspose( false );
00303     if (info != 0) { return Undefined; }
00304   }
00305   
00306   if (info != 0) { return Error; }
00307   return Ok;  
00308 }
00309 

Generated on Thu Sep 18 12:30:12 2008 for Belos by doxygen 1.3.9.1