Thyra  Version of the Day
Thyra_MultiVectorDefaultBase_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_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
43 #define THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
44 
45 
46 #include "Thyra_MultiVectorDefaultBase_decl.hpp"
47 #include "Thyra_LinearOpDefaultBase.hpp"
48 #include "Thyra_MultiVectorStdOps.hpp"
49 #include "Thyra_VectorSpaceFactoryBase.hpp"
50 #include "Thyra_VectorSpaceBase.hpp"
51 #include "Thyra_VectorBase.hpp"
52 #include "Thyra_AssertOp.hpp"
53 #include "Thyra_DefaultColumnwiseMultiVector.hpp"
54 #include "RTOpPack_TOpAssignScalar.hpp"
55 #include "Teuchos_Workspace.hpp"
56 #include "Teuchos_Assert.hpp"
57 #include "Teuchos_as.hpp"
58 
59 
60 namespace Thyra {
61 
62 
63 // Overridden public member functions from MultiVectorBase
64 
65 
66 template<class Scalar>
67 RCP<MultiVectorBase<Scalar> >
69 {
71  &l_domain = *this->domain(),
72  &l_range = *this->range();
73  RCP<MultiVectorBase<Scalar> >
74  copy = createMembers(l_range,l_domain.dim());
75  ::Thyra::assign<Scalar>(copy.ptr(), *this);
76  return copy;
77 }
78 
79 
80 // protected
81 
82 
83 // Overridden protected member functions from MultiVectorBase
84 
85 template<class Scalar>
87 {
88  using Teuchos::tuple; using Teuchos::null;
89  RTOpPack::TOpAssignScalar<Scalar> assign_scalar_op(alpha);
90  Thyra::applyOp<Scalar>(assign_scalar_op,
91  ArrayView<Ptr<const MultiVectorBase<Scalar> > >(null),
92  tuple<Ptr<MultiVectorBase<Scalar> > >(ptr(this)), null);
93 }
94 
95 
96 template<class Scalar>
97 RCP<const MultiVectorBase<Scalar> >
99 {
100  using Teuchos::Workspace;
101  using Teuchos::as;
102  Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
103  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
104  const VectorSpaceBase<Scalar> &l_range = *this->range();
105  const Ordinal dimDomain = l_domain.dim();
106  const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
107  if( colRng.lbound() == 0 && as<Ordinal>(colRng.ubound()) == dimDomain-1 )
108  return Teuchos::rcp(this,false); // Takes all of the columns!
109  if( colRng.size() ) {
110  // We have to create a view of a subset of the columns
111  Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.size());
112  for( Ordinal j = colRng.lbound(); j <= colRng.ubound(); ++j )
113  col_vecs[j-colRng.lbound()] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(j));
114  return Teuchos::rcp(
116  this->range(),l_range.smallVecSpcFcty()->createVecSpc(colRng.size()),col_vecs
117  )
118  );
119  }
120  return Teuchos::null; // There was an empty set in colRng_in!
121 }
122 
123 
124 template<class Scalar>
125 RCP<MultiVectorBase<Scalar> >
127 {
128  using Teuchos::Workspace;
129  using Teuchos::as;
130  Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
131  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
132  const VectorSpaceBase<Scalar> &l_range = *this->range();
133  const Ordinal dimDomain = l_domain.dim();
134  const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
135  if( colRng.lbound() == 0 && as<Ordinal>(colRng.ubound()) == dimDomain-1 )
136  return Teuchos::rcp(this,false); // Takes all of the columns!
137  if( colRng.size() ) {
138  // We have to create a view of a subset of the columns
139  Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.size());
140  for( Ordinal j = colRng.lbound(); j <= colRng.ubound(); ++j )
141  col_vecs[j-colRng.lbound()] = this->col(j);
142  return Teuchos::rcp(
144  this->range(),l_range.smallVecSpcFcty()->createVecSpc(colRng.size()),col_vecs
145  )
146  );
147  }
148  return Teuchos::null; // There was an empty set in colRng_in!
149 }
150 
151 
152 template<class Scalar>
153 RCP<const MultiVectorBase<Scalar> >
155  const ArrayView<const int> &cols
156  ) const
157 {
158  using Teuchos::Workspace;
159  Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
160  const VectorSpaceBase<Scalar> &l_range = *this->range();
161  const int numCols = cols.size();
162 #ifdef TEUCHOS_DEBUG
163  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
164  const Ordinal dimDomain = l_domain.dim();
165  const char msg_err[] = "MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
166  TEUCHOS_TEST_FOR_EXCEPTION( numCols < 1 || dimDomain < numCols, std::invalid_argument, msg_err );
167 #endif
168  // We have to create a view of a subset of the columns
169  Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
170  for( int k = 0; k < numCols; ++k ) {
171  const int col_k = cols[k];
172 #ifdef TEUCHOS_DEBUG
173  TEUCHOS_TEST_FOR_EXCEPTION(
174  !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
175  ,msg_err << " col["<<k<<"] = " << col_k << " is not in the range [0,"<<(dimDomain-1)<<"]!"
176  );
177 #endif
178  col_vecs[k] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(col_k));
179  }
180  return Teuchos::rcp(
182  this->range(), l_range.smallVecSpcFcty()->createVecSpc(numCols), col_vecs
183  )
184  );
185 }
186 
187 
188 template<class Scalar>
189 RCP<MultiVectorBase<Scalar> >
191  const ArrayView<const int> &cols
192  )
193 {
194  using Teuchos::Workspace;
195  Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
196  const VectorSpaceBase<Scalar> &l_range = *this->range();
197  const int numCols = cols.size();
198 #ifdef TEUCHOS_DEBUG
199  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
200  const Ordinal dimDomain = l_domain.dim();
201  const char msg_err[] = "MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
202  TEUCHOS_TEST_FOR_EXCEPTION( numCols < 1 || dimDomain < numCols, std::invalid_argument, msg_err );
203 #endif
204  // We have to create a view of a subset of the columns
205  Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
206  for( int k = 0; k < numCols; ++k ) {
207  const int col_k = cols[k];
208 #ifdef TEUCHOS_DEBUG
209  TEUCHOS_TEST_FOR_EXCEPTION(
210  !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
211  ,msg_err << " col["<<k<<"] = " << col_k << " is not in the range [0,"<<(dimDomain-1)<<"]!"
212  );
213 #endif
214  col_vecs[k] = this->col(col_k);
215  }
216  return Teuchos::rcp(
218  this->range(), l_range.smallVecSpcFcty()->createVecSpc(numCols), col_vecs
219  )
220  );
221 }
222 
223 
224 template<class Scalar>
226  const RTOpPack::RTOpT<Scalar> &prim_op,
227  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs,
228  const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs,
229  const ArrayView<const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
230  const Ordinal prim_global_offset_in
231  ) const
232 {
233 
234  using Teuchos::Workspace;
235  using Teuchos::as;
236  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
237 
238  const int num_multi_vecs = multi_vecs.size();
239  const int num_targ_multi_vecs = targ_multi_vecs.size();
240 
241  // ToDo: Validate the input!
242 
243  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
244 
245  // Get the primary and secondary dimensions.
246 
247  const Ordinal sec_dim = l_domain.dim();
248 
249  //
250  // Apply the reduction/transformation operator and transform the
251  // target vectors and reduce each of the reduction objects.
252  //
253 
254  Workspace<RCP<const VectorBase<Scalar> > > vecs_s(wss, num_multi_vecs);
255  Workspace<Ptr<const VectorBase<Scalar> > > vecs(wss, num_multi_vecs);
256  Workspace<RCP<VectorBase<Scalar> > > targ_vecs_s(wss, num_targ_multi_vecs);
257  Workspace<Ptr<VectorBase<Scalar> > > targ_vecs(wss, num_targ_multi_vecs);
258 
259  for(Ordinal j = 0; j < sec_dim; ++j) {
260  // Fill the arrays of vector arguments
261  {for(Ordinal k = 0; k < as<Ordinal>(num_multi_vecs); ++k) {
262  vecs_s[k] = multi_vecs[k]->col(j);
263  vecs[k] = vecs_s[k].ptr();
264  }}
265  {for(Ordinal k = 0; k < as<Ordinal>(num_targ_multi_vecs); ++k) {
266  targ_vecs_s[k] = targ_multi_vecs[k]->col(j);
267  targ_vecs[k] = targ_vecs_s[k].ptr();
268  }}
269  // Apply the reduction/transformation operator
270  Thyra::applyOp(
271  prim_op,
272  vecs().getConst(),
273  targ_vecs().getConst(),
274  reduct_objs.size() ? reduct_objs[j] : Ptr<RTOpPack::ReductTarget>(),
275  prim_global_offset_in);
276  }
277  // At this point all of the designated targ vectors in the target multi-vectors have
278  // been transformed and all the reduction objects in reduct_obj[] have accumulated
279  // the reductions.
280 }
281 
282 
283 template<class Scalar>
285  const RTOpPack::RTOpT<Scalar> &prim_op,
286  const RTOpPack::RTOpT<Scalar> &sec_op,
287  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs,
288  const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs,
289  const Ptr<RTOpPack::ReductTarget> &reduct_obj,
290  const Ordinal prim_global_offset_in
291  ) const
292 {
293 
294  using Teuchos::Workspace;
295  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
296 
297  // ToDo: Validate the input!
298 
299  const VectorSpaceBase<Scalar> &l_domain = *this->domain();
300 
301  // Get the primary and secondary dimensions.
302  const Ordinal sec_dim = l_domain.dim();
303 
304  // Create a temporary buffer for the reduction objects of the primary reduction
305  // so that we can call the companion version of this method.
306  const int reduct_objs_size = (!is_null(reduct_obj) ? sec_dim : 0);
307  Workspace<RCP<RTOpPack::ReductTarget> > rcp_reduct_objs(wss, reduct_objs_size);
308  Workspace<Ptr<RTOpPack::ReductTarget> > reduct_objs(wss, reduct_objs_size);
309  if (!is_null(reduct_obj)) {
310  for(Ordinal k = 0; k < sec_dim; ++k) {
311  rcp_reduct_objs[k] = prim_op.reduct_obj_create();
312  reduct_objs[k] = rcp_reduct_objs[k].ptr();
313  }
314  }
315 
316  // Call the companion version that accepts an array of reduction objects
317  this->applyOp(
318  prim_op, multi_vecs, targ_multi_vecs, reduct_objs,
319  prim_global_offset_in);
320 
321  // Reduce all the reduction objects using the secondary reduction operator
322  // into one reduction object and free the intermediate reduction objects.
323  if (!is_null(reduct_obj)) {
324  for (Ordinal k = 0; k < sec_dim; ++k) {
325  sec_op.reduce_reduct_objs( *reduct_objs[k], reduct_obj );
326  }
327  }
328 }
329 
330 
331 template<class Scalar>
333  const Range1D &rowRng_in,
334  const Range1D &colRng_in,
336  ) const
337 {
338  const Ordinal
339  rangeDim = this->range()->dim(),
340  domainDim = this->domain()->dim();
341  const Range1D
342  rowRng = rowRng_in.full_range() ? Range1D(0,rangeDim-1) : rowRng_in,
343  colRng = colRng_in.full_range() ? Range1D(0,domainDim-1) : colRng_in;
344 #ifdef TEUCHOS_DEBUG
345  TEUCHOS_TEST_FOR_EXCEPTION(
346  !(rowRng.ubound() < rangeDim), std::out_of_range
347  ,"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, rowRng = ["
348  <<rowRng.lbound()<<","<<rowRng.ubound()<<"] is not in the range = [0,"
349  <<(rangeDim-1)<<"]!"
350  );
351  TEUCHOS_TEST_FOR_EXCEPTION(
352  !(colRng.ubound() < domainDim), std::out_of_range
353  ,"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, colRng = ["
354  <<colRng.lbound()<<","<<colRng.ubound()<<"] is not in the range = [0,"
355  <<(domainDim-1)<<"]!"
356  );
357 #endif
358  // Allocate storage for the multi-vector (stored column-major)
359  const ArrayRCP<Scalar> values = Teuchos::arcp<Scalar>(rowRng.size() * colRng.size());
360  // Extract multi-vector values column by column
361  RTOpPack::ConstSubVectorView<Scalar> sv; // uninitialized by default
362  for( int k = colRng.lbound(); k <= colRng.ubound(); ++k ) {
363  RCP<const VectorBase<Scalar> > col_k = this->col(k);
364  col_k->acquireDetachedView( rowRng, &sv );
365  for( int i = 0; i < rowRng.size(); ++i )
366  values[ i + k*rowRng.size() ] = sv[i];
367  col_k->releaseDetachedView( &sv );
368  }
369  // Initialize the multi-vector view object
370  sub_mv->initialize(
371  rowRng.lbound(), // globalOffset
372  rowRng.size(), // subDim
373  colRng.lbound(), // colOffset
374  colRng.size(), // numSubCols
375  values, // values
376  rowRng.size() // leadingDim
377  );
378 }
379 
380 
381 template<class Scalar>
384  ) const
385 {
386  // Here we just need to free the view and that is it!
387  sub_mv->uninitialize();
388 }
389 
390 
391 template<class Scalar>
393  const Range1D &rowRng,
394  const Range1D &colRng,
396  )
397 {
398  using Teuchos::as;
399  // Use the non-const implementation since it does exactly the
400  // correct thing in this case also!
402  rowRng, colRng,
404  // This cast will work as long as SubMultiVectorView
405  // maintains no extra state over ConstSubMultiVectorView (which it
406  // currently does not) but this is something that I should
407  // technically check for some how.
408  );
409 }
410 
411 
412 template<class Scalar>
415  )
416 {
417 #ifdef TEUCHOS_DEBUG
418  TEUCHOS_TEST_FOR_EXCEPTION(
419  sub_mv==NULL, std::logic_error,
420  "MultiVectorDefaultBase<Scalar>::commitNonconstDetachedMultiVectorViewImpl(...): Error!"
421  );
422 #endif
423  // Set back the multi-vector values column by column
424  const Range1D rowRng(sub_mv->globalOffset(),sub_mv->globalOffset()+sub_mv->subDim()-1);
425  RTOpPack::SubVectorView<Scalar> msv; // uninitialized by default
426  for( int k = sub_mv->colOffset(); k < sub_mv->numSubCols(); ++k ) {
427  RCP<VectorBase<Scalar> > col_k = this->col(k);
428  col_k->acquireDetachedView( rowRng, &msv );
429  for( int i = 0; i < rowRng.size(); ++i )
430  msv[i] = sub_mv->values()[ i + k*rowRng.size() ];
431  col_k->commitDetachedView( &msv );
432  }
433  // Zero out the view
434  sub_mv->uninitialize();
435 }
436 
437 
438 } // end namespace Thyra
439 
440 
441 #endif // THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
virtual void mvMultiReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const ArrayView< const Ptr< RTOpPack::ReductTarget > > &reduct_objs, const Ordinal primary_global_offset) const
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
virtual void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void mvSingleReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const RTOpPack::RTOpT< Scalar > &secondary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const Ptr< RTOpPack::ReductTarget > &reduct_obj, const Ordinal primary_global_offset) const
Abstract interface for objects that represent a space for vectors.
virtual void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
virtual void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
Teuchos::RCP< ReductTarget > reduct_obj_create() const
virtual RCP< MultiVectorBase< Scalar > > clone_mv() const
virtual RCP< const VectorSpaceFactoryBase< Scalar > > smallVecSpcFcty() const =0
Return a VectorSpaceFactoryBase object for the creation of (usually serial) vector spaces with a smal...
Abstract interface for finite-dimensional dense vectors.
Default subclass for MultiVectorBase implemented using columns of separate abstract vectors...
virtual void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
const ArrayRCP< Scalar > values() const
void initialize(Ordinal globalOffset_in, Ordinal subDim_in, Ordinal colOffset_in, Ordinal numSubCols_in, const ArrayRCP< const Scalar > &values_in, Ordinal leadingDim_in)
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols) const
void reduce_reduct_objs(const ReductTarget &in_reduct_obj, const Ptr< ReductTarget > &inout_reduct_obj) const
virtual Ordinal dim() const =0
Return the dimension of the vector space.
Teuchos::Range1D Range1D