Thyra Version of the Day
ExampleTridiagSpmdLinearOp.hpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
00043 #define THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
00044 
00045 #include "Thyra_LinearOpDefaultBase.hpp"
00046 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00047 #include "Thyra_DetachedSpmdVectorView.hpp"
00048 #include "Teuchos_DefaultComm.hpp"
00049 #include "Teuchos_CommHelpers.hpp"
00050 
00051 
00142 template<class Scalar>
00143 class ExampleTridiagSpmdLinearOp : public Thyra::LinearOpDefaultBase<Scalar> {
00144 public:
00145 
00147   ExampleTridiagSpmdLinearOp() {}
00148 
00150   ExampleTridiagSpmdLinearOp(
00151     const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00152     const Thyra::Ordinal localDim,
00153     const Teuchos::ArrayView<const Scalar> &lower,
00154     const Teuchos::ArrayView<const Scalar> &diag,
00155     const Teuchos::ArrayView<const Scalar> &upper
00156     )
00157     { this->initialize(comm, localDim, lower, diag, upper); }
00158   
00180   void initialize(
00181     const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00182     const Thyra::Ordinal localDim,
00183     const Teuchos::ArrayView<const Scalar> &lower,
00184     const Teuchos::ArrayView<const Scalar> &diag,
00185     const Teuchos::ArrayView<const Scalar> &upper
00186     );
00187 
00188 protected:
00189 
00190   // Overridden from LinearOpBase
00191 
00193   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > range() const
00194     { return space_; }
00195 
00197   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > domain() const
00198     { return space_; }
00199 
00201   bool opSupportedImpl(Thyra::EOpTransp M_trans) const
00202     {
00203       return (M_trans == Thyra::NOTRANS);
00204     }
00205 
00207   void applyImpl(
00208     const Thyra::EOpTransp M_trans,
00209     const Thyra::MultiVectorBase<Scalar> &X_in,
00210     const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00211     const Scalar alpha,
00212     const Scalar beta
00213     ) const;
00214 
00215 private:
00216 
00217   Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > comm_;
00218   Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > space_;
00219   Teuchos::Array<Scalar> lower_; // size = ( procRank == 0 ? localDim - 1 : localDim ) 
00220   Teuchos::Array<Scalar> diag_; // size = localDim
00221   Teuchos::Array<Scalar> upper_; // size = ( procRank == numProc-1 ? localDim - 1 : localDim )
00222 
00223   void communicate( const bool first, const bool last,
00224     const Teuchos::ArrayView<const Scalar> &x,
00225     const Teuchos::Ptr<Scalar> &x_km1,
00226     const Teuchos::Ptr<Scalar> &x_kp1 ) const;
00227 
00228 };  // end class ExampleTridiagSpmdLinearOp
00229 
00230 
00231 template<class Scalar>
00232 void ExampleTridiagSpmdLinearOp<Scalar>::initialize(
00233   const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm,
00234   const Thyra::Ordinal localDim,
00235   const Teuchos::ArrayView<const Scalar> &lower,
00236   const Teuchos::ArrayView<const Scalar> &diag,
00237   const Teuchos::ArrayView<const Scalar> &upper
00238   )
00239 {
00240   TEUCHOS_TEST_FOR_EXCEPT( localDim < 2 );
00241   comm_ = comm;
00242   space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm, localDim, -1);
00243   lower_ = lower;
00244   diag_ = diag;
00245   upper_ = upper;
00246 }
00247 
00248 
00249 template<class Scalar>
00250 void ExampleTridiagSpmdLinearOp<Scalar>::applyImpl(
00251   const Thyra::EOpTransp M_trans,
00252   const Thyra::MultiVectorBase<Scalar> &X_in,
00253   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00254   const Scalar alpha,
00255   const Scalar beta
00256   ) const
00257 {
00258 
00259   typedef Teuchos::ScalarTraits<Scalar> ST;
00260   typedef Thyra::Ordinal Ordinal;
00261   using Teuchos::outArg;
00262 
00263   TEUCHOS_ASSERT(this->opSupported(M_trans));
00264 
00265   // Get constants
00266   const Scalar zero = ST::zero();
00267   const Ordinal localDim = space_->localSubDim();
00268   const Ordinal procRank = comm_->getRank();
00269   const Ordinal numProcs = comm_->getSize();
00270   const Ordinal m = X_in.domain()->dim();
00271       
00272   // Loop over the input columns
00273 
00274   for (Ordinal col_j = 0; col_j < m; ++col_j) {
00275     
00276     // Get access the the elements of column j
00277     Thyra::ConstDetachedSpmdVectorView<Scalar> x_vec(X_in.col(col_j));
00278     Thyra::DetachedSpmdVectorView<Scalar> y_vec(Y_inout->col(col_j));
00279     const Teuchos::ArrayRCP<const Scalar> x = x_vec.sv().values();
00280     const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
00281 
00282     // Determine what process we are
00283     const bool first = (procRank == 0), last = ( procRank == numProcs-1 );
00284     
00285     // Communicate ghost elements
00286     Scalar x_km1, x_kp1;
00287     communicate( first, last, x(), outArg(x_km1), outArg(x_kp1) );
00288 
00289     // Perform operation (if beta==0 then we must be careful since y could
00290     // be uninitialized on input!)
00291     Thyra::Ordinal k = 0, lk = 0;
00292     if( beta == zero ) {
00293       // First local row
00294       y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
00295         + upper_[k]*x[k+1] );
00296       if(!first) ++lk;
00297       // Middle local rows
00298       for( k = 1; k < localDim - 1; ++lk, ++k )
00299         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
00300       // Last local row
00301       y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
00302         + (last?zero:upper_[k]*x_kp1) );
00303     }
00304     else {
00305       // First local row
00306       y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
00307         + upper_[k]*x[k+1] ) + beta*y[k];
00308       if(!first) ++lk;
00309       // Middle local rows
00310       for( k = 1; k < localDim - 1; ++lk, ++k )
00311         y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] )
00312           + beta*y[k];
00313       // Last local row
00314       y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
00315         + (last?zero:upper_[k]*x_kp1) ) + beta*y[k];
00316     }
00317     
00318   }
00319 
00320 }
00321 
00322 
00323 template<class Scalar>
00324 void ExampleTridiagSpmdLinearOp<Scalar>::communicate(
00325   const bool first, const bool last,
00326   const Teuchos::ArrayView<const Scalar> &x,
00327     const Teuchos::Ptr<Scalar> &x_km1,
00328     const Teuchos::Ptr<Scalar> &x_kp1
00329   ) const
00330 {
00331   typedef Thyra::Ordinal Ordinal;
00332   const Ordinal localDim = space_->localSubDim();
00333   const Ordinal procRank = comm_->getRank();
00334   const Ordinal numProcs = comm_->getSize();
00335   const bool isEven = (procRank % 2 == 0);
00336   const bool isOdd = !isEven;
00337   // 1) Send x[localDim-1] from each even process forward to the next odd
00338   // process where it is received in x_km1.
00339   if(isEven && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
00340   if(isOdd && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
00341   // 2) Send x[0] from each odd process backward to the previous even
00342   // process where it is received in x_kp1.
00343   if( isOdd && procRank-1 >= 0 )         send(*comm_, x[0], procRank-1);
00344   if( isEven && procRank+1 < numProcs ) receive(*comm_, procRank+1, &*x_kp1);
00345   // 3) Send x[localDim-1] from each odd process forward to the next even
00346   // process where it is received in x_km1.
00347   if (isOdd && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
00348   if (isEven && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
00349   // 4) Send x[0] from each even process backward to the previous odd
00350   // process where it is received in x_kp1.
00351   if (isEven && procRank-1 >= 0) send(*comm_, x[0], procRank-1);
00352   if (isOdd && procRank+1 < numProcs) receive(*comm_, procRank+1, &*x_kp1);
00353 }
00354 
00355 
00356 #endif  // THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines