Xpetra Version of the Day
Xpetra_BlockedMultiVector_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Xpetra: A linear algebra interface package
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Tobias Wiesner (tawiesn@sandia.gov)
42// Ray Tuminaro (rstumin@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
47#ifndef XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
48#define XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
49
51
52#include "Xpetra_MultiVectorFactory.hpp"
53#include "Xpetra_BlockedVector.hpp"
54#include "Xpetra_MapExtractor.hpp"
55
56
57namespace Xpetra {
58
59
60template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63 size_t NumVectors,
64 bool zeroOut)
65 : map_(map)
66{
67 numVectors_ = NumVectors;
68
69 vv_.reserve(map->getNumMaps());
70
71 // add CrsMatrix objects in row,column order
72 for(size_t r = 0; r < map->getNumMaps(); ++r)
73 vv_.push_back(
74 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map->getMap(r, map_->getThyraMode()), NumVectors, zeroOut));
75}
76
77
89template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93{
94 XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(),
96 "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has "
97 << bmap->getMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements()
98 << ".");
99 XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(),
101 "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has "
102 << bmap->getMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements()
103 << ".");
104 // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError,
105 // "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " <<
106 // bmap->getFullMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements() << ".");
107 // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError,
108 // "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has " <<
109 // bmap->getFullMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements() << ".");
110
111 numVectors_ = v->getNumVectors();
112
113 map_ = bmap;
114
115 // Extract vector
116 vv_.reserve(bmap->getNumMaps());
117
118 // add CrsMatrix objects in row,column order
119 for(size_t r = 0; r < bmap->getNumMaps(); ++r)
120 vv_.push_back(this->ExtractVector(v, r, bmap->getThyraMode()));
121}
122
123
135template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
139{
140 numVectors_ = v->getNumVectors();
141
142 // create blocked map out of MapExtractor
143 std::vector<RCP<const Map>> maps;
144 maps.reserve(mapExtractor->NumMaps());
145 for(size_t r = 0; r < mapExtractor->NumMaps(); ++r)
146 maps.push_back(mapExtractor->getMap(r, mapExtractor->getThyraMode()));
147 map_ = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapExtractor->getFullMap(), maps, mapExtractor->getThyraMode()));
148
149 // Extract vector
150 vv_.reserve(mapExtractor->NumMaps());
151
152 // add CrsMatrix objects in row,column order
153 for(size_t r = 0; r < mapExtractor->NumMaps(); ++r)
154 vv_.push_back(this->ExtractVector(v, r, mapExtractor->getThyraMode()));
155}
156
157
158template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
161 std::vector<Teuchos::RCP<MultiVector>>& vin)
162{
163 numVectors_ = vin[ 0 ]->getNumVectors();
164 map_ = map;
165 vv_.resize(vin.size());
166 for(size_t i = 0; i < vv_.size(); i++)
167 vv_[ i ] = vin[ i ];
169
170
171template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175 for(size_t r = 0; r < vv_.size(); ++r)
176 {
177 vv_[ r ] = Teuchos::null;
178 }
179
180 map_ = Teuchos::null;
181 numVectors_ = 0;
182}
183
184
185template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
188operator=(const MultiVector& rhs)
190 assign(rhs); // dispatch to protected virtual method
191 return *this;
192}
194
195template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
196void
198replaceGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */)
199{
200 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceGlobalValue: Not (yet) supported by BlockedMultiVector.");
202
203
204template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205void
207sumIntoGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */)
208{
209 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoGlobalValue: Not (yet) supported by BlockedMultiVector.");
211
212template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213void
215replaceLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */)
216{
217 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceLocalValue: Not supported by BlockedMultiVector.");
219
220
221template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
222void
224sumIntoLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */)
225{
226 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoLocalValue:Not (yet) supported by BlockedMultiVector.");
227}
228
229
231template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
232void
234putScalar(const Scalar& value)
235{
236 XPETRA_MONITOR("BlockedMultiVector::putScalar");
237 for(size_t r = 0; r < map_->getNumMaps(); r++)
239 getMultiVector(r)->putScalar(value);
240 }
241}
243
244template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
247getVector(size_t j) const
248{
251
252 for(size_t r = 0; r < getBlockedMap()->getNumMaps(); r++)
253 {
255 this->getMultiVector(r, this->getBlockedMap()->getThyraMode())->getVector(j);
256 ret->setMultiVector(r, subvec, this->getBlockedMap()->getThyraMode());
257 }
258 return ret;
259}
260
261
262template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
265getVectorNonConst(size_t j)
266{
268 Teuchos::rcp_const_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(getVector(j));
269 return ret;
270}
272
273template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
276getData(size_t j) const
277{
278 if(map_->getNumMaps() == 1)
280 return vv_[ 0 ]->getData(j);
281 }
282 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getData: Not (yet) supported by BlockedMultiVector.");
283 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
284}
285
286
287template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290getDataNonConst(size_t j)
291{
292 if(map_->getNumMaps() == 1)
293 {
294 return vv_[ 0 ]->getDataNonConst(j);
296 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getDataNonConst: Not (yet) supported by BlockedMultiVector.");
297 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
298}
300
301template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
302void
304dot(const MultiVector& /* A */, const Teuchos::ArrayView<Scalar>& /* dots */) const
305{
306 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::dot: Not (yet) supported by BlockedMultiVector.");
308
309
310template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
311void
313abs(const MultiVector& /* A */)
314{
315 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::abs: Not (yet) supported by BlockedMultiVector.");
316}
317
318
319template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
320void
322reciprocal(const MultiVector& /* A */)
324 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::reciprocal: Not (yet) supported by BlockedMultiVector.");
325}
326
328template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329void
331scale(const Scalar& alpha)
332{
333 XPETRA_MONITOR("BlockedMultiVector::scale (Scalar)");
334 for(size_t r = 0; r < map_->getNumMaps(); ++r)
335 {
336 if(getMultiVector(r) != Teuchos::null)
337 {
338 getMultiVector(r)->scale(alpha);
339 }
340 }
341}
343
344template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
345void
348{
349 XPETRA_MONITOR("BlockedMultiVector::scale (Array)");
350 for(size_t r = 0; r < map_->getNumMaps(); ++r)
351 {
352 if(getMultiVector(r) != Teuchos::null)
353 {
354 getMultiVector(r)->scale(alpha);
355 }
356 }
357}
358
359
360template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
361void
363update(const Scalar& alpha, const MultiVector& A, const Scalar& beta)
364{
365 XPETRA_MONITOR("BlockedMultiVector::update");
366 Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
367 Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
368 TEUCHOS_TEST_FOR_EXCEPTION(numVectors_ != rcpA->getNumVectors(),
370 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector).");
371 if(bA != Teuchos::null)
372 {
373 // A is a BlockedMultiVector (and compatible with this)
374 // Call update recursively on all sub vectors
375 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(),
377 "BlockedMultiVector::update: update with incompatible vector (different thyra mode).");
378 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(),
380 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors).");
381 for(size_t r = 0; r < map_->getNumMaps(); r++)
384 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != bA->getMultiVector(r)->getMap()->getNodeNumElements(),
386 "BlockedMultiVector::update: in subvector "
387 << r << ": Cannot add a vector of (local) length " << bA->getMultiVector(r)->getMap()->getNodeNumElements()
388 << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
389 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != bA->getMultiVector(r)->getMap()->getGlobalNumElements(),
391 "BlockedMultiVector::update: in subvector "
392 << r << ": Cannot add a vector of length " << bA->getMultiVector(r)->getMap()->getGlobalNumElements()
393 << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
394
395 // TAW: 12/6 We basically want to do something like:
396 // getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta);
397 // But if the left hand side is a standard MultiVector and bA->getMultiVector(r) is
398 // a blocked MultiVector this will not work, as the update implementation of the standard
399 // multivector cannot deal with blocked multivectors.
400 // The only use case where this could happen is, if the rhs vector is a single block multivector
401 Teuchos::RCP<MultiVector> lmv = getMultiVector(r);
402 Teuchos::RCP<MultiVector> rmv = bA->getMultiVector(r);
403
404 // check whether lmv/rmv are blocked or not
405 Teuchos::RCP<BlockedMultiVector> blmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(lmv);
406 Teuchos::RCP<BlockedMultiVector> brmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rmv);
407
408 if(blmv.is_null() == true && brmv.is_null() == false)
409 {
410 // special case: lmv is standard MultiVector but rmv is BlockedMultiVector
411 TEUCHOS_TEST_FOR_EXCEPTION(brmv->getBlockedMap()->getNumMaps() > 1,
413 "BlockedMultiVector::update: Standard MultiVector object does not accept BlockedMultVector object as "
414 "parameter in update call.");
415 lmv->update(alpha, *(brmv->getMultiVector(0)), beta);
416 }
417 else
418 lmv->update(alpha, *rmv, beta);
419 }
420 }
421 else
422 {
423 // A is a MultiVector
424 // If this is a BlockedMultiVector with only one sub-vector of same length we can just update
425 // Otherwise, A is not compatible with this as BlockedMultiVector and we have to extract the vector data from A
426
427 if(getBlockedMap()->getNumMaps() == 1)
428 {
429 // Actually call "update" on the underlying MultiVector (Epetra or Tpetra).
430 // The maps have to be compatible.
432 getMultiVector(0)->getMap()->isSameAs(*(rcpA->getMap())) == false,
434 "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor).");
435 getMultiVector(0)->update(alpha, *rcpA, beta);
436 }
437 else
438 {
439 // general case: A has to be splitted and subparts have to be extracted and stored in this BlockedMultiVector
440 // XPETRA_TEST_FOR_EXCEPTION(map_->getFullMap()->isSameAs(*(rcpA->getMap()))==false, Xpetra::Exceptions::RuntimeError,
441 // "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor). - Note:
442 // This test might be too strict and can probably be relaxed!");
443 for(size_t r = 0; r < map_->getNumMaps(); r++)
444 {
445 // Call "update" on the subvector. Note, that getMultiVector(r) could return another BlockedMultiVector.
446 // That is, in Thyra mode the maps could differ (local Xpetra versus Thyra style gids)
447 Teuchos::RCP<const MultiVector> part = this->ExtractVector(rcpA, r, map_->getThyraMode());
448 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != part->getMap()->getNodeNumElements(),
450 "BlockedMultiVector::update: in subvector "
451 << r << ": Cannot add a vector of (local) length " << part->getMap()->getNodeNumElements()
452 << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
453 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != part->getMap()->getGlobalNumElements(),
455 "BlockedMultiVector::update: in subvector "
456 << r << ": Cannot add a vector of length " << part->getMap()->getGlobalNumElements()
457 << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
458 getMultiVector(r)->update(alpha, *part, beta);
459 }
460 }
461 }
462}
463
464
465template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
466void
468update(const Scalar& alpha, const MultiVector& A, const Scalar& beta, const MultiVector& B, const Scalar& gamma)
469{
470 XPETRA_MONITOR("BlockedMultiVector::update2");
471 Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
472 Teuchos::RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
473 Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
474 Teuchos::RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
475 if(bA != Teuchos::null && bB != Teuchos::null)
476 {
477 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(),
479 "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector A).");
480 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(),
482 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector A).");
484 numVectors_ != bA->getNumVectors(),
486 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector A).");
487 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bB->getBlockedMap()->getThyraMode(),
489 "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector B).");
490 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bB->getBlockedMap()->getNumMaps(),
492 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector B).");
494 numVectors_ != bB->getNumVectors(),
496 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector B).");
497
498 for(size_t r = 0; r < map_->getNumMaps(); r++)
499 {
500 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->isSameAs(*(bA->getMultiVector(r)->getMap())) == false,
502 "BlockedMultiVector::update: update with incompatible vector (different maps in partial vector " << r << ").");
503 getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta, *(bB->getMultiVector(r)), gamma);
504 }
505 return;
506 }
507 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::update: only supports update with other BlockedMultiVector.");
508}
509
510
511template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
512void
515{
516 XPETRA_MONITOR("BlockedMultiVector::norm1");
517 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
518 Array<Magnitude> temp_norms(getNumVectors());
519 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
520 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
521 for(size_t r = 0; r < map_->getNumMaps(); ++r)
522 {
523 if(getMultiVector(r) != Teuchos::null)
524 {
525 getMultiVector(r)->norm1(temp_norms);
526 for(size_t c = 0; c < getNumVectors(); ++c)
527 norms[ c ] += temp_norms[ c ];
528 }
529 }
530}
531
532
533
534template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
535void
538{
539 XPETRA_MONITOR("BlockedMultiVector::norm2");
540 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
541 Array<Magnitude> results(getNumVectors());
542 Array<Magnitude> temp_norms(getNumVectors());
543 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
544 std::fill(results.begin(), results.end(), ScalarTraits<Magnitude>::zero());
545 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
546 for(size_t r = 0; r < map_->getNumMaps(); ++r)
547 {
548 if(getMultiVector(r) != Teuchos::null)
549 {
550 getMultiVector(r)->norm2(temp_norms);
551 for(size_t c = 0; c < getNumVectors(); ++c)
552 results[ c ] += temp_norms[ c ] * temp_norms[ c ];
553 }
554 }
555 for(size_t c = 0; c < getNumVectors(); ++c)
556 norms[ c ] = Teuchos::ScalarTraits<Magnitude>::squareroot(results[ c ]);
557}
558
559
560template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
561void
564{
565 XPETRA_MONITOR("BlockedMultiVector::normInf");
566 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
567 Array<Magnitude> temp_norms(getNumVectors());
568 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
569 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
570 for(size_t r = 0; r < map_->getNumMaps(); ++r)
571 {
572 if(getMultiVector(r) != Teuchos::null)
573 {
574 getMultiVector(r)->normInf(temp_norms);
575 for(size_t c = 0; c < getNumVectors(); ++c)
576 norms[ c ] = std::max(norms[ c ], temp_norms[ c ]);
577 }
578 }
579}
580
581
582template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
583void
585meanValue(const Teuchos::ArrayView<Scalar>& /* means */) const
586{
587 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::meanValue: Not (yet) supported by BlockedMultiVector.");
588}
589
590
591template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
592void
594multiply(Teuchos::ETransp /* transA */,
595 Teuchos::ETransp /* transB */,
596 const Scalar& /* alpha */,
597 const MultiVector& /* A */,
598 const MultiVector& /* B */,
599 const Scalar& /* beta */)
600{
601 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::multiply: Not (yet) supported by BlockedMultiVector.");
602}
603
604
605template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
606void
608elementWiseMultiply(Scalar scalarAB,
610 const MultiVector& B,
611 Scalar scalarThis)
612{
613 XPETRA_MONITOR("BlockedMultiVector::elementWiseMultiply");
614 XPETRA_TEST_FOR_EXCEPTION(B.getMap()->isSameAs(*(this->getMap())) == false,
616 "BlockedMultiVector::elementWiseMultipy: B must have same blocked map than this.");
617 // XPETRA_TEST_FOR_EXCEPTION(A.getMap()->isSameAs(*(this->getMap()))==false, Xpetra::Exceptions::RuntimeError,
618 // "BlockedMultiVector::elementWiseMultipy: A must have same blocked map than this.");
619 TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getNodeNumElements() != B.getMap()->getNodeNumElements(),
621 "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getNodeNumElements() << " elements, B has "
622 << B.getMap()->getNodeNumElements() << ".");
623 TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getGlobalNumElements() != B.getMap()->getGlobalNumElements(),
625 "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getGlobalNumElements() << " elements, B has "
626 << B.getMap()->getGlobalNumElements() << ".");
627
628 RCP<const BlockedMap> bmap = getBlockedMap();
631 Teuchos::rcp_dynamic_cast<const Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(rcpA);
632
633 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
634 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
636 bB.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: B must be a BlockedMultiVector.");
637
638 // RCP<Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > me = Teuchos::rcp(new
639 // Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(bmap));
640
641 for(size_t m = 0; m < bmap->getNumMaps(); m++)
642 {
643 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partA = bA->getMultiVector(m, bmap->getThyraMode())->getVector(0);
644 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partB = bB->getMultiVector(m, bmap->getThyraMode());
645 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thisPart = this->getMultiVector(m, bmap->getThyraMode());
646
647 thisPart->elementWiseMultiply(scalarAB, *partA, *partB, scalarThis);
648 }
649}
650
651
652
653template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
654size_t
656getNumVectors() const
657{
658 return numVectors_;
659}
660
661
662template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
663size_t
665getLocalLength() const
666{
667 XPETRA_MONITOR("BlockedMultiVector::getLocalLength()");
668 return map_->getFullMap()->getNodeNumElements();
669}
670
671
672template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
675getGlobalLength() const
676{
677 XPETRA_MONITOR("BlockedMultiVector::getGlobalLength()");
678 return map_->getFullMap()->getGlobalNumElements();
679}
680
681
682template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
683bool
686{
687 const BlockedMultiVector* Vb = dynamic_cast<const BlockedMultiVector*>(&vec);
688 if(!Vb)
689 return false;
690 for(size_t r = 0; r < map_->getNumMaps(); ++r)
691 {
692 RCP<const MultiVector> a = getMultiVector(r);
694 if((a == Teuchos::null && b != Teuchos::null) || (a != Teuchos::null && b == Teuchos::null))
695 return false;
696 if(a != Teuchos::null && b != Teuchos::null && !a->isSameSize(*b))
697 return false;
698 }
699 return true;
700}
701
702
703template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
704std::string
706description() const
707{
708 return std::string("BlockedMultiVector");
709}
710
711
712template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
713void
716{
717 out << description() << std::endl;
718 for(size_t r = 0; r < map_->getNumMaps(); r++)
719 getMultiVector(r)->describe(out, verbLevel);
720}
721
722
723template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
724void
726replaceMap(const RCP<const Map>& map)
727{
728 XPETRA_MONITOR("BlockedMultiVector::replaceMap");
729 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
730 if(bmap.is_null() == true)
731 {
732 // if this has more than 1 sub blocks but "map" is not a blocked map, they are very likely not compatible
733 if(this->getBlockedMap()->getNumMaps() > 1)
734 {
736 "BlockedMultiVector::replaceMap: map is not of type BlockedMap. General implementation not available, yet.");
738 }
739 // special case: this is a blocked map with only one block
740 // TODO add more debug check (especially for Thyra mode)
741 std::vector<Teuchos::RCP<const Map>> subMaps(1, map);
742 map_ = Teuchos::rcp(new BlockedMap(map, subMaps, this->getBlockedMap()->getThyraMode()));
743 this->getMultiVector(0)->replaceMap(map);
744 return;
745 }
746 RCP<const BlockedMap> mybmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map_);
747
749 mybmap->getThyraMode() != bmap->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::replaceMap: inconsistent Thyra mode");
750 map_ = bmap;
751 for(size_t r = 0; r < map_->getNumMaps(); r++)
752 getMultiVector(r)->replaceMap(bmap->getMap(r, map_->getThyraMode()));
753}
754
755
756template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
757void
760 const Import& /* importer */,
761 CombineMode /* CM */)
762{
763 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
764}
765
766
767template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
768void
771 const Import& /* importer */,
772 CombineMode /* CM */)
773{
774 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
775}
776
777
778template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
779void
782 const Export& /* exporter */,
783 CombineMode /* CM */)
784{
785 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
786}
787
788
789template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
790void
793 const Export& /* exporter */,
794 CombineMode /* CM */)
795{
796 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
797}
798
799
800template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
801void
803setSeed(unsigned int seed)
804{
805 for(size_t r = 0; r < map_->getNumMaps(); ++r)
806 {
807 getMultiVector(r)->setSeed(seed);
808 }
809}
810
811
812template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
813void
815randomize(bool bUseXpetraImplementation)
816{
817 for(size_t r = 0; r < map_->getNumMaps(); ++r)
818 {
819 getMultiVector(r)->randomize(bUseXpetraImplementation);
820 }
821}
822
823
824template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
825void
828{
830}
831
832
833template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
836getMap() const
837{
838 XPETRA_MONITOR("BlockedMultiVector::getMap");
839 return map_;
840}
841
842
843template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
846getBlockedMap() const
847{
848 XPETRA_MONITOR("BlockedMultiVector::getBlockedMap");
849 return map_;
850}
851
852
853template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
856getMultiVector(size_t r) const
857{
858 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r)");
859 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
860 std::out_of_range,
861 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
862 << " partial blocks.");
863 return vv_[ r ];
864}
865
866
867template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
870getMultiVector(size_t r, bool bThyraMode) const
871{
872 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r,bThyraMode)");
873 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
874 std::out_of_range,
875 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
876 << " partial blocks.");
878 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::getMultiVector: inconsistent Thyra mode");
879 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
880 return vv_[ r ];
881}
882
883
884template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
885void
887setMultiVector(size_t r,
889 bool bThyraMode)
890{
891 // The map of the MultiVector should be the same as the stored submap
892 // In thyra mode the vectors should live on the thyra maps
893 // in xpetra mode the should live in the xpetra maps
894 // that should be also ok in the nested case for thyra (if the vectors are distributed accordingly)
895 XPETRA_MONITOR("BlockedMultiVector::setMultiVector");
896 XPETRA_TEST_FOR_EXCEPTION(r >= map_->getNumMaps(),
897 std::out_of_range,
898 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
899 XPETRA_TEST_FOR_EXCEPTION(numVectors_ != v->getNumVectors(),
901 "The BlockedMultiVectors expects " << getNumVectors() << " vectors. The provided partial multivector has "
902 << v->getNumVectors() << " vectors.");
904 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::setMultiVector: inconsistent Thyra mode");
905 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
906 Teuchos::RCP<MultiVector> vv = Teuchos::rcp_const_cast<MultiVector>(v);
907 TEUCHOS_TEST_FOR_EXCEPTION(vv == Teuchos::null, Xpetra::Exceptions::RuntimeError, "Partial vector must not be Teuchos::null");
908 vv_[ r ] = vv;
909}
910
911
912template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
915Merge() const
916{
917 XPETRA_MONITOR("BlockedMultiVector::Merge");
918 using Teuchos::RCP;
919
921 for(size_t r = 0; r < map_->getNumMaps(); ++r)
922 {
923 RCP<MultiVector> vi = getMultiVector(r);
924 RCP<BlockedMultiVector> bvi = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vi);
925 if(bvi.is_null() == true)
926 {
927 this->InsertVector(vi, r, v, map_->getThyraMode());
928 }
929 else
930 {
931 RCP<MultiVector> mvi = bvi->Merge();
932 this->InsertVector(mvi, r, v, map_->getThyraMode());
933 }
934 }
935
936 // TODO plausibility checks
937
938 return v;
939}
940
941
942template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
943void
946{
947 Teuchos::RCP<const MultiVector> rcpRhs = Teuchos::rcpFromRef(rhs);
948 Teuchos::RCP<const BlockedMultiVector> bRhs = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpRhs);
949 if(bRhs == Teuchos::null)
950 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: argument is not a blocked multi vector.");
951
952 map_ = Teuchos::rcp(new BlockedMap(*(bRhs->getBlockedMap())));
953 vv_ = std::vector<Teuchos::RCP<MultiVector>>(map_->getNumMaps());
954 for(size_t r = 0; r < map_->getNumMaps(); ++r)
955 {
956 // extract source vector (is of type MultiVector or BlockedMultiVector)
957 RCP<MultiVector> src = bRhs->getMultiVector(r, map_->getThyraMode());
958
959 // create new (empty) multivector (is of type MultiVector or BlockedMultiVector)
961 map_->getMap(r, bRhs->getBlockedMap()->getThyraMode()), rcpRhs->getNumVectors(), true);
962
963 // check type of source and target vector
964 RCP<BlockedMultiVector> bsrc = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(src);
965 RCP<BlockedMultiVector> bvv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vv);
966 if(bsrc.is_null() == true && bvv.is_null() == true)
967 {
968 *vv = *src; // deep copy
969 }
970 else if(bsrc.is_null() == true && bvv.is_null() == false)
971 {
972 // special case: source vector is a merged MultiVector but target vector is blocked
973 *vv = *src; // deep copy (is this a problem???)
974 }
975 else if(bsrc.is_null() == false && bvv.is_null() == true)
976 {
977 // special case: source vector is blocked but target vector is a merged MultiVector
978 // This is a problem and only works if bsrc has only one block
979 if(bsrc->getBlockedMap()->getNumMaps() > 1)
980 {
981 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: source vector is of type BlockedMultiVector (with more than "
982 "1 blocks) and target is a MultiVector.");
984 }
985 RCP<MultiVector> ssrc = bsrc->getMultiVector(0, map_->getThyraMode());
986 XPETRA_TEST_FOR_EXCEPTION(ssrc.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: cannot extract vector");
987 XPETRA_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<BlockedMultiVector>(ssrc) != Teuchos::null,
989 "BlockedMultiVector::assign: sub block must not be of type BlockedMultiVector.");
990 *vv = *ssrc;
991 }
992 else
993 {
994 // this should work (no exception necessary, i guess)
995 XPETRA_TEST_FOR_EXCEPTION(bsrc->getBlockedMap()->getNumMaps() != bvv->getBlockedMap()->getNumMaps(),
997 "BlockedMultiVector::assign: source and target are BlockedMultiVectors with a different number of submaps.");
998 *vv = *src; // deep copy
999 }
1000 vv_[ r ] = vv;
1001 }
1002 numVectors_ = rcpRhs->getNumVectors();
1003}
1004
1005
1006template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1007void
1010 size_t block,
1012{
1013 ExtractVector(*full, block, *partial);
1014}
1015
1016
1017template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1018void
1021 size_t block,
1023{
1024 ExtractVector(*full, block, *partial);
1025}
1026
1027
1028template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1032 size_t block,
1033 bool bThyraMode) const
1034{
1035 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1036 std::out_of_range,
1037 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1038 << " partial blocks.");
1040 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1041 RCP<const BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(full);
1042 if(bfull.is_null() == true)
1043 {
1044 // standard case: full is not of type BlockedMultiVector
1045 // first extract partial vector from full vector (using xpetra style GIDs)
1046 const RCP<MultiVector> vv =
1047 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1048 // if(bThyraMode == false) {
1049 // ExtractVector(*full, block, *vv);
1050 // return vv;
1051 //} else {
1052 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1053 RCP<MultiVector> rcpNonConstFull = Teuchos::rcp_const_cast<MultiVector>(full);
1054 rcpNonConstFull->replaceMap(map_->getImporter(block)->getSourceMap());
1055 ExtractVector(*rcpNonConstFull, block, *vv);
1056 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1058 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1059 "created using Thyra-style numbered submaps.");
1060 if(bThyraMode == true)
1061 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1062 rcpNonConstFull->replaceMap(oldThyMapFull);
1063 return vv;
1064 //}
1065 }
1066 else
1067 {
1068 // special case: full is of type BlockedMultiVector
1069 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1071 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1072 << bfull->getBlockedMap()->getNumMaps()
1073 << " (number of blocks in BlockedMultiVector)");
1074 return bfull->getMultiVector(block, bThyraMode);
1075 }
1076}
1077
1078
1079template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1083 size_t block,
1084 bool bThyraMode) const
1085{
1086 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1087 std::out_of_range,
1088 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1089 << " partial blocks.");
1091 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1092 RCP<BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(full);
1093 if(bfull.is_null() == true)
1094 {
1095 // standard case: full is not of type BlockedMultiVector
1096 // first extract partial vector from full vector (using xpetra style GIDs)
1097 const RCP<MultiVector> vv =
1098 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1099 // if(bThyraMode == false) {
1100 // ExtractVector(*full, block, *vv);
1101 // return vv;
1102 //} else {
1103 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1104 full->replaceMap(map_->getImporter(block)->getSourceMap());
1105 ExtractVector(*full, block, *vv);
1106 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1108 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1109 "created using Thyra-style numbered submaps.");
1110 if(bThyraMode == true)
1111 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1112 full->replaceMap(oldThyMapFull);
1113 return vv;
1114 //}
1115 }
1116 else
1117 {
1118 // special case: full is of type BlockedMultiVector
1119 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1121 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1122 << bfull->getBlockedMap()->getNumMaps()
1123 << " (number of blocks in BlockedMultiVector)");
1124 return bfull->getMultiVector(block, bThyraMode);
1125 }
1126}
1127
1128
1129template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1130void
1133 size_t block,
1135{
1136 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1137 std::out_of_range,
1138 "ExtractVector: Error, block = " << block << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
1139 << " partial blocks.");
1140 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: maps_[" << block << "] is null");
1141 partial.doImport(full, *(map_->getImporter(block)), Xpetra::INSERT);
1142}
1143
1144
1145template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1146void
1149 size_t block,
1151 bool bThyraMode) const
1152{
1153 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1154 std::out_of_range,
1155 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1156 << " partial blocks.");
1158 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1159 XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1161 "MapExtractor::InsertVector: InsertVector in Thyra-style numbering only possible if MapExtractor has been created "
1162 "using Thyra-style numbered submaps.");
1163 if(bThyraMode)
1164 {
1165 // NOTE: the importer objects in the BlockedMap are always using Xpetra GIDs (or Thyra style Xpetra GIDs)
1166 // The source map corresponds to the full map (in Xpetra GIDs) starting with GIDs from zero. The GIDs are consecutive in Thyra mode
1167 // The target map is the partial map (in the corresponding Xpetra GIDs)
1168
1169 // TODO can we skip the Export call in special cases (i.e. Src = Target map, same length, etc...)
1170
1171 // store original GIDs (could be Thyra GIDs)
1172 RCP<const MultiVector> rcpPartial = Teuchos::rcpFromRef(partial);
1173 RCP<MultiVector> rcpNonConstPartial = Teuchos::rcp_const_cast<MultiVector>(rcpPartial);
1174 RCP<const Map> oldThyMapPartial = rcpNonConstPartial->getMap(); // temporarely store map of partial
1175 RCP<const Map> oldThyMapFull = full.getMap(); // temporarely store map of full
1176
1177 // check whether getMap(block,false) is identical to target map of importer
1178 // XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false)->isSameAs(*(map_->getImporter(block)->getTargetMap()))==false,
1179 // Xpetra::Exceptions::RuntimeError,
1180 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of partial vector are not identical to target Map
1181 // of Importer. This should not be.");
1182
1183 // XPETRA_TEST_FOR_EXCEPTION(full.getMap()->isSameAs(*(map_->getImporter(block)->getSourceMap()))==false,
1184 // Xpetra::Exceptions::RuntimeError,
1185 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of full vector are not identical to source Map of
1186 // Importer. This should not be.");
1187
1188 rcpNonConstPartial->replaceMap(map_->getMap(block, false)); // temporarely switch to xpetra-style map
1189 full.replaceMap(map_->getImporter(block)->getSourceMap()); // temporarely switch to Xpetra GIDs
1190
1191 // do the Export
1192 full.doExport(*rcpNonConstPartial, *(map_->getImporter(block)), Xpetra::INSERT);
1193
1194 // switch back to original maps
1195 full.replaceMap(oldThyMapFull); // reset original map (Thyra GIDs)
1196 rcpNonConstPartial->replaceMap(oldThyMapPartial); // change map back to original map
1197 }
1198 else
1199 {
1200 // Xpetra style numbering
1201 full.doExport(partial, *(map_->getImporter(block)), Xpetra::INSERT);
1202 }
1203}
1204
1205
1206template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1207void
1210 size_t block,
1212 bool bThyraMode) const
1213{
1215 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
1216 if(bfull.is_null() == true)
1217 InsertVector(*partial, block, *full, bThyraMode);
1218 else
1219 {
1220 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1221
1222#if 0
1223 // WCMCLEN - ETI: MultiVector doesn't have a setMultiVector method,
1224 // WCMCLEN - ETI: but BlockedMultiVector does... should this be bfull->...?
1225 full->setMultiVector(block, partial, bThyraMode);
1226#else
1227 throw Xpetra::Exceptions::RuntimeError("MultiVector::setMultiVector() is not implemented.");
1228#endif
1229 }
1230}
1231
1232
1233template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1234void
1237 size_t block,
1239 bool bThyraMode) const
1240{
1242 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
1243 if(bfull.is_null() == true)
1244 InsertVector(*partial, block, *full, bThyraMode);
1245 else
1246 {
1247 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1248 bfull->setMultiVector(block, partial, bThyraMode);
1249 }
1250}
1251
1252
1253} // namespace Xpetra
1254
1255
1256#endif // XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
#define XPETRA_MONITOR(funcName)
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
iterator begin()
iterator end()
bool is_null() const
virtual void meanValue(const Teuchos::ArrayView< Scalar > &) const
Compute mean (average) value of each vector in multi-vector. The outcome of this routine is undefined...
virtual void elementWiseMultiply(Scalar scalarAB, const Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const MultiVector &B, Scalar scalarThis)
Element-wise multiply of a Vector A with a MultiVector B.
Teuchos::RCP< const BlockedMap > map_
blocked map containing the sub block maps (either thyra or xpetra mode)
virtual void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with the given verbosity level to a FancyOStream.
virtual void sumIntoLocalValue(LocalOrdinal, size_t, const Scalar &)
Add value to existing value, using local (row) index.
Teuchos::RCP< const Map > getMap() const
Access function for the underlying Map this DistObject was constructed with.
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const
Const view of the local values in a particular vector of this multivector.
BlockedMultiVector(const Teuchos::RCP< const BlockedMap > &map, size_t NumVectors, bool zeroOut=true)
Constructor.
Teuchos::RCP< MultiVector > getMultiVector(size_t r) const
return partial multivector associated with block row r
void setMultiVector(size_t r, Teuchos::RCP< const MultiVector > v, bool bThyraMode)
set partial multivector associated with block row r
virtual global_size_t getGlobalLength() const
Global number of rows in the multivector.
virtual void doImport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &, const Import &, CombineMode)
Import.
virtual Teuchos::RCP< const Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVector(size_t j) const
Return a Vector which is a const view of column j.
virtual Teuchos::ArrayRCP< Scalar > getDataNonConst(size_t j)
View of the local values in a particular vector of this multivector.
virtual size_t getNumVectors() const
Number of columns in the multivector.
virtual void dot(const MultiVector &, const Teuchos::ArrayView< Scalar > &) const
Compute dot product of each corresponding pair of vectors, dots[i] = this[i].dot(A[i]).
virtual void reciprocal(const MultiVector &)
Put element-wise reciprocal values of input Multi-vector in target, this(i,j) = 1/A(i,...
virtual void sumIntoGlobalValue(GlobalOrdinal, size_t, const Scalar &)
Add value to existing value, using global (row) index.
virtual void replaceLocalValue(LocalOrdinal, size_t, const Scalar &)
Replace value, using local (row) index.
virtual bool isSameSize(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec) const
Local number of rows on the calling process.
Teuchos::RCP< MultiVector > Merge() const
merge BlockedMultiVector blocks to a single MultiVector
virtual Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVectorNonConst(size_t j)
Return a Vector which is a nonconst view of column j.
void ExtractVector(RCP< const MultiVector > &full, size_t block, RCP< MultiVector > &partial) const
void InsertVector(const MultiVector &partial, size_t block, MultiVector &full, bool bThyraMode=false) const
virtual void randomize(bool bUseXpetraImplementation=false)
Set multi-vector values to random numbers.
virtual void replaceGlobalValue(GlobalOrdinal, size_t, const Scalar &)
Replace value, using global (row) index.
virtual void norm2(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void doExport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &, const Import &, CombineMode)
Export.
virtual void setSeed(unsigned int seed)
Set seed for Random function.
size_t numVectors_
number of vectors (columns in multi vector)
virtual void assign(const MultiVector &rhs)
Implementation of the assignment operator (operator=); does a deep copy.
virtual std::string description() const
A simple one-line description of this object.
virtual void scale(const Scalar &alpha)
Scale the current values of a multi-vector, this = alpha*this.
virtual void normInf(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const
Compute Inf-norm of each vector in multi-vector.
Teuchos::RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > getBlockedMap() const
Access function for the underlying Map this DistObject was constructed with.
std::vector< Teuchos::RCP< MultiVector > > vv_
array containing RCPs of the partial vectors
virtual void update(const Scalar &alpha, const MultiVector &A, const Scalar &beta)
Update multi-vector values with scaled values of A, this = beta*this + alpha*A.
virtual void Xpetra_randomize()
Set multi-vector values to random numbers. XPetra implementation.
virtual void abs(const MultiVector &)
Put element-wise absolute values of input Multi-vector in target: A = abs(this).
virtual void norm1(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const
Compute 1-norm of each vector in multi-vector.
BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & operator=(const MultiVector &rhs)
Assignment operator: Does a deep copy.
virtual void multiply(Teuchos::ETransp, Teuchos::ETransp, const Scalar &, const MultiVector &, const MultiVector &, const Scalar &)
Matrix-matrix multiplication: this = beta*this + alpha*op(A)*op(B).
virtual size_t getLocalLength() const
Local number of rows on the calling process.
virtual void replaceMap(const RCP< const Map > &map)
virtual void doExport(const DistObject< Packet, LocalOrdinal, GlobalOrdinal, Node > &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)=0
Export data into this object using an Export object ("forward mode").
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
virtual void doImport(const DistObject< Packet, LocalOrdinal, GlobalOrdinal, Node > &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM)=0
Import data into this object using an Import object ("forward mode").
Exception throws to report errors in the internal logical of the program.
Factory for any type of Xpetra::MultiVector and its derived classes.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
virtual void replaceMap(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map)=0
virtual void Xpetra_randomize()
Set multi-vector values to random numbers. XPetra implementation.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
size_t global_size_t
Global size_t object.
CombineMode
Xpetra::Combine Mode enumerable type.