Thyra_DefaultClusteredSpmdProductVector.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_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP
00030 #define THYRA_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP
00031 
00032 #include "Thyra_DefaultClusteredSpmdProductVectorDecl.hpp"
00033 #include "Thyra_DefaultClusteredSpmdProductVectorSpace.hpp"
00034 #include "Thyra_SpmdVectorBase.hpp"
00035 #include "RTOp_parallel_helpers.h"
00036 #include "RTOpPack_SPMD_apply_op.hpp"
00037 #include "Teuchos_Workspace.hpp"
00038 #include "Teuchos_dyn_cast.hpp"
00039 #include "Teuchos_arrayArg.hpp"
00040 
00041 namespace Thyra {
00042 
00043 // Constructors/initializers/accessors
00044 
00045 template <class Scalar>
00046 DefaultClusteredSpmdProductVector<Scalar>::DefaultClusteredSpmdProductVector()
00047 {
00048   uninitialize();
00049 }
00050 
00051 template <class Scalar>
00052 DefaultClusteredSpmdProductVector<Scalar>::DefaultClusteredSpmdProductVector(
00053   const Teuchos::RefCountPtr<const DefaultClusteredSpmdProductVectorSpace<Scalar> >  &productSpace
00054   ,const Teuchos::RefCountPtr<VectorBase<Scalar> >                                   vecs[]
00055   )
00056 {
00057   initialize(productSpace,vecs);
00058 }
00059 
00060 template <class Scalar>
00061 void DefaultClusteredSpmdProductVector<Scalar>::initialize(
00062   const Teuchos::RefCountPtr<const DefaultClusteredSpmdProductVectorSpace<Scalar> >  &productSpace
00063   ,const Teuchos::RefCountPtr<VectorBase<Scalar> >                                   vecs[]
00064   )
00065 {
00066   // ToDo: Validate input!
00067   productSpace_ = productSpace;
00068   const int numBlocks = productSpace_->numBlocks();
00069   vecs_.resize(numBlocks);
00070   if(vecs) {
00071     std::copy( vecs, vecs + numBlocks, &vecs_[0] );
00072   }
00073   else {
00074     for( int k = 0; k < numBlocks; ++k )
00075       vecs_[k] = createMember(productSpace->getBlock(k));
00076   }
00077 }
00078 
00079 template <class Scalar>
00080 void DefaultClusteredSpmdProductVector<Scalar>::uninitialize(
00081   Teuchos::RefCountPtr<const DefaultClusteredSpmdProductVectorSpace<Scalar> >  *productSpace
00082   ,Teuchos::RefCountPtr<VectorBase<Scalar> >                                   vecs[]
00083   )
00084 {
00085   const int numBlocks = vecs_.size();
00086   if(productSpace) *productSpace = productSpace_;
00087   if(vecs) std::copy( &vecs_[0], &vecs_[0]+numBlocks, vecs );
00088   productSpace_ = Teuchos::null;
00089   vecs_.resize(0);
00090 }
00091 
00092 // Overridden from ProductVectorBase
00093 
00094 template <class Scalar>
00095 Teuchos::RefCountPtr<VectorBase<Scalar> >
00096 DefaultClusteredSpmdProductVector<Scalar>::getNonconstVectorBlock(const int k)
00097 {
00098   TEST_FOR_EXCEPT( ! ( 0 <= k && k < vecs_.size() ) );
00099   return vecs_[k];
00100 }
00101 
00102 template <class Scalar>
00103 Teuchos::RefCountPtr<const VectorBase<Scalar> >
00104 DefaultClusteredSpmdProductVector<Scalar>::getVectorBlock(const int k) const
00105 {
00106   TEST_FOR_EXCEPT( ! ( 0 <= k && k < vecs_.size() ) );
00107   return vecs_[k];
00108 }
00109 
00110 // Overridden from ProductVectorBase
00111 
00112 template <class Scalar>
00113 Teuchos::RefCountPtr<const ProductVectorSpaceBase<Scalar> >
00114 DefaultClusteredSpmdProductVector<Scalar>::productSpace() const
00115 {
00116   return productSpace_;
00117 }
00118 
00119 template <class Scalar>
00120 bool DefaultClusteredSpmdProductVector<Scalar>::blockIsConst(const int k) const
00121 {
00122   TEST_FOR_EXCEPT( ! ( 0 <= k && k < vecs_.size() ) );
00123   return false;
00124 }
00125 
00126 template <class Scalar>
00127 Teuchos::RefCountPtr<MultiVectorBase<Scalar> >
00128 DefaultClusteredSpmdProductVector<Scalar>::getNonconstMultiVectorBlock(const int k)
00129 {
00130   return getNonconstVectorBlock(k);
00131 }
00132 
00133 template <class Scalar>
00134 Teuchos::RefCountPtr<const MultiVectorBase<Scalar> >
00135 DefaultClusteredSpmdProductVector<Scalar>::getMultiVectorBlock(const int k) const
00136 {
00137   return getVectorBlock(k);
00138 }
00139 
00140 // Overridden from VectorBase
00141 
00142 template <class Scalar>
00143 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00144 DefaultClusteredSpmdProductVector<Scalar>::space() const
00145 {
00146   return productSpace_;
00147 }
00148 
00149 template <class Scalar>
00150 void DefaultClusteredSpmdProductVector<Scalar>::applyOp(
00151   const RTOpPack::RTOpT<Scalar>    &op
00152   ,const int                       num_vecs
00153   ,const VectorBase<Scalar>*const  vecs[]
00154   ,const int                       num_targ_vecs
00155   ,VectorBase<Scalar>*const        targ_vecs[]
00156   ,RTOpPack::ReductTarget          *reduct_obj
00157   ,const Index                     first_ele_offset_in
00158   ,const Index                     sub_dim_in
00159   ,const Index                     global_offset_in
00160   ) const
00161 {
00162   using Teuchos::dyn_cast;
00163   using Teuchos::Workspace;
00164   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00165   //
00166   const Index n = productSpace_->dim();
00167   // Validate input
00168 #ifdef TEUCHOS_DEBUG
00169   TEST_FOR_EXCEPTION(
00170     !(0 <= first_ele_offset_in && first_ele_offset_in < n), std::out_of_range
00171     ,"DefaultClusteredSpmdProductVector::applyOp(...): Error, "
00172     "first_ele_offset_in = " << first_ele_offset_in << " is not in range [0,"<<(n-1)<<"]" );
00173   TEST_FOR_EXCEPTION(
00174     global_offset_in < 0, std::invalid_argument
00175     ,"DefaultClusteredSpmdProductVector::applyOp(...): Error, "
00176     "global_offset_in = " << global_offset_in << " is not acceptable" );
00177   TEST_FOR_EXCEPTION(
00178     sub_dim_in > 0 && sub_dim_in - first_ele_offset_in > n, std::length_error
00179     ,"DefaultClusteredSpmdProductVector::applyOp(...): Error, "
00180     "global_offset_in = " << global_offset_in << ", sub_dim_in = " << sub_dim_in
00181     << "first_ele_offset_in = " << first_ele_offset_in << " and n = " << n
00182     << " are not compatible" );
00183   bool test_failed;
00184   for(int k = 0; k < num_vecs; ++k) {
00185     test_failed = !this->space()->isCompatible(*vecs[k]->space());
00186     TEST_FOR_EXCEPTION(
00187       test_failed, Exceptions::IncompatibleVectorSpaces
00188       ,"DefaultClusteredSpmdProductVector::applyOp(...): Error vecs["<<k<<"]->space() "
00189       <<"of type \'"<<typeid(*vecs[k]->space()).name()<<"\' is not compatible with this "
00190       <<"\'VectorSpaceBlocked\' vector space!"
00191       );
00192   }
00193   for(int k = 0; k < num_targ_vecs; ++k) {
00194     test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
00195     TEST_FOR_EXCEPTION(
00196       test_failed, Exceptions::IncompatibleVectorSpaces
00197       ,"DefaultClusteredSpmdProductVector::applyOp(...): Error targ_vecs["<<k<<"]->space() "
00198       <<"of type \'"<<typeid(*vecs[k]->space()).name()<<"\' is not compatible with this "
00199       <<"\'VectorSpaceBlocked\' vector space!"
00200       );
00201   }
00202 #endif
00203   //
00204   // Cast all of the vector arguments to DefaultClusteredSpmdProductVector and
00205   // make sure that they are all compatible!
00206   //
00207   Workspace<const DefaultClusteredSpmdProductVector<Scalar>*> cl_vecs(wss,num_vecs);
00208   for( int k = 0; k < num_vecs; ++k ) {
00209 #ifdef TEUCHOS_DEBUG
00210     TEST_FOR_EXCEPT(vecs[k]==NULL);
00211 #endif
00212     cl_vecs[k] = &dyn_cast<const DefaultClusteredSpmdProductVector<Scalar> >(*vecs[k]);
00213   }
00214   Workspace<DefaultClusteredSpmdProductVector<Scalar>*> cl_targ_vecs(wss,num_targ_vecs);
00215   for( int k = 0; k < num_targ_vecs; ++k ) {
00216 #ifdef TEUCHOS_DEBUG
00217     TEST_FOR_EXCEPT(targ_vecs[k]==NULL);
00218 #endif
00219     cl_targ_vecs[k] = &dyn_cast<DefaultClusteredSpmdProductVector<Scalar> >(*targ_vecs[k]);
00220   }
00221   //
00222   // Get the overlap of the element for this cluster that will participate in
00223   // the RTOp operation.
00224   //
00225   const Teuchos::RefCountPtr<const Teuchos::Comm<Index> >
00226     intraClusterComm = productSpace_->intraClusterComm(),
00227     interClusterComm = productSpace_->interClusterComm();
00228   const Index
00229     clusterSubDim = productSpace_->clusterSubDim(),
00230     clusterOffset = productSpace_->clusterOffset(),
00231     globalDim = productSpace_->dim();
00232   Index  overlap_first_cluster_ele_off  = 0;
00233   Index  overlap_cluster_sub_dim        = 0;
00234   Index  overlap_global_off             = 0;
00235   if(clusterSubDim) {
00236     RTOp_parallel_calc_overlap(
00237       globalDim,clusterSubDim,clusterOffset,first_ele_offset_in,sub_dim_in
00238       ,global_offset_in
00239       ,&overlap_first_cluster_ele_off,&overlap_cluster_sub_dim,&overlap_global_off
00240       );
00241   }
00242   //
00243   // Perform the RTOp for each set of block vectors just within this cluster
00244   // of processes.
00245   //
00246   Teuchos::RefCountPtr<RTOpPack::ReductTarget> i_reduct_obj;
00247   if(reduct_obj) i_reduct_obj = op.reduct_obj_create();
00248   // Note: i_reduct_obj will accumulate the reduction within this cluster of
00249   // processes.
00250   const int numBlocks = vecs_.size();
00251   if( overlap_first_cluster_ele_off >=0 ) {
00252     //
00253     // There is overlap with at least one element in one block
00254     // vector for this cluster
00255     //
00256     Workspace<const VectorBase<Scalar>*>  v_vecs(wss,num_vecs);
00257     Workspace<VectorBase<Scalar>*>        v_targ_vecs(wss,num_targ_vecs);
00258     Index overall_global_offset = overlap_global_off;
00259     for( int j = 0; j < numBlocks; ++j ) {
00260       const VectorBase<Scalar>
00261         &v = *vecs_[j];
00262       const VectorSpaceBase<Scalar>
00263         &v_space = *v.space();
00264       // Load up the constutuent block SPMD vectors
00265       for( int k = 0; k < num_vecs ; ++k )
00266         v_vecs[k] = &*cl_vecs[k]->vecs_[j];
00267       for( int k = 0; k < num_targ_vecs ; ++k )
00268         v_targ_vecs[k] = &*cl_targ_vecs[k]->vecs_[j];
00269       TEST_FOR_EXCEPTION(
00270         numBlocks > 1, std::logic_error
00271         ,"Error, Have not implemented general support for numBlocks > 1!"
00272         ); // ToDo: Fix the below code for numBlocks_ > 1!
00273       Index
00274         v_global_offset     = overall_global_offset,
00275         v_first_ele_offset  = overlap_first_cluster_ele_off,
00276         v_sub_dim           = overlap_cluster_sub_dim;
00277       // Apply RTOp on just this cluster
00278       Thyra::applyOp(
00279         op
00280         ,num_vecs, num_vecs ? &v_vecs[0] : NULL
00281         ,num_targ_vecs, num_targ_vecs ? &v_targ_vecs[0] : NULL
00282         ,i_reduct_obj.get() ? &*i_reduct_obj : NULL
00283         ,v_first_ele_offset,v_sub_dim,v_global_offset
00284         );
00285       //
00286       overall_global_offset += v_space.dim();
00287     }
00288   }
00289   //
00290   // Perform the global reduction across all of the root processes in each of
00291   // the clusters and then move the global reduction out to each of the
00292   // processes within the cluster.
00293   //
00294   if( reduct_obj ) {
00295     Teuchos::RefCountPtr<RTOpPack::ReductTarget>
00296       icl_reduct_obj = op.reduct_obj_create();
00297     // First, accumulate the global reduction across all of the elements by
00298     // just performing the global reduction involving the root processes of
00299     // each cluster.
00300     if(interClusterComm.get()) {
00301       RTOpPack::SPMD_all_reduce(
00302         &*interClusterComm
00303         ,op
00304         ,1
00305         ,Teuchos::arrayArg<const RTOpPack::ReductTarget*>(&*i_reduct_obj)()
00306         ,Teuchos::arrayArg<RTOpPack::ReductTarget*>(&*icl_reduct_obj)()
00307         );
00308     }
00309     // Now the root processes in each cluster have the full global reduction
00310     // across all elements stored in *icl_reduct_obj and the other processes
00311     // in each cluster have empty reductions in *icl_reduct_obj.  The last
00312     // thing to do is to just perform the reduction within each cluster of
00313     // processes and to add into the in/out *reduct_obj.
00314     RTOpPack::SPMD_all_reduce(
00315       &*intraClusterComm
00316       ,op
00317       ,1
00318       ,Teuchos::arrayArg<const RTOpPack::ReductTarget*>(&*icl_reduct_obj)()
00319       ,Teuchos::arrayArg<RTOpPack::ReductTarget*>(reduct_obj)()
00320       );
00321     // ToDo: Replace the above operation with a reduction across clustes into
00322     // reduct_obj in the root processes and then broadcast the result to all
00323     // of the slave processes.
00324   }
00325 }
00326 
00327 } // namespace Thyra
00328 
00329 #endif // THYRA_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP

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