ExampleTridiagSpmdLinearOp.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
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 #ifndef THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
00030 #define THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
00031 
00032 #include "Thyra_SpmdLinearOpBase.hpp"
00033 #include "Teuchos_DefaultComm.hpp"
00034 
00132 template<class Scalar>
00133 class ExampleTridiagSpmdLinearOp : public Thyra::SpmdLinearOpBase<Scalar> {
00134 private:
00135 
00136   Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >  comm_;
00137   int                  procRank_;
00138   int                  numProcs_;
00139   Thyra::Index         localDim_;
00140   std::vector<Scalar>  lower_;   // size = ( procRank == 0         ? localDim - 1 : localDim )    
00141   std::vector<Scalar>  diag_;    // size = localDim
00142   std::vector<Scalar>  upper_;   // size = ( procRank == numProc-1 ? localDim - 1 : localDim )
00143 
00144   void communicate( const bool first, const bool last, const Scalar x[], Scalar *x_km1, Scalar *x_kp1 ) const;
00145 
00146 public:
00147 
00149   using Thyra::SpmdLinearOpBase<Scalar>::euclideanApply;
00150 
00152   ExampleTridiagSpmdLinearOp() : procRank_(0), numProcs_(0) {}
00153 
00155   ExampleTridiagSpmdLinearOp(
00156     const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> > &comm
00157     ,const Thyra::Index localDim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
00158     { this->initialize(comm,localDim,lower,diag,upper); }
00159   
00179   void initialize(
00180     const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> > &comm
00181     ,const Thyra::Index             localDim // >= 2
00182     ,const Scalar                   lower[]  // size == ( procRank == 0         ? localDim - 1 : localDim )
00183     ,const Scalar                   diag[]   // size == localDim
00184     ,const Scalar                   upper[]  // size == ( procRank == numProc-1 ? localDim - 1 : localDim )
00185     )
00186     {
00187       TEST_FOR_EXCEPT( localDim < 2 );
00188       this->setLocalDimensions(comm,localDim,localDim); // Needed to set up range() and domain()
00189       comm_  = Teuchos::DefaultComm<Thyra::Index>::getDefaultSerialComm(comm);
00190       localDim_ = localDim;
00191       numProcs_ = comm->getSize();
00192       procRank_ = comm->getRank();
00193       const Thyra::Index
00194         lowerDim = ( procRank_ == 0           ? localDim - 1 : localDim ),
00195         upperDim = ( procRank_ == numProcs_-1 ? localDim - 1 : localDim );
00196       lower_.resize(lowerDim);  for( int k = 0; k < lowerDim; ++k ) lower_[k] = lower[k];
00197       diag_.resize(localDim);   for( int k = 0; k < localDim; ++k ) diag_[k]  = diag[k];
00198       upper_.resize(upperDim);  for( int k = 0; k < upperDim; ++k ) upper_[k] = upper[k];
00199     }
00200 
00201   // Overridden form Teuchos::Describable */
00202 
00203   std::string description() const
00204     {
00205       return (std::string("ExampleTridiagSpmdLinearOp<") + Teuchos::ScalarTraits<Scalar>::name() + std::string(">"));
00206     }
00207 
00208 protected:
00209 
00210 
00211   // Overridden from SingleScalarEuclideanLinearOpBase
00212 
00213   bool opSupported( Thyra::ETransp M_trans ) const
00214     {
00215       typedef Teuchos::ScalarTraits<Scalar> ST;
00216       return (M_trans == Thyra::NOTRANS || (!ST::isComplex && M_trans == Thyra::CONJ) );
00217     }
00218 
00219   // Overridden from SerialLinearOpBase
00220 
00221   void euclideanApply(
00222     const Thyra::ETransp                         M_trans
00223     ,const RTOpPack::ConstSubVectorView<Scalar>  &local_x_in
00224     ,const RTOpPack::SubVectorView<Scalar>       *local_y_out
00225     ,const Scalar                                alpha
00226     ,const Scalar                                beta
00227     ) const
00228     {
00229       typedef Teuchos::ScalarTraits<Scalar> ST;
00230       TEST_FOR_EXCEPTION( M_trans != Thyra::NOTRANS, std::logic_error, "Error, can not handle transpose!" );
00231       // Get constants
00232       const Scalar zero = ST::zero();
00233       // Get raw pointers to vector data to make me feel better!
00234       const Scalar *x = local_x_in.values();
00235       Scalar       *y = local_y_out->values();
00236       // Determine what process we are
00237       const bool first = ( procRank_ == 0 ), last = ( procRank_ == numProcs_-1 );
00238       // Communicate ghost elements
00239       Scalar x_km1, x_kp1;
00240       communicate( first, last, x, &x_km1, &x_kp1 );
00241       // Perform operation (if beta==0 then we must be careful since y could be uninitialized on input!)
00242       Thyra::Index k = 0, lk = 0;
00243       if( beta == zero ) {
00244         y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k] + upper_[k]*x[k+1] ); if(!first) ++lk;             // First local row
00245         for( k = 1; k < localDim_ - 1; ++lk, ++k )
00246           y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );                                        // Middle local rows
00247         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + (last?zero:upper_[k]*x_kp1) );                               // Last local row
00248       }
00249       else {
00250         y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k] + upper_[k]*x[k+1] ) + beta*y[k]; if(!first) ++lk; // First local row
00251         for( k = 1; k < localDim_ - 1; ++lk, ++k )
00252           y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] ) + beta*y[k];                            // Middle local rows
00253         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + (last?zero:upper_[k]*x_kp1) ) + beta*y[k];                   // Last local row
00254       }
00255       //std::cout << "\ny = ["; for(k=0;k<localDim_;++k) { std::cout << y[k]; if(k<localDim_-1) std::cout << ","; } std::cout << "]\n";
00256     }
00257 
00258 };  // end class ExampleTridiagSpmdLinearOp
00259 
00260 // private
00261 
00262 template<class Scalar>
00263 void ExampleTridiagSpmdLinearOp<Scalar>::communicate(
00264   const bool first, const bool last, const Scalar x[], Scalar *x_km1, Scalar *x_kp1
00265   ) const
00266 {
00267   const bool isEven = (procRank_ % 2 == 0);
00268   const bool isOdd = !isEven;
00269   // 1) Send x[localDim_-1] from each even process forward to the next odd
00270   // process where it is received in x_km1.
00271   if( isEven && procRank_+1 < numProcs_ ) send(*comm_,x[localDim_-1],procRank_+1);
00272   if( isOdd && procRank_-1 >= 0 )         receive(*comm_,procRank_-1,x_km1);
00273   // 2) Send x[0] from each odd process backward to the previous even
00274   // process where it is received in x_kp1.
00275   if( isOdd && procRank_-1 >= 0 )         send(*comm_,x[0],procRank_-1);
00276   if( isEven && procRank_+1 < numProcs_ ) receive(*comm_,procRank_+1,x_kp1);
00277   // 3) Send x[localDim_-1] from each odd process forward to the next even
00278   // process where it is received in x_km1.
00279   if( isOdd && procRank_+1 < numProcs_ )  send(*comm_,x[localDim_-1],procRank_+1);
00280   if( isEven && procRank_-1 >= 0 )        receive(*comm_,procRank_-1,x_km1);
00281   // 4) Send x[0] from each even process backward to the previous odd
00282   // process where it is received in x_kp1.
00283   if(isEven && procRank_-1 >= 0 )         send(*comm_,x[0],procRank_-1);
00284   if(isOdd && procRank_+1 < numProcs_ )   receive(*comm_,procRank_+1,x_kp1);
00285 }
00286 
00287 #endif  // THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP

Generated on Thu Sep 18 12:33:01 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1