Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultMultiVectorProductVector_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_MULTI_VECTOR_PRODUCT_VECTOR_HPP
43#define THYRA_DEFAULT_MULTI_VECTOR_PRODUCT_VECTOR_HPP
44
45
46#include "Thyra_DefaultMultiVectorProductVector_decl.hpp"
47#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
48#include "Thyra_AssertOp.hpp"
49#include "Teuchos_Assert.hpp"
50
51
52namespace Thyra {
53
54
55// Constructors/initializers/accessors
56
57
58template <class Scalar>
63
64
65template <class Scalar>
67 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &productSpace_in,
68 const RCP<MultiVectorBase<Scalar> > &multiVec
69 )
70{
71#ifdef TEUCHOS_DEBUG
72 TEUCHOS_TEST_FOR_EXCEPT(is_null(productSpace_in));
73 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVec));
75 "DefaultMultiVectorProductVector<Scalar>::initialize(productSpace,multiVec)",
76 *multiVec->range(), *productSpace_in->getBlock(0)
77 );
78 TEUCHOS_ASSERT_EQUALITY( multiVec->domain()->dim(), productSpace_in->numBlocks());
79#endif
80
81 numBlocks_ = productSpace_in->numBlocks();
82
83 productSpace_ = productSpace_in;
84
85 multiVec_ = multiVec;
86
87}
88
89
90template <class Scalar>
92 const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &productSpace_in,
93 const RCP<const MultiVectorBase<Scalar> > &multiVec
94 )
95{
96#ifdef TEUCHOS_DEBUG
97 TEUCHOS_TEST_FOR_EXCEPT(is_null(productSpace_in));
98 TEUCHOS_TEST_FOR_EXCEPT(is_null(multiVec));
100 "DefaultMultiVectorProductVector<Scalar>::initialize(productSpace_in,multiVec)",
101 *multiVec->range(), *productSpace_in->getBlock(0)
102 );
103 TEUCHOS_ASSERT_EQUALITY( multiVec->domain()->dim(), productSpace_in->numBlocks() );
104#endif
105
106 numBlocks_ = productSpace_in->numBlocks();
107
108 productSpace_ = productSpace_in;
109
110 multiVec_ = multiVec;
111
112}
113
114
115template <class Scalar>
118{
119 return multiVec_.getNonconstObj();
120}
121
122
123template <class Scalar>
126{
127 return multiVec_.getConstObj();
128}
129
130
131template <class Scalar>
133{
134 numBlocks_ = 0;
135 productSpace_ = Teuchos::null;
136 multiVec_.uninitialize();
137}
138
139
140// Overridden from Teuchos::Describable
141
142
143template<class Scalar>
145{
146 std::ostringstream oss;
147 oss
149 << "{"
150 << "dim="<<this->space()->dim()
151 << ",numColumns = "<<numBlocks_
152 << "}";
153 return oss.str();
154}
155
156template<class Scalar>
158 Teuchos::FancyOStream &out_arg,
159 const Teuchos::EVerbosityLevel verbLevel
160 ) const
161{
162 using Teuchos::OSTab;
163 using Teuchos::describe;
164 RCP<FancyOStream> out = rcp(&out_arg,false);
165 OSTab tab(out);
166 switch(verbLevel) {
169 *out << this->description() << std::endl;
170 break;
174 {
175 *out
177 << "dim=" << this->space()->dim()
178 << "}\n";
179 OSTab tab2(out);
180 *out << "multiVec = " << Teuchos::describe(*multiVec_.getConstObj(),verbLevel);
181 break;
182 }
183 default:
184 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
185 }
186}
187
188
189// Overridden from ProductVectorBase
190
191
192template <class Scalar>
195{
196#ifdef TEUCHOS_DEBUG
197 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( k, 0, numBlocks_ );
198#endif
199 return multiVec_.getNonconstObj()->col(k);
200}
201
202
203template <class Scalar>
206{
207#ifdef TEUCHOS_DEBUG
208 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( k, 0, numBlocks_ );
209#endif
210 return multiVec_.getConstObj()->col(k);
211}
212
213
214// Overridden from ProductMultiVectorBase
215
216
217template <class Scalar>
220{
221 return productSpace_;
222}
223
224
225template <class Scalar>
227{
228#ifdef TEUCHOS_DEBUG
229 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( k, 0, numBlocks_ );
230#else
231 (void)k;
232#endif
233 return multiVec_.isConst();
234}
235
236
237template <class Scalar>
240{
241 return getNonconstVectorBlock(k);
242}
243
244
245template <class Scalar>
248{
249 return getVectorBlock(k);
250}
251
252
253// Overridden public functions from VectorBase
254
255
256template <class Scalar>
259{
260 return productSpace_;
261}
262
263
264// protected
265
266
267// Overridden protected functions from VectorBase
268
269
270template <class Scalar>
272 Scalar l,
273 Scalar u
274 )
275{
276 for(int k = 0; k < numBlocks_; ++k) {
277 multiVec_.getNonconstObj()->col(k)->randomize(l, u);
278 }
279}
280
281
282template <class Scalar>
284 const VectorBase<Scalar>& x
285 )
286{
287#ifdef TEUCHOS_DEBUG
288 TEUCHOS_ASSERT(productSpace_->isCompatible(*x.space()));
289#endif
290 // Apply operation block-by-block if mv is also a ProductVector
293 if (nonnull(pv)) {
294 for(int k = 0; k < numBlocks_; ++k) {
295 multiVec_.getNonconstObj()->col(k)->abs(*pv->getVectorBlock(k));
296 }
297 } else {
299 }
300}
301
302
303template <class Scalar>
305 const VectorBase<Scalar>& x
306 )
307{
308#ifdef TEUCHOS_DEBUG
309 TEUCHOS_ASSERT(productSpace_->isCompatible(*x.space()));
310#endif
311 // Apply operation block-by-block if mv is also a ProductVector
314 if (nonnull(pv)) {
315 for(int k = 0; k < numBlocks_; ++k) {
316 multiVec_.getNonconstObj()->col(k)->reciprocal(*pv->getVectorBlock(k));
317 }
318 } else {
320 }
321}
322
323
324template <class Scalar>
326 const VectorBase<Scalar>& x
327 )
328{
329#ifdef TEUCHOS_DEBUG
330 TEUCHOS_ASSERT(productSpace_->isCompatible(*x.space()));
331#endif
332 // Apply operation block-by-block if mv is also a ProductVector
335 if (nonnull(pv)) {
336 for(int k = 0; k < numBlocks_; ++k) {
337 multiVec_.getNonconstObj()->col(k)->ele_wise_scale(*pv->getVectorBlock(k));
338 }
339 } else {
341 }
342}
343
344
345template <class Scalar>
348 const VectorBase<Scalar>& x
349 ) const
350{
351#ifdef TEUCHOS_DEBUG
352 TEUCHOS_ASSERT(productSpace_->isCompatible(*x.space()));
353#endif
354 // Apply operation block-by-block if mv is also a ProductVector
355 typedef ScalarTraits<Scalar> ST;
358 if (nonnull(pv)) {
359 typename ST::magnitudeType norm = ST::magnitude(ST::zero());
360 for(int k = 0; k < numBlocks_; ++k) {
361 typename ST::magnitudeType subNorm
362 = multiVec_.getConstObj()->col(k)->norm_2(*pv->getVectorBlock(k));
363 norm += subNorm*subNorm;
364 }
365 return ST::magnitude(ST::squareroot(norm));
366 } else {
368 }
369}
370
371
372template <class Scalar>
374 const RTOpPack::RTOpT<Scalar> &op,
375 const ArrayView<const Ptr<const VectorBase<Scalar> > > &vecs,
376 const ArrayView<const Ptr<VectorBase<Scalar> > > &targ_vecs,
377 const Ptr<RTOpPack::ReductTarget> &reduct_obj,
378 const Ordinal global_offset
379 ) const
380{
381 this->getDefaultProductVector()->applyOp(
382 op, vecs, targ_vecs, reduct_obj, global_offset );
383}
384
385
386template <class Scalar>
388 const Range1D& rng_in, RTOpPack::ConstSubVectorView<Scalar>* sub_vec
389 ) const
390{
391 this->getDefaultProductVector()->acquireDetachedView(rng_in,sub_vec);
392}
393
394
395template <class Scalar>
398 ) const
399{
400 this->getDefaultProductVector()->releaseDetachedView(sub_vec);
401}
402
403
404template <class Scalar>
406 const Range1D& /* rng_in */, RTOpPack::SubVectorView<Scalar>* /* sub_vec */
407 )
408{
409 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement DefaultMultiVectorProductVector<Scalar>::acquireNonconstDetachedVectorViewImpl(...)!");
410}
411
412
413template <class Scalar>
416 )
417{
418 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement DefaultMultiVectorProductVector<Scalar>::commitNonconstDetachedVectorViewImpl(...)!");
419}
420
421
422template <class Scalar>
424 const RTOpPack::SparseSubVectorT<Scalar>& /* sub_vec */
425 )
426{
427 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement DefaultMultiVectorProductVector<Scalar>::setSubVector(...)!");
428}
429
430
431// Overridden protected functions from VectorBase
432
433
434template <class Scalar>
436{
437 multiVec_.getNonconstObj()->assign(alpha);
438}
439
440
441template <class Scalar>
444 )
445{
446#ifdef TEUCHOS_DEBUG
447 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
448 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
449#endif
452 if (nonnull(pv)) {
453 for(int k = 0; k < numBlocks_; ++k) {
454 multiVec_.getNonconstObj()->col(k)->assign(*pv->getMultiVectorBlock(k));
455 }
456 } else {
458 }
459}
460
461
462template <class Scalar>
464{
465 multiVec_.getNonconstObj()->scale(alpha);
466}
467
468
469template <class Scalar>
471 Scalar alpha,
473 )
474{
475#ifdef TEUCHOS_DEBUG
476 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
477 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
478#endif
481 if (nonnull(pv)) {
482 for(int k = 0; k < numBlocks_; ++k) {
483 multiVec_.getNonconstObj()->col(k)->update(alpha, *pv->getMultiVectorBlock(k));
484 }
485 } else {
487 }
488}
489
490
491template <class Scalar>
493 const ArrayView<const Scalar>& alpha,
494 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
495 const Scalar& beta
496 )
497{
498#ifdef TEUCHOS_DEBUG
499 TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
500 for (Ordinal i = 0; i < mv.size(); ++i) {
501 TEUCHOS_ASSERT_EQUALITY(mv[i]->domain()->dim(), 1);
502 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv[i]->range()));
503 }
504#endif
505
506 // Apply operation block-by-block if each element of mv is also a ProductMultiVector
507 bool allCastsSuccessful = true;
509 for (Ordinal i = 0; i < mv.size(); ++i) {
511 if (pvs[i].is_null()) {
512 allCastsSuccessful = false;
513 break;
514 }
515 }
516
517 if (allCastsSuccessful) {
518 Array<RCP<const MultiVectorBase<Scalar> > > blocks_rcp(pvs.size());
520 for (int k = 0; k < numBlocks_; ++k) {
521 for (Ordinal i = 0; i < pvs.size(); ++i) {
522 blocks_rcp[i] = pvs[i]->getMultiVectorBlock(k);
523 blocks[i] = blocks_rcp[i].ptr();
524 }
525 multiVec_.getNonconstObj()->col(k)->linear_combination(alpha, blocks(), beta);
526 }
527 } else {
529 }
530}
531
532
533template <class Scalar>
535 const MultiVectorBase<Scalar>& mv,
536 const ArrayView<Scalar>& prods
537 ) const
538{
539#ifdef TEUCHOS_DEBUG
540 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
541 TEUCHOS_ASSERT_EQUALITY(prods.size(), 1);
542 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
543#endif
546 if (nonnull(pv)) {
547 prods[0] = ScalarTraits<Scalar>::zero();
548 for(int k = 0; k < numBlocks_; ++k) {
549 Scalar prod;
550 multiVec_.getConstObj()->col(k)->dots(*pv->getMultiVectorBlock(k), Teuchos::arrayView(&prod, 1));
551 prods[0] += prod;
552 }
553 } else {
555 }
556}
557
558
559template <class Scalar>
562 ) const
563{
564#ifdef TEUCHOS_DEBUG
565 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
566#endif
567 typedef ScalarTraits<Scalar> ST;
568 norms[0] = ST::magnitude(ST::zero());
569 for(int k = 0; k < numBlocks_; ++k) {
570 norms[0] += multiVec_.getConstObj()->col(k)->norm_1();
571 }
572}
573
574
575template <class Scalar>
578 ) const
579{
580#ifdef TEUCHOS_DEBUG
581 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
582#endif
583 typedef ScalarTraits<Scalar> ST;
584 norms[0] = ST::magnitude(ST::zero());
585 for(int k = 0; k < numBlocks_; ++k) {
586 typename ST::magnitudeType subNorm = multiVec_.getConstObj()->col(k)->norm_2();
587 norms[0] += subNorm*subNorm;
588 }
589 norms[0] = ST::magnitude(ST::squareroot(norms[0]));
590}
591
592
593template <class Scalar>
596 ) const
597{
598#ifdef TEUCHOS_DEBUG
599 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
600#endif
601 typedef ScalarTraits<Scalar> ST;
602 norms[0] = ST::magnitude(ST::zero());
603 for(int k = 0; k < numBlocks_; ++k) {
604 typename ST::magnitudeType subNorm = multiVec_.getConstObj()->col(k)->norm_inf();
605 if (subNorm > norms[0])
606 norms[0] = subNorm;
607 }
608}
609
610
611// private
612
613
614template <class Scalar>
617{
618
619 // This function exists since in general we can not create views of a column
620 // vectors and expect the changes to be mirrored in the mulit-vector
621 // automatically. Later, we might be able to change this once we have a
622 // Thyra::MultiVectorBase::hasDirectColumnVectorView() function and it
623 // returns true. Until then, this is the safe way to do this ...
624
626 for ( int k = 0; k < numBlocks_; ++k) {
627 vecArray.push_back(multiVec_.getConstObj()->col(k));
628 }
629
630 return Thyra::defaultProductVector<Scalar>(
631 productSpace_->getDefaultProductVectorSpace(),
632 vecArray()
633 );
634
635}
636
637
638} // namespace Thyra
639
640
641#endif // THYRA_DEFAULT_MULTI_VECTOR_PRODUCT_VECTOR_HPP
size_type size() const
size_type size() const
void push_back(const value_type &x)
virtual std::string description() const
Standard concrete implementation of a product vector space that creates product vectors fromed implic...
Concrete implementation of a product vector which is really composed out of the columns of a multi-ve...
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
void initialize(const RCP< const DefaultMultiVectorProductVectorSpace< Scalar > > &productSpace, const RCP< MultiVectorBase< Scalar > > &multiVec)
Initialize with a non-const multi-vector.
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
virtual void eleWiseScaleImpl(const VectorBase< Scalar > &x)
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)
RCP< const VectorSpaceBase< Scalar > > space() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const
virtual void reciprocalImpl(const VectorBase< Scalar > &x)
void releaseDetachedVectorViewImpl(RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
void setSubVectorImpl(const RTOpPack::SparseSubVectorT< Scalar > &sub_vec)
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
RCP< VectorBase< Scalar > > getNonconstVectorBlock(const int k)
void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
virtual Teuchos::ScalarTraits< Scalar >::magnitudeType norm2WeightedImpl(const VectorBase< Scalar > &x) const
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Default implementation of dots using RTOps.
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Default implementation of linear_combination using RTOps.
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
Default implementation of assign(MV) using RTOps.
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Default implementation of update using RTOps.
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.
virtual Teuchos::ScalarTraits< Scalar >::magnitudeType norm2WeightedImpl(const VectorBase< Scalar > &x) const
Default implementation of norm_2 (weighted) using RTOps.
virtual void reciprocalImpl(const VectorBase< Scalar > &x)
Default implementation of reciprocal using RTOps.
virtual void absImpl(const VectorBase< Scalar > &x)
Default implementation of abs using RTOps.
virtual void eleWiseScaleImpl(const VectorBase< Scalar > &x)
Default implementation of ele_wise_scale using RTOps.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
T_To & dyn_cast(T_From &from)