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