Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultClusteredSpmdProductVector_def.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP
43#define THYRA_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP
44
45#include "Thyra_DefaultClusteredSpmdProductVector_decl.hpp"
46#include "Thyra_DefaultClusteredSpmdProductVectorSpace.hpp"
47#include "Thyra_SpmdVectorBase.hpp"
48#include "RTOp_parallel_helpers.h"
49#include "RTOpPack_SPMD_apply_op.hpp"
50#include "Teuchos_Workspace.hpp"
51#include "Teuchos_dyn_cast.hpp"
52#include "Teuchos_implicit_cast.hpp"
53
54
55namespace Thyra {
56
57
58// Constructors/initializers/accessors
59
60
61template <class Scalar>
66
67
68template <class Scalar>
76
77
78template <class Scalar>
81 ,const Teuchos::RCP<VectorBase<Scalar> > vecs[]
82 )
83{
84 // ToDo: Validate input!
85 productSpace_ = productSpace_in;
86 const int numBlocks = productSpace_->numBlocks();
87 vecs_.resize(numBlocks);
88 if(vecs) {
89 std::copy( vecs, vecs + numBlocks, &vecs_[0] );
90 }
91 else {
92 for( int k = 0; k < numBlocks; ++k )
93 vecs_[k] = createMember(productSpace_->getBlock(k));
94 }
95}
96
97
98template <class Scalar>
102 )
103{
104 const int numBlocks = vecs_.size();
105 if(productSpace_in) *productSpace_in = productSpace_;
106 if(vecs) std::copy( &vecs_[0], &vecs_[0]+numBlocks, vecs );
107 productSpace_ = Teuchos::null;
108 vecs_.resize(0);
109}
110
111
112// Overridden from ProductVectorBase
113
114
115template <class Scalar>
118{
120 TEUCHOS_TEST_FOR_EXCEPT( ! ( 0 <= k && k < implicit_cast<int>(vecs_.size()) ) );
121 return vecs_[k];
122}
123
124
125template <class Scalar>
128{
130 TEUCHOS_TEST_FOR_EXCEPT( ! ( 0 <= k && k < implicit_cast<int>(vecs_.size()) ) );
131 return vecs_[k];
132}
133
134
135// Overridden from ProductVectorBase
136
137
138template <class Scalar>
141{
142 return productSpace_;
143}
144
145
146template <class Scalar>
148{
150 TEUCHOS_TEST_FOR_EXCEPT( ! ( 0 <= k && k < implicit_cast<int>(vecs_.size()) ) );
151 return false;
152}
153
154
155template <class Scalar>
158{
159 return getNonconstVectorBlock(k);
160}
161
162
163template <class Scalar>
166{
167 return getVectorBlock(k);
168}
169
170
171// Overridden from VectorBase
172
173
174template <class Scalar>
177{
178 return productSpace_;
179}
180
181
182// Overridden protected members from VectorBase
183
184
185template <class Scalar>
187 const RTOpPack::RTOpT<Scalar> &op,
188 const ArrayView<const Ptr<const VectorBase<Scalar> > > &vecs,
189 const ArrayView<const Ptr<VectorBase<Scalar> > > &targ_vecs,
190 const Ptr<RTOpPack::ReductTarget> &reduct_obj,
191 const Ordinal global_offset_in
192 ) const
193{
194
195 const Ordinal first_ele_offset_in = 0;
196 const Ordinal sub_dim_in = -1;
197
198 using Teuchos::null;
199 using Teuchos::ptr_dynamic_cast;
200 using Teuchos::tuple;
201
202 const int num_vecs = vecs.size();
203 const int num_targ_vecs = targ_vecs.size();
204
205 // Validate input
206#ifdef TEUCHOS_DEBUG
208 global_offset_in < 0, std::invalid_argument,
209 "DefaultClusteredSpmdProductVector::applyOp(...): Error, "
210 "global_offset_in = " << global_offset_in << " is not acceptable" );
211 bool test_failed;
212 for (int k = 0; k < num_vecs; ++k) {
213 test_failed = !this->space()->isCompatible(*vecs[k]->space());
216 "DefaultClusteredSpmdProductVector::applyOp(...): Error vecs["<<k<<"]->space() "
217 <<"of type \'"<<typeName(*vecs[k]->space())<<"\' is not compatible with this "
218 <<"\'VectorSpaceBlocked\' vector space!"
219 );
220 }
221 for (int k = 0; k < num_targ_vecs; ++k) {
222 test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
225 ,"DefaultClusteredSpmdProductVector::applyOp(...): Error targ_vecs["<<k<<"]->space() "
226 <<"of type \'"<<typeName(*vecs[k]->space())<<"\' is not compatible with this "
227 <<"\'VectorSpaceBlocked\' vector space!"
228 );
229 }
230#endif
231
232 //
233 // Cast all of the vector arguments to DefaultClusteredSpmdProductVector and
234 // make sure that they are all compatible!
235 //
237 for ( int k = 0; k < num_vecs; ++k ) {
238#ifdef TEUCHOS_DEBUG
239 TEUCHOS_TEST_FOR_EXCEPT(vecs[k]==null);
240#endif
241 cl_vecs[k] = ptr_dynamic_cast<const DefaultClusteredSpmdProductVector<Scalar> >(vecs[k],true);
242 }
243 Array<Ptr<DefaultClusteredSpmdProductVector<Scalar> > > cl_targ_vecs(num_targ_vecs);
244 for ( int k = 0; k < num_targ_vecs; ++k ) {
245#ifdef TEUCHOS_DEBUG
246 TEUCHOS_TEST_FOR_EXCEPT(targ_vecs[k]==null);
247#endif
248 cl_targ_vecs[k] = ptr_dynamic_cast<DefaultClusteredSpmdProductVector<Scalar> >(targ_vecs[k],true);
249 }
250
251 //
252 // Get the overlap of the element for this cluster that will participate in
253 // the RTOp operation.
254 //
256 intraClusterComm = productSpace_->intraClusterComm(),
257 interClusterComm = productSpace_->interClusterComm();
258 const Ordinal
259 clusterSubDim = productSpace_->clusterSubDim(),
260 clusterOffset = productSpace_->clusterOffset(),
261 globalDim = productSpace_->dim();
262 Ordinal overlap_first_cluster_ele_off = 0;
263 Ordinal overlap_cluster_sub_dim = 0;
264 Ordinal overlap_global_off = 0;
265 if (clusterSubDim) {
266 RTOp_parallel_calc_overlap(
267 globalDim,clusterSubDim,clusterOffset,first_ele_offset_in,sub_dim_in
268 ,global_offset_in
269 ,&overlap_first_cluster_ele_off,&overlap_cluster_sub_dim,&overlap_global_off
270 );
271 }
272
273 //
274 // Perform the RTOp for each set of block vectors just within this cluster
275 // of processes.
276 //
278 if (!is_null(reduct_obj))
279 i_reduct_obj = op.reduct_obj_create();
280 // Note: i_reduct_obj will accumulate the reduction within this cluster of
281 // processes.
282 const int numBlocks = vecs_.size();
283 if (overlap_first_cluster_ele_off >=0) {
284
285 //
286 // There is overlap with at least one element in one block
287 // vector for this cluster
288 //
289 Array<Ptr<const VectorBase<Scalar> > > v_vecs(num_vecs);
290 Array<Ptr<VectorBase<Scalar> > > v_targ_vecs(num_targ_vecs);
291 Ordinal overall_global_offset = overlap_global_off;
292 for( int j = 0; j < numBlocks; ++j ) {
294 &v = *vecs_[j];
296 &v_space = *v.space();
297 // Load up the constutuent block SPMD vectors
298 for( int k = 0; k < num_vecs ; ++k )
299 v_vecs[k] = cl_vecs[k]->vecs_[j].ptr();
300 for( int k = 0; k < num_targ_vecs ; ++k )
301 v_targ_vecs[k] = cl_targ_vecs[k]->vecs_[j].ptr();
303 numBlocks > 1, std::logic_error
304 ,"Error, Have not implemented general support for numBlocks > 1!"
305 ); // ToDo: Fix the below code for numBlocks_ > 1!
306 Ordinal v_global_offset = overall_global_offset;
307 // Apply RTOp on just this cluster
308 Thyra::applyOp<Scalar>(
309 op, v_vecs(), v_targ_vecs(), i_reduct_obj.ptr(),
310 v_global_offset);
311 //
312 overall_global_offset += v_space.dim();
313 }
314
315 }
316
317 //
318 // Perform the global reduction across all of the root processes in each of
319 // the clusters and then move the global reduction out to each of the
320 // processes within the cluster.
321 //
322 if (!is_null(reduct_obj)) {
324 icl_reduct_obj = op.reduct_obj_create();
325 // First, accumulate the global reduction across all of the elements by
326 // just performing the global reduction involving the root processes of
327 // each cluster.
328 if (interClusterComm.get()) {
329 RTOpPack::SPMD_all_reduce(
330 &*interClusterComm,
331 op,
332 1,
333 tuple<const RTOpPack::ReductTarget*>(&*i_reduct_obj).getRawPtr(),
334 tuple<RTOpPack::ReductTarget*>(&*icl_reduct_obj).getRawPtr()
335 );
336 }
337 // Now the root processes in each cluster have the full global reduction
338 // across all elements stored in *icl_reduct_obj and the other processes
339 // in each cluster have empty reductions in *icl_reduct_obj. The last
340 // thing to do is to just perform the reduction within each cluster of
341 // processes and to add into the in/out *reduct_obj.
342 RTOpPack::SPMD_all_reduce(
343 &*intraClusterComm,
344 op,
345 1,
346 tuple<const RTOpPack::ReductTarget*>(&*icl_reduct_obj).getRawPtr(),
347 tuple<RTOpPack::ReductTarget*>(&*reduct_obj).getRawPtr()
348 );
349 // ToDo: Replace the above operation with a reduction across clustere into
350 // reduct_obj in the root processes and then broadcast the result to all
351 // of the slave processes.
352 }
353
354}
355
356
357} // namespace Thyra
358
359
360#endif // THYRA_DEFAULT_CLUSTERED_SPMD_PRODUCT_VECTOR_HPP
Teuchos::RCP< ReductTarget > reduct_obj_create() const
Ptr< T > ptr() const
void applyOpImpl(const RTOpPack::RTOpT< Scalar > &op, const ArrayView< const Ptr< const VectorBase< Scalar > > > &vecs, const ArrayView< const Ptr< VectorBase< Scalar > > > &targ_vecs, const Ptr< RTOpPack::ReductTarget > &reduct_obj, const Ordinal global_offset) const
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
void uninitialize(Teuchos::RCP< const DefaultClusteredSpmdProductVectorSpace< Scalar > > *productSpace=NULL, Teuchos::RCP< VectorBase< Scalar > > vecs[]=NULL)
Uninitialize.
Teuchos::RCP< const VectorSpaceBase< Scalar > > space() const
void initialize(const Teuchos::RCP< const DefaultClusteredSpmdProductVectorSpace< Scalar > > &productSpace, const Teuchos::RCP< VectorBase< Scalar > > vecs[])
Initialize.
Teuchos::RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
Teuchos::RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const
Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
Teuchos::RCP< VectorBase< Scalar > > getNonconstVectorBlock(const int k)
DefaultClusteredSpmdProductVector()
Constructs to uninitialized (see postconditions for uninitialize()).
Thrown if vector spaces are incompatible.
Abstract interface for finite-dimensional dense vectors.
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
Return a smart pointer to the vector space that this vector belongs to.
Abstract interface for objects that represent a space for vectors.
virtual Ordinal dim() const =0
Return the dimension of the vector space.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
TypeTo implicit_cast(const TypeFrom &t)
T_To & dyn_cast(T_From &from)