Thyra Package Browser (Single Doxygen Collection) Version of the Day
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_LinearOpDefaultBase.hpp"
00033 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00034 #include "Thyra_DetachedSpmdVectorView.hpp"
00035 #include "Teuchos_DefaultComm.hpp"
00036 #include "Teuchos_CommHelpers.hpp"
00037 
00038 
00129 template<class Scalar>
00130 class ExampleTridiagSpmdLinearOp : public Thyra::LinearOpDefaultBase<Scalar> {
00131 public:
00132 
00134   ExampleTridiagSpmdLinearOp() {}
00135 
00137   ExampleTridiagSpmdLinearOp(
00138     const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00139     const Thyra::Ordinal localDim,
00140     const Teuchos::ArrayView<const Scalar> &lower,
00141     const Teuchos::ArrayView<const Scalar> &diag,
00142     const Teuchos::ArrayView<const Scalar> &upper
00143     )
00144     { this->initialize(comm, localDim, lower, diag, upper); }
00145   
00167   void initialize(
00168     const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00169     const Thyra::Ordinal localDim,
00170     const Teuchos::ArrayView<const Scalar> &lower,
00171     const Teuchos::ArrayView<const Scalar> &diag,
00172     const Teuchos::ArrayView<const Scalar> &upper
00173     );
00174 
00175 protected:
00176 
00177   // Overridden from LinearOpBase
00178 
00180   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > range() const
00181     { return space_; }
00182 
00184   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > domain() const
00185     { return space_; }
00186 
00188   bool opSupportedImpl(Thyra::EOpTransp M_trans) const
00189     {
00190       return (M_trans == Thyra::NOTRANS);
00191     }
00192 
00194   void applyImpl(
00195     const Thyra::EOpTransp M_trans,
00196     const Thyra::MultiVectorBase<Scalar> &X_in,
00197     const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00198     const Scalar alpha,
00199     const Scalar beta
00200     ) const;
00201 
00202 private:
00203 
00204   Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > comm_;
00205   Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > space_;
00206   Teuchos::Array<Scalar> lower_; // size = ( procRank == 0 ? localDim - 1 : localDim ) 
00207   Teuchos::Array<Scalar> diag_; // size = localDim
00208   Teuchos::Array<Scalar> upper_; // size = ( procRank == numProc-1 ? localDim - 1 : localDim )
00209 
00210   void communicate( const bool first, const bool last,
00211     const Teuchos::ArrayView<const Scalar> &x,
00212     const Teuchos::Ptr<Scalar> &x_km1,
00213     const Teuchos::Ptr<Scalar> &x_kp1 ) const;
00214 
00215 };  // end class ExampleTridiagSpmdLinearOp
00216 
00217 
00218 template<class Scalar>
00219 void ExampleTridiagSpmdLinearOp<Scalar>::initialize(
00220   const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00221   const Thyra::Ordinal localDim,
00222   const Teuchos::ArrayView<const Scalar> &lower,
00223   const Teuchos::ArrayView<const Scalar> &diag,
00224   const Teuchos::ArrayView<const Scalar> &upper
00225   )
00226 {
00227   TEST_FOR_EXCEPT( localDim < 2 );
00228   comm_ = comm;
00229   space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm, localDim, -1);
00230   lower_ = lower;
00231   diag_ = diag;
00232   upper_ = upper;
00233 }
00234 
00235 
00236 template<class Scalar>
00237 void ExampleTridiagSpmdLinearOp<Scalar>::applyImpl(
00238   const Thyra::EOpTransp M_trans,
00239   const Thyra::MultiVectorBase<Scalar> &X_in,
00240   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00241   const Scalar alpha,
00242   const Scalar beta
00243   ) const
00244 {
00245 
00246   typedef Teuchos::ScalarTraits<Scalar> ST;
00247   typedef Thyra::Ordinal Ordinal;
00248   using Teuchos::outArg;
00249 
00250   TEUCHOS_ASSERT(this->opSupported(M_trans));
00251 
00252   // Get constants
00253   const Scalar zero = ST::zero();
00254   const Ordinal localDim = space_->localSubDim();
00255   const Ordinal procRank = comm_->getRank();
00256   const Ordinal numProcs = comm_->getSize();
00257   const Ordinal m = X_in.domain()->dim();
00258       
00259   // Loop over the input columns
00260 
00261   for (Ordinal col_j = 0; col_j < m; ++col_j) {
00262     
00263     // Get access the the elements of column j
00264     Thyra::ConstDetachedSpmdVectorView<Scalar> x_vec(X_in.col(col_j));
00265     Thyra::DetachedSpmdVectorView<Scalar> y_vec(Y_inout->col(col_j));
00266     const Teuchos::ArrayRCP<const Scalar> x = x_vec.sv().values();
00267     const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
00268 
00269     // Determine what process we are
00270     const bool first = (procRank == 0), last = ( procRank == numProcs-1 );
00271     
00272     // Communicate ghost elements
00273     Scalar x_km1, x_kp1;
00274     communicate( first, last, x(), outArg(x_km1), outArg(x_kp1) );
00275 
00276     // Perform operation (if beta==0 then we must be careful since y could
00277     // be uninitialized on input!)
00278     Thyra::Ordinal k = 0, lk = 0;
00279     if( beta == zero ) {
00280       // First local row
00281       y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
00282         + upper_[k]*x[k+1] );
00283       if(!first) ++lk;
00284       // Middle local rows
00285       for( k = 1; k < localDim - 1; ++lk, ++k )
00286         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
00287       // Last local row
00288       y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
00289         + (last?zero:upper_[k]*x_kp1) );
00290     }
00291     else {
00292       // First local row
00293       y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
00294         + upper_[k]*x[k+1] ) + beta*y[k];
00295       if(!first) ++lk;
00296       // Middle local rows
00297       for( k = 1; k < localDim - 1; ++lk, ++k )
00298         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] )
00299           + beta*y[k];
00300       // Last local row
00301       y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
00302         + (last?zero:upper_[k]*x_kp1) ) + beta*y[k];
00303     }
00304     
00305   }
00306 
00307 }
00308 
00309 
00310 template<class Scalar>
00311 void ExampleTridiagSpmdLinearOp<Scalar>::communicate(
00312   const bool first, const bool last,
00313   const Teuchos::ArrayView<const Scalar> &x,
00314     const Teuchos::Ptr<Scalar> &x_km1,
00315     const Teuchos::Ptr<Scalar> &x_kp1
00316   ) const
00317 {
00318   typedef Thyra::Ordinal Ordinal;
00319   const Ordinal localDim = space_->localSubDim();
00320   const Ordinal procRank = comm_->getRank();
00321   const Ordinal numProcs = comm_->getSize();
00322   const bool isEven = (procRank % 2 == 0);
00323   const bool isOdd = !isEven;
00324   // 1) Send x[localDim-1] from each even process forward to the next odd
00325   // process where it is received in x_km1.
00326   if(isEven && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
00327   if(isOdd && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
00328   // 2) Send x[0] from each odd process backward to the previous even
00329   // process where it is received in x_kp1.
00330   if( isOdd && procRank-1 >= 0 )         send(*comm_, x[0], procRank-1);
00331   if( isEven && procRank+1 < numProcs ) receive(*comm_, procRank+1, &*x_kp1);
00332   // 3) Send x[localDim-1] from each odd process forward to the next even
00333   // process where it is received in x_km1.
00334   if (isOdd && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
00335   if (isEven && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
00336   // 4) Send x[0] from each even process backward to the previous odd
00337   // process where it is received in x_kp1.
00338   if (isEven && procRank-1 >= 0) send(*comm_, x[0], procRank-1);
00339   if (isOdd && procRank+1 < numProcs) receive(*comm_, procRank+1, &*x_kp1);
00340 }
00341 
00342 
00343 #endif  // THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines