Xpetra Version of the Day
Loading...
Searching...
No Matches
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()->getLocalNumElements() != v->getMap()->getLocalNumElements(),
96 "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has "
97 << bmap->getMap()->getLocalNumElements() << " local elements. The vector has " << v->getMap()->getLocalNumElements()
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()->getLocalNumElements() != v->getMap()->getLocalNumElements(), Xpetra::Exceptions::RuntimeError,
105 // "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " <<
106 // bmap->getFullMap()->getLocalNumElements() << " local elements. The vector has " << v->getMap()->getLocalNumElements() << ".");
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.");
326
327
328template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329void
331scale(const Scalar& alpha)
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);
340 }
341}
342
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)
352 if(getMultiVector(r) != Teuchos::null)
353 {
354 getMultiVector(r)->scale(alpha);
356 }
357}
358
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++)
383
384 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getLocalNumElements() != bA->getMultiVector(r)->getMap()->getLocalNumElements(),
386 "BlockedMultiVector::update: in subvector "
387 << r << ": Cannot add a vector of (local) length " << bA->getMultiVector(r)->getMap()->getLocalNumElements()
388 << " to the existing vector with " << getMultiVector(r)->getMap()->getLocalNumElements() << " 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()->getLocalNumElements() != part->getMap()->getLocalNumElements(),
450 "BlockedMultiVector::update: in subvector "
451 << r << ": Cannot add a vector of (local) length " << part->getMap()->getLocalNumElements()
452 << " to the existing vector with " << getMultiVector(r)->getMap()->getLocalNumElements() << " 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()->getLocalNumElements() != B.getMap()->getLocalNumElements(),
621 "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getLocalNumElements() << " elements, B has "
622 << B.getMap()->getLocalNumElements() << ".");
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()->getLocalNumElements();
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
827randomize(const Scalar& minVal, const Scalar& maxVal, bool bUseXpetraImplementation)
828{
829 for(size_t r = 0; r < map_->getNumMaps(); ++r)
830 {
831 getMultiVector(r)->randomize(minVal, maxVal, bUseXpetraImplementation);
832 }
833}
834
835
836template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
837void
840{
842}
843
844
845template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
846void
848Xpetra_randomize(const Scalar& minVal, const Scalar& maxVal)
849{
851}
852
853
854template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
857getMap() const
858{
859 XPETRA_MONITOR("BlockedMultiVector::getMap");
860 return map_;
861}
862
863
864template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
867getBlockedMap() const
868{
869 XPETRA_MONITOR("BlockedMultiVector::getBlockedMap");
870 return map_;
871}
872
873
874template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
877getMultiVector(size_t r) const
878{
879 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r)");
880 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
881 std::out_of_range,
882 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
883 << " partial blocks.");
884 return vv_[ r ];
885}
886
887
888template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
891getMultiVector(size_t r, bool bThyraMode) const
892{
893 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r,bThyraMode)");
894 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
895 std::out_of_range,
896 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
897 << " partial blocks.");
899 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::getMultiVector: inconsistent Thyra mode");
900 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
901 return vv_[ r ];
902}
903
904
905template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
906void
908setMultiVector(size_t r,
910 bool bThyraMode)
911{
912 // The map of the MultiVector should be the same as the stored submap
913 // In thyra mode the vectors should live on the thyra maps
914 // in xpetra mode the should live in the xpetra maps
915 // that should be also ok in the nested case for thyra (if the vectors are distributed accordingly)
916 XPETRA_MONITOR("BlockedMultiVector::setMultiVector");
917 XPETRA_TEST_FOR_EXCEPTION(r >= map_->getNumMaps(),
918 std::out_of_range,
919 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
920 XPETRA_TEST_FOR_EXCEPTION(numVectors_ != v->getNumVectors(),
922 "The BlockedMultiVectors expects " << getNumVectors() << " vectors. The provided partial multivector has "
923 << v->getNumVectors() << " vectors.");
925 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::setMultiVector: inconsistent Thyra mode");
926 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
927 Teuchos::RCP<MultiVector> vv = Teuchos::rcp_const_cast<MultiVector>(v);
928 TEUCHOS_TEST_FOR_EXCEPTION(vv == Teuchos::null, Xpetra::Exceptions::RuntimeError, "Partial vector must not be Teuchos::null");
929 vv_[ r ] = vv;
930}
931
932
933template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
936Merge() const
937{
938 XPETRA_MONITOR("BlockedMultiVector::Merge");
939 using Teuchos::RCP;
940
942 for(size_t r = 0; r < map_->getNumMaps(); ++r)
943 {
944 RCP<MultiVector> vi = getMultiVector(r);
945 RCP<BlockedMultiVector> bvi = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vi);
946 if(bvi.is_null() == true)
947 {
948 this->InsertVector(vi, r, v, map_->getThyraMode());
949 }
950 else
951 {
952 RCP<MultiVector> mvi = bvi->Merge();
953 this->InsertVector(mvi, r, v, map_->getThyraMode());
954 }
955 }
956
957 // TODO plausibility checks
958
959 return v;
960}
961
962
963template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
964void
967{
968 Teuchos::RCP<const MultiVector> rcpRhs = Teuchos::rcpFromRef(rhs);
969 Teuchos::RCP<const BlockedMultiVector> bRhs = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpRhs);
970 if(bRhs == Teuchos::null)
971 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: argument is not a blocked multi vector.");
972
973 map_ = Teuchos::rcp(new BlockedMap(*(bRhs->getBlockedMap())));
974 vv_ = std::vector<Teuchos::RCP<MultiVector>>(map_->getNumMaps());
975 for(size_t r = 0; r < map_->getNumMaps(); ++r)
976 {
977 // extract source vector (is of type MultiVector or BlockedMultiVector)
978 RCP<MultiVector> src = bRhs->getMultiVector(r, map_->getThyraMode());
979
980 // create new (empty) multivector (is of type MultiVector or BlockedMultiVector)
982 map_->getMap(r, bRhs->getBlockedMap()->getThyraMode()), rcpRhs->getNumVectors(), true);
983
984 // check type of source and target vector
985 RCP<BlockedMultiVector> bsrc = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(src);
986 RCP<BlockedMultiVector> bvv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vv);
987 if(bsrc.is_null() == true && bvv.is_null() == true)
988 {
989 *vv = *src; // deep copy
990 }
991 else if(bsrc.is_null() == true && bvv.is_null() == false)
992 {
993 // special case: source vector is a merged MultiVector but target vector is blocked
994 *vv = *src; // deep copy (is this a problem???)
995 }
996 else if(bsrc.is_null() == false && bvv.is_null() == true)
997 {
998 // special case: source vector is blocked but target vector is a merged MultiVector
999 // This is a problem and only works if bsrc has only one block
1000 if(bsrc->getBlockedMap()->getNumMaps() > 1)
1001 {
1002 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: source vector is of type BlockedMultiVector (with more than "
1003 "1 blocks) and target is a MultiVector.");
1005 }
1006 RCP<MultiVector> ssrc = bsrc->getMultiVector(0, map_->getThyraMode());
1007 XPETRA_TEST_FOR_EXCEPTION(ssrc.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: cannot extract vector");
1008 XPETRA_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<BlockedMultiVector>(ssrc) != Teuchos::null,
1010 "BlockedMultiVector::assign: sub block must not be of type BlockedMultiVector.");
1011 *vv = *ssrc;
1012 }
1013 else
1014 {
1015 // this should work (no exception necessary, i guess)
1016 XPETRA_TEST_FOR_EXCEPTION(bsrc->getBlockedMap()->getNumMaps() != bvv->getBlockedMap()->getNumMaps(),
1018 "BlockedMultiVector::assign: source and target are BlockedMultiVectors with a different number of submaps.");
1019 *vv = *src; // deep copy
1020 }
1021 vv_[ r ] = vv;
1022 }
1023 numVectors_ = rcpRhs->getNumVectors();
1024}
1025
1026
1027template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1028void
1031 size_t block,
1033{
1034 ExtractVector(*full, block, *partial);
1035}
1036
1037
1038template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1039void
1042 size_t block,
1044{
1045 ExtractVector(*full, block, *partial);
1046}
1047
1048
1049template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1053 size_t block,
1054 bool bThyraMode) const
1055{
1056 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1057 std::out_of_range,
1058 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1059 << " partial blocks.");
1061 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1062 RCP<const BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(full);
1063 if(bfull.is_null() == true)
1064 {
1065 // standard case: full is not of type BlockedMultiVector
1066 // first extract partial vector from full vector (using xpetra style GIDs)
1067 const RCP<MultiVector> vv =
1068 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1069 // if(bThyraMode == false) {
1070 // ExtractVector(*full, block, *vv);
1071 // return vv;
1072 //} else {
1073 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1074 RCP<MultiVector> rcpNonConstFull = Teuchos::rcp_const_cast<MultiVector>(full);
1075 rcpNonConstFull->replaceMap(map_->getImporter(block)->getSourceMap());
1076 ExtractVector(*rcpNonConstFull, block, *vv);
1077 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1079 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1080 "created using Thyra-style numbered submaps.");
1081 if(bThyraMode == true)
1082 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1083 rcpNonConstFull->replaceMap(oldThyMapFull);
1084 return vv;
1085 //}
1086 }
1087 else
1088 {
1089 // special case: full is of type BlockedMultiVector
1090 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1092 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1093 << bfull->getBlockedMap()->getNumMaps()
1094 << " (number of blocks in BlockedMultiVector)");
1095 return bfull->getMultiVector(block, bThyraMode);
1096 }
1097}
1098
1099
1100template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1104 size_t block,
1105 bool bThyraMode) const
1106{
1107 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1108 std::out_of_range,
1109 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1110 << " partial blocks.");
1112 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1113 RCP<BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(full);
1114 if(bfull.is_null() == true)
1115 {
1116 // standard case: full is not of type BlockedMultiVector
1117 // first extract partial vector from full vector (using xpetra style GIDs)
1118 const RCP<MultiVector> vv =
1119 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1120 // if(bThyraMode == false) {
1121 // ExtractVector(*full, block, *vv);
1122 // return vv;
1123 //} else {
1124 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1125 full->replaceMap(map_->getImporter(block)->getSourceMap());
1126 ExtractVector(*full, block, *vv);
1127 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1129 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1130 "created using Thyra-style numbered submaps.");
1131 if(bThyraMode == true)
1132 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1133 full->replaceMap(oldThyMapFull);
1134 return vv;
1135 //}
1136 }
1137 else
1138 {
1139 // special case: full is of type BlockedMultiVector
1140 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1142 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1143 << bfull->getBlockedMap()->getNumMaps()
1144 << " (number of blocks in BlockedMultiVector)");
1145 return bfull->getMultiVector(block, bThyraMode);
1146 }
1147}
1148
1149
1150template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1151void
1154 size_t block,
1156{
1157 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1158 std::out_of_range,
1159 "ExtractVector: Error, block = " << block << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
1160 << " partial blocks.");
1161 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: maps_[" << block << "] is null");
1162 partial.doImport(full, *(map_->getImporter(block)), Xpetra::INSERT);
1163}
1164
1165
1166template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1167void
1170 size_t block,
1172 bool bThyraMode) const
1173{
1174 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1175 std::out_of_range,
1176 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1177 << " partial blocks.");
1179 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1180 XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1182 "MapExtractor::InsertVector: InsertVector in Thyra-style numbering only possible if MapExtractor has been created "
1183 "using Thyra-style numbered submaps.");
1184 if(bThyraMode)
1185 {
1186 // NOTE: the importer objects in the BlockedMap are always using Xpetra GIDs (or Thyra style Xpetra GIDs)
1187 // The source map corresponds to the full map (in Xpetra GIDs) starting with GIDs from zero. The GIDs are consecutive in Thyra mode
1188 // The target map is the partial map (in the corresponding Xpetra GIDs)
1189
1190 // TODO can we skip the Export call in special cases (i.e. Src = Target map, same length, etc...)
1191
1192 // store original GIDs (could be Thyra GIDs)
1193 RCP<const MultiVector> rcpPartial = Teuchos::rcpFromRef(partial);
1194 RCP<MultiVector> rcpNonConstPartial = Teuchos::rcp_const_cast<MultiVector>(rcpPartial);
1195 RCP<const Map> oldThyMapPartial = rcpNonConstPartial->getMap(); // temporarely store map of partial
1196 RCP<const Map> oldThyMapFull = full.getMap(); // temporarely store map of full
1197
1198 // check whether getMap(block,false) is identical to target map of importer
1199 // XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false)->isSameAs(*(map_->getImporter(block)->getTargetMap()))==false,
1200 // Xpetra::Exceptions::RuntimeError,
1201 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of partial vector are not identical to target Map
1202 // of Importer. This should not be.");
1203
1204 // XPETRA_TEST_FOR_EXCEPTION(full.getMap()->isSameAs(*(map_->getImporter(block)->getSourceMap()))==false,
1205 // Xpetra::Exceptions::RuntimeError,
1206 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of full vector are not identical to source Map of
1207 // Importer. This should not be.");
1208
1209 rcpNonConstPartial->replaceMap(map_->getMap(block, false)); // temporarely switch to xpetra-style map
1210 full.replaceMap(map_->getImporter(block)->getSourceMap()); // temporarely switch to Xpetra GIDs
1211
1212 // do the Export
1213 full.doExport(*rcpNonConstPartial, *(map_->getImporter(block)), Xpetra::INSERT);
1214
1215 // switch back to original maps
1216 full.replaceMap(oldThyMapFull); // reset original map (Thyra GIDs)
1217 rcpNonConstPartial->replaceMap(oldThyMapPartial); // change map back to original map
1218 }
1219 else
1220 {
1221 // Xpetra style numbering
1222 full.doExport(partial, *(map_->getImporter(block)), Xpetra::INSERT);
1223 }
1224}
1225
1226
1227template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1228void
1231 size_t block,
1233 bool bThyraMode) const
1234{
1236 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
1237 if(bfull.is_null() == true)
1238 InsertVector(*partial, block, *full, bThyraMode);
1239 else
1240 {
1241 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1242
1243#if 0
1244 // WCMCLEN - ETI: MultiVector doesn't have a setMultiVector method,
1245 // WCMCLEN - ETI: but BlockedMultiVector does... should this be bfull->...?
1246 full->setMultiVector(block, partial, bThyraMode);
1247#else
1248 throw Xpetra::Exceptions::RuntimeError("MultiVector::setMultiVector() is not implemented.");
1249#endif
1250 }
1251}
1252
1253
1254template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1255void
1258 size_t block,
1260 bool bThyraMode) const
1261{
1263 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
1264 if(bfull.is_null() == true)
1265 InsertVector(*partial, block, *full, bThyraMode);
1266 else
1267 {
1268 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1269 bfull->setMultiVector(block, partial, bThyraMode);
1270 }
1271}
1272
1273
1274} // namespace Xpetra
1275
1276
1277#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.