Anasazi Version of the Day
Tsqr_RMessenger.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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 __TSQR_RMessenger_hpp
00030 #define __TSQR_RMessenger_hpp
00031 
00032 #include <Tsqr_MatView.hpp>
00033 #include <Tsqr_MessengerBase.hpp>
00034 
00035 #include <Teuchos_RCP.hpp>
00036 
00037 #include <algorithm>
00038 #include <vector>
00039 
00042 
00043 namespace TSQR {
00049   template< class Ordinal, class Scalar >
00050   class RMessenger {
00051   public:
00052     typedef Scalar scalar_type;
00053     typedef Ordinal ordinal_type;
00054     typedef MessengerBase< Scalar > messenger_type;
00055     typedef Teuchos::RCP< messenger_type > messenger_ptr;
00056 
00058     RMessenger (const messenger_ptr& messenger) :
00059       messenger_ (messenger) {}
00060 
00061     template< class ConstMatrixViewType >
00062     void
00063     send (const ConstMatrixViewType& R, const int destProc)
00064     {
00065       pack (R);
00066       messenger_->send (&buffer_[0], buffer_.size(), destProc, 0);
00067     }
00068 
00069     template< class MatrixViewType >
00070     void
00071     recv (MatrixViewType& R, const int srcProc)
00072     {
00073       const typename MatrixViewType::ordinal_type ncols = R.ncols();
00074       const Ordinal buflen = buffer_length (ncols);
00075       buffer_.resize (buflen);
00076       messenger_->recv (&buffer_[0], buflen, srcProc, 0);
00077       unpack (R);
00078     }
00079 
00080     template< class MatrixViewType >
00081     void
00082     broadcast (MatrixViewType& R, const int rootProc)
00083     {
00084       const int myRank = messenger_->rank();
00085       if (myRank == rootProc)
00086   pack (R);
00087       messenger_->broadcast (&buffer_[0], buffer_length (R.ncols()), rootProc);
00088       if (myRank != rootProc)
00089   unpack (R);
00090     }
00091 
00094     RMessenger (const RMessenger& rhs) :
00095       messenger_ (rhs.messenger_), buffer_ (0) // don't need to copy the buffer
00096     {}
00097 
00100     RMessenger& operator= (const RMessenger& rhs) {
00101       if (this != &rhs)
00102   {
00103     this->messenger_ = rhs.messenger_;
00104     // Don't need to do anything to this->buffer_; the various
00105     // operations such as pack() will resize it as necessary.
00106   }
00107       return *this;
00108     }
00109 
00110 
00111   private:
00112     messenger_ptr messenger_;
00113     std::vector< Scalar > buffer_;
00114 
00115     // Default construction doesn't make sense, so we forbid it.
00116     RMessenger ();
00117 
00122     Ordinal buffer_length (const Ordinal ncols) const {
00123       return (ncols * (ncols + Ordinal(1))) / Ordinal(2);
00124     }
00125 
00126     template< class ConstMatrixViewType >
00127     void
00128     pack (const ConstMatrixViewType& R)
00129     {
00130       typedef typename ConstMatrixViewType::scalar_type view_scalar_type;
00131       typedef typename ConstMatrixViewType::ordinal_type view_ordinal_type;
00132       typedef typename std::vector< Scalar >::iterator iter_type;
00133 
00134       const view_ordinal_type ncols = R.ncols();
00135       const Ordinal buf_length = buffer_length (ncols);
00136       buffer_.resize (buf_length);
00137       iter_type iter = buffer_.begin();
00138       for (view_ordinal_type j = 0; j < ncols; ++j)
00139   {
00140     const view_scalar_type* const R_j = &R(0,j);
00141     std::copy (R_j, R_j + (j+1), iter);
00142     iter += (j+1);
00143   }
00144     }
00145 
00146     template< class MatrixViewType >
00147     void
00148     unpack (MatrixViewType& R)
00149     {
00150       typedef typename MatrixViewType::ordinal_type view_ordinal_type;
00151       typedef typename std::vector< Scalar >::const_iterator const_iter_type;
00152 
00153       const view_ordinal_type ncols = R.ncols();
00154       const_iter_type iter = buffer_.begin();
00155       for (view_ordinal_type j = 0; j < ncols; ++j)
00156   {
00157     std::copy (iter, iter + (j+1), &R(0,j));
00158     iter += (j+1);
00159   }
00160     }
00161   };
00162 
00163 
00175   template< class MatrixViewType, class ConstMatrixViewType >
00176   void
00177   scatterStack (const ConstMatrixViewType& R_stack, 
00178     MatrixViewType& R_local,
00179     const Teuchos::RCP< MessengerBase< typename MatrixViewType::scalar_type > >& messenger)
00180   {
00181     typedef typename MatrixViewType::ordinal_type ordinal_type;
00182     typedef typename MatrixViewType::scalar_type scalar_type;
00183     typedef ConstMatView< ordinal_type, scalar_type > const_view_type;
00184 
00185     const int nprocs = messenger->size();
00186     const int my_rank = messenger->rank();
00187 
00188     if (my_rank == 0)
00189       {
00190   const ordinal_type ncols = R_stack.ncols();
00191 
00192   // Copy data from top ncols x ncols block of R_stack into R_local.
00193   const_view_type R_stack_view_first (ncols, ncols, R_stack.get(), R_stack.lda());
00194   R_local.copy (R_stack_view_first);
00195 
00196   // Loop through all other processors, sending each the next
00197   // ncols x ncols block of R_stack.
00198   RMessenger< ordinal_type, scalar_type > sender (messenger);
00199   for (int destProc = 1; destProc < nprocs; ++destProc)
00200     {
00201       const scalar_type* const R_ptr = R_stack.get() + destProc*ncols;
00202       const_view_type R_stack_view_cur (ncols, ncols, R_ptr, R_stack.lda());
00203       sender.send (R_stack_view_cur, destProc);
00204     }
00205       }
00206     else
00207       {
00208   const int srcProc = 0;
00209   R_local.fill (scalar_type(0));
00210   RMessenger< ordinal_type, scalar_type > receiver (messenger);
00211   receiver.recv (R_local, srcProc);
00212       }
00213   }
00214 
00215 
00216 
00217 
00218   template< class MatrixViewType, class ConstMatrixViewType >
00219   void
00220   gatherStack (MatrixViewType& R_stack, 
00221          ConstMatrixViewType& R_local,
00222          const Teuchos::RCP< MessengerBase< typename MatrixViewType::scalar_type > >& messenger)
00223   {
00224     typedef typename MatrixViewType::ordinal_type ordinal_type;
00225     typedef typename MatrixViewType::scalar_type scalar_type;
00226     typedef MatView< ordinal_type, scalar_type > matrix_view_type;
00227 
00228     const int nprocs = messenger->size();
00229     const int my_rank = messenger->rank();
00230 
00231     if (my_rank == 0)
00232       {
00233   const ordinal_type ncols = R_stack.ncols();
00234 
00235   // Copy data from R_local into top ncols x ncols block of R_stack.
00236   matrix_view_type R_stack_view_first (ncols, ncols, R_stack.get(), R_stack.lda());
00237   R_stack_view_first.copy (R_local);
00238 
00239   // Loop through all other processors, fetching their matrix data.
00240   RMessenger< ordinal_type, scalar_type > receiver (messenger);
00241   for (int srcProc = 1; srcProc < nprocs; ++srcProc)
00242     {
00243       const scalar_type* const R_ptr = R_stack.get() + srcProc*ncols;
00244       matrix_view_type R_stack_view_cur (ncols, ncols, R_ptr, R_stack.lda());
00245       // Fill (the lower triangle) with zeros, since
00246       // RMessenger::recv() only writes to the upper triangle.
00247       R_stack_view_cur.fill (scalar_type (0));
00248       receiver.recv (R_stack_view_cur, srcProc);
00249     }
00250       }
00251     else
00252       {
00253   // We only read R_stack on Proc 0, not on this proc.
00254   // Send data from R_local to Proc 0.
00255   const int destProc = 0;
00256   RMessenger< ordinal_type, scalar_type > sender (messenger);
00257   sender.send (R_local, destProc);
00258       }
00259     messenger->barrier();
00260   }
00261 
00262 } // namespace TSQR
00263 
00264 #endif // __TSQR_RMessenger_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends