Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_BlockedMultiVector_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Xpetra: A linear algebra interface package
4//
5// Copyright 2012 NTESS and the Xpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
11#define XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
12
13#include "Xpetra_BlockedMultiVector_decl.hpp"
14
15#include "Xpetra_MultiVectorFactory.hpp"
16#include "Xpetra_BlockedVector.hpp"
17#include "Xpetra_MapExtractor.hpp"
18
19namespace Xpetra {
20
21template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
22BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
23 BlockedMultiVector(const Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>>& map,
24 size_t NumVectors,
25 bool zeroOut)
26 : map_(map) {
27 numVectors_ = NumVectors;
28
29 vv_.reserve(map->getNumMaps());
30
31 // add CrsMatrix objects in row,column order
32 for (size_t r = 0; r < map->getNumMaps(); ++r)
33 vv_.push_back(
34 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map->getMap(r, map_->getThyraMode()), NumVectors, zeroOut));
35}
36
48template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
49BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
50 BlockedMultiVector(Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>> bmap,
51 Teuchos::RCP<const MultiVector> v) {
52 XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getLocalNumElements() != v->getMap()->getLocalNumElements(),
53 Xpetra::Exceptions::RuntimeError,
54 "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has "
55 << bmap->getMap()->getLocalNumElements() << " local elements. The vector has " << v->getMap()->getLocalNumElements()
56 << ".");
57 XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(),
58 Xpetra::Exceptions::RuntimeError,
59 "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has "
60 << bmap->getMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements()
61 << ".");
62 // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getLocalNumElements() != v->getMap()->getLocalNumElements(), Xpetra::Exceptions::RuntimeError,
63 // "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " <<
64 // bmap->getFullMap()->getLocalNumElements() << " local elements. The vector has " << v->getMap()->getLocalNumElements() << ".");
65 // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError,
66 // "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has " <<
67 // bmap->getFullMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements() << ".");
68
69 numVectors_ = v->getNumVectors();
70
71 map_ = bmap;
72
73 // Extract vector
74 vv_.reserve(bmap->getNumMaps());
75
76 // add CrsMatrix objects in row,column order
77 for (size_t r = 0; r < bmap->getNumMaps(); ++r)
78 vv_.push_back(this->ExtractVector(v, r, bmap->getThyraMode()));
79}
80
92template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
94 BlockedMultiVector(Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>> mapExtractor,
95 Teuchos::RCP<const MultiVector> v) {
96 numVectors_ = v->getNumVectors();
97
98 // create blocked map out of MapExtractor
99 std::vector<RCP<const Map>> maps;
100 maps.reserve(mapExtractor->NumMaps());
101 for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
102 maps.push_back(mapExtractor->getMap(r, mapExtractor->getThyraMode()));
103 map_ = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapExtractor->getFullMap(), maps, mapExtractor->getThyraMode()));
104
105 // Extract vector
106 vv_.reserve(mapExtractor->NumMaps());
107
108 // add CrsMatrix objects in row,column order
109 for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
110 vv_.push_back(this->ExtractVector(v, r, mapExtractor->getThyraMode()));
111}
112
113template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
115 BlockedMultiVector(const Teuchos::RCP<const BlockedMap>& map,
116 std::vector<Teuchos::RCP<MultiVector>>& vin) {
117 numVectors_ = vin[0]->getNumVectors();
118 map_ = map;
119 vv_.resize(vin.size());
120 for (size_t i = 0; i < vv_.size(); i++)
121 vv_[i] = vin[i];
122}
123
124template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
126 ~BlockedMultiVector() {
127 for (size_t r = 0; r < vv_.size(); ++r) {
128 vv_[r] = Teuchos::null;
129 }
130
131 map_ = Teuchos::null;
132 numVectors_ = 0;
133}
134
135template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
137BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
138operator=(const MultiVector& rhs) {
139 assign(rhs); // dispatch to protected virtual method
140 return *this;
141}
142
143template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
144void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
145 replaceGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */) {
146 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceGlobalValue: Not (yet) supported by BlockedMultiVector.");
147}
148
149template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
150void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
151 sumIntoGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */) {
152 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoGlobalValue: Not (yet) supported by BlockedMultiVector.");
153}
154
155template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
156void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
157 replaceLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */) {
158 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceLocalValue: Not supported by BlockedMultiVector.");
159}
160
161template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
162void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
163 sumIntoLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */) {
164 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoLocalValue:Not (yet) supported by BlockedMultiVector.");
165}
166
167template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
168void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
169 putScalar(const Scalar& value) {
170 XPETRA_MONITOR("BlockedMultiVector::putScalar");
171 for (size_t r = 0; r < map_->getNumMaps(); r++) {
172 getMultiVector(r)->putScalar(value);
173 }
174}
175
176template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
177Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
178BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
179 getVector(size_t j) const {
180 RCP<Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
181 Teuchos::rcp(new Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(this->getBlockedMap(), false));
182
183 for (size_t r = 0; r < getBlockedMap()->getNumMaps(); r++) {
184 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> subvec =
185 this->getMultiVector(r, this->getBlockedMap()->getThyraMode())->getVector(j);
186 ret->setMultiVector(r, subvec, this->getBlockedMap()->getThyraMode());
187 }
188 return ret;
189}
190
191template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
192Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
193BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
194 getVectorNonConst(size_t j) {
195 RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
196 Teuchos::rcp_const_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(getVector(j));
197 return ret;
198}
199
200template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
201Teuchos::ArrayRCP<const Scalar>
202BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
203 getData(size_t j) const {
204 if (map_->getNumMaps() == 1) {
205 return vv_[0]->getData(j);
206 }
207 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getData: Not (yet) supported by BlockedMultiVector.");
208 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
209}
210
211template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
212Teuchos::ArrayRCP<Scalar>
213BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
214 getDataNonConst(size_t j) {
215 if (map_->getNumMaps() == 1) {
216 return vv_[0]->getDataNonConst(j);
217 }
218 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getDataNonConst: Not (yet) supported by BlockedMultiVector.");
219 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
220}
221
222template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
223void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
224 dot(const MultiVector& /* A */, const Teuchos::ArrayView<Scalar>& /* dots */) const {
225 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::dot: Not (yet) supported by BlockedMultiVector.");
226}
227
228template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
230 abs(const MultiVector& /* A */) {
231 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::abs: Not (yet) supported by BlockedMultiVector.");
232}
233
234template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
235void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
236 reciprocal(const MultiVector& /* A */) {
237 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::reciprocal: Not (yet) supported by BlockedMultiVector.");
238}
239
240template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
241void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
242 scale(const Scalar& alpha) {
243 XPETRA_MONITOR("BlockedMultiVector::scale (Scalar)");
244 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
245 if (getMultiVector(r) != Teuchos::null) {
246 getMultiVector(r)->scale(alpha);
247 }
248 }
249}
250
251template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
252void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
253 scale(Teuchos::ArrayView<const Scalar> alpha) {
254 XPETRA_MONITOR("BlockedMultiVector::scale (Array)");
255 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
256 if (getMultiVector(r) != Teuchos::null) {
257 getMultiVector(r)->scale(alpha);
258 }
259 }
260}
261
262template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
263void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
264 update(const Scalar& alpha, const MultiVector& A, const Scalar& beta) {
265 XPETRA_MONITOR("BlockedMultiVector::update");
266 Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
267 Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
268 TEUCHOS_TEST_FOR_EXCEPTION(numVectors_ != rcpA->getNumVectors(),
269 Xpetra::Exceptions::RuntimeError,
270 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector).");
271 if (bA != Teuchos::null) {
272 // A is a BlockedMultiVector (and compatible with this)
273 // Call update recursively on all sub vectors
274 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(),
275 Xpetra::Exceptions::RuntimeError,
276 "BlockedMultiVector::update: update with incompatible vector (different thyra mode).");
277 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(),
278 Xpetra::Exceptions::RuntimeError,
279 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors).");
280 for (size_t r = 0; r < map_->getNumMaps(); r++) {
281 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getLocalNumElements() != bA->getMultiVector(r)->getMap()->getLocalNumElements(),
282 Xpetra::Exceptions::RuntimeError,
283 "BlockedMultiVector::update: in subvector "
284 << r << ": Cannot add a vector of (local) length " << bA->getMultiVector(r)->getMap()->getLocalNumElements()
285 << " to the existing vector with " << getMultiVector(r)->getMap()->getLocalNumElements() << " entries.");
286 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != bA->getMultiVector(r)->getMap()->getGlobalNumElements(),
287 Xpetra::Exceptions::RuntimeError,
288 "BlockedMultiVector::update: in subvector "
289 << r << ": Cannot add a vector of length " << bA->getMultiVector(r)->getMap()->getGlobalNumElements()
290 << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
291
292 // TAW: 12/6 We basically want to do something like:
293 // getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta);
294 // But if the left hand side is a standard MultiVector and bA->getMultiVector(r) is
295 // a blocked MultiVector this will not work, as the update implementation of the standard
296 // multivector cannot deal with blocked multivectors.
297 // The only use case where this could happen is, if the rhs vector is a single block multivector
298 Teuchos::RCP<MultiVector> lmv = getMultiVector(r);
299 Teuchos::RCP<MultiVector> rmv = bA->getMultiVector(r);
300
301 // check whether lmv/rmv are blocked or not
302 Teuchos::RCP<BlockedMultiVector> blmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(lmv);
303 Teuchos::RCP<BlockedMultiVector> brmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rmv);
304
305 if (blmv.is_null() == true && brmv.is_null() == false) {
306 // special case: lmv is standard MultiVector but rmv is BlockedMultiVector
307 TEUCHOS_TEST_FOR_EXCEPTION(brmv->getBlockedMap()->getNumMaps() > 1,
308 Xpetra::Exceptions::RuntimeError,
309 "BlockedMultiVector::update: Standard MultiVector object does not accept BlockedMultVector object as "
310 "parameter in update call.");
311 lmv->update(alpha, *(brmv->getMultiVector(0)), beta);
312 } else
313 lmv->update(alpha, *rmv, beta);
314 }
315 } else {
316 // A is a MultiVector
317 // If this is a BlockedMultiVector with only one sub-vector of same length we can just update
318 // Otherwise, A is not compatible with this as BlockedMultiVector and we have to extract the vector data from A
319
320 if (getBlockedMap()->getNumMaps() == 1) {
321 // Actually call "update" on the underlying MultiVector (Epetra or Tpetra).
322 // The maps have to be compatible.
323 XPETRA_TEST_FOR_EXCEPTION(
324 getMultiVector(0)->getMap()->isSameAs(*(rcpA->getMap())) == false,
325 Xpetra::Exceptions::RuntimeError,
326 "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor).");
327 getMultiVector(0)->update(alpha, *rcpA, beta);
328 } else {
329 // general case: A has to be splitted and subparts have to be extracted and stored in this BlockedMultiVector
330 // XPETRA_TEST_FOR_EXCEPTION(map_->getFullMap()->isSameAs(*(rcpA->getMap()))==false, Xpetra::Exceptions::RuntimeError,
331 // "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor). - Note:
332 // This test might be too strict and can probably be relaxed!");
333 for (size_t r = 0; r < map_->getNumMaps(); r++) {
334 // Call "update" on the subvector. Note, that getMultiVector(r) could return another BlockedMultiVector.
335 // That is, in Thyra mode the maps could differ (local Xpetra versus Thyra style gids)
336 Teuchos::RCP<const MultiVector> part = this->ExtractVector(rcpA, r, map_->getThyraMode());
337 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getLocalNumElements() != part->getMap()->getLocalNumElements(),
338 Xpetra::Exceptions::RuntimeError,
339 "BlockedMultiVector::update: in subvector "
340 << r << ": Cannot add a vector of (local) length " << part->getMap()->getLocalNumElements()
341 << " to the existing vector with " << getMultiVector(r)->getMap()->getLocalNumElements() << " entries.");
342 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != part->getMap()->getGlobalNumElements(),
343 Xpetra::Exceptions::RuntimeError,
344 "BlockedMultiVector::update: in subvector "
345 << r << ": Cannot add a vector of length " << part->getMap()->getGlobalNumElements()
346 << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
347 getMultiVector(r)->update(alpha, *part, beta);
348 }
349 }
350 }
351}
352
353template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
354void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
355 update(const Scalar& alpha, const MultiVector& A, const Scalar& beta, const MultiVector& B, const Scalar& gamma) {
356 XPETRA_MONITOR("BlockedMultiVector::update2");
357 Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
358 Teuchos::RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
359 Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
360 Teuchos::RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
361 if (bA != Teuchos::null && bB != Teuchos::null) {
362 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(),
363 Xpetra::Exceptions::RuntimeError,
364 "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector A).");
365 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(),
366 Xpetra::Exceptions::RuntimeError,
367 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector A).");
368 TEUCHOS_TEST_FOR_EXCEPTION(
369 numVectors_ != bA->getNumVectors(),
370 Xpetra::Exceptions::RuntimeError,
371 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector A).");
372 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bB->getBlockedMap()->getThyraMode(),
373 Xpetra::Exceptions::RuntimeError,
374 "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector B).");
375 TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bB->getBlockedMap()->getNumMaps(),
376 Xpetra::Exceptions::RuntimeError,
377 "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector B).");
378 TEUCHOS_TEST_FOR_EXCEPTION(
379 numVectors_ != bB->getNumVectors(),
380 Xpetra::Exceptions::RuntimeError,
381 "BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector B).");
382
383 for (size_t r = 0; r < map_->getNumMaps(); r++) {
384 XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->isSameAs(*(bA->getMultiVector(r)->getMap())) == false,
385 Xpetra::Exceptions::RuntimeError,
386 "BlockedMultiVector::update: update with incompatible vector (different maps in partial vector " << r << ").");
387 getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta, *(bB->getMultiVector(r)), gamma);
388 }
389 return;
390 }
391 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::update: only supports update with other BlockedMultiVector.");
392}
393
394template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
395void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
396 norm1(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const {
397 XPETRA_MONITOR("BlockedMultiVector::norm1");
398 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
399 Array<Magnitude> temp_norms(getNumVectors());
400 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
401 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
402 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
403 if (getMultiVector(r) != Teuchos::null) {
404 getMultiVector(r)->norm1(temp_norms);
405 for (size_t c = 0; c < getNumVectors(); ++c)
406 norms[c] += temp_norms[c];
407 }
408 }
409}
410
411template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
412void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
413 norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const {
414 XPETRA_MONITOR("BlockedMultiVector::norm2");
415 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
416 Array<Magnitude> results(getNumVectors());
417 Array<Magnitude> temp_norms(getNumVectors());
418 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
419 std::fill(results.begin(), results.end(), ScalarTraits<Magnitude>::zero());
420 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
421 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
422 if (getMultiVector(r) != Teuchos::null) {
423 getMultiVector(r)->norm2(temp_norms);
424 for (size_t c = 0; c < getNumVectors(); ++c)
425 results[c] += temp_norms[c] * temp_norms[c];
426 }
427 }
428 for (size_t c = 0; c < getNumVectors(); ++c)
429 norms[c] = Teuchos::ScalarTraits<Magnitude>::squareroot(results[c]);
430}
431
432template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
433void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
434 normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const {
435 XPETRA_MONITOR("BlockedMultiVector::normInf");
436 typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
437 Array<Magnitude> temp_norms(getNumVectors());
438 std::fill(norms.begin(), norms.end(), ScalarTraits<Magnitude>::zero());
439 std::fill(temp_norms.begin(), temp_norms.end(), ScalarTraits<Magnitude>::zero());
440 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
441 if (getMultiVector(r) != Teuchos::null) {
442 getMultiVector(r)->normInf(temp_norms);
443 for (size_t c = 0; c < getNumVectors(); ++c)
444 norms[c] = std::max(norms[c], temp_norms[c]);
445 }
446 }
447}
448
449template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
450void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
451 meanValue(const Teuchos::ArrayView<Scalar>& /* means */) const {
452 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::meanValue: Not (yet) supported by BlockedMultiVector.");
453}
454
455template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
456void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
457 multiply(Teuchos::ETransp /* transA */,
458 Teuchos::ETransp /* transB */,
459 const Scalar& /* alpha */,
460 const MultiVector& /* A */,
461 const MultiVector& /* B */,
462 const Scalar& /* beta */) {
463 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::multiply: Not (yet) supported by BlockedMultiVector.");
464}
465
466template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
467void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
468 elementWiseMultiply(Scalar scalarAB,
469 const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
470 const MultiVector& B,
471 Scalar scalarThis) {
472 XPETRA_MONITOR("BlockedMultiVector::elementWiseMultiply");
473 XPETRA_TEST_FOR_EXCEPTION(B.getMap()->isSameAs(*(this->getMap())) == false,
474 Xpetra::Exceptions::RuntimeError,
475 "BlockedMultiVector::elementWiseMultipy: B must have same blocked map than this.");
476 // XPETRA_TEST_FOR_EXCEPTION(A.getMap()->isSameAs(*(this->getMap()))==false, Xpetra::Exceptions::RuntimeError,
477 // "BlockedMultiVector::elementWiseMultipy: A must have same blocked map than this.");
478 TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getLocalNumElements() != B.getMap()->getLocalNumElements(),
479 Xpetra::Exceptions::RuntimeError,
480 "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getLocalNumElements() << " elements, B has "
481 << B.getMap()->getLocalNumElements() << ".");
482 TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getGlobalNumElements() != B.getMap()->getGlobalNumElements(),
483 Xpetra::Exceptions::RuntimeError,
484 "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getGlobalNumElements() << " elements, B has "
485 << B.getMap()->getGlobalNumElements() << ".");
486
487 RCP<const BlockedMap> bmap = getBlockedMap();
488 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rcpA = Teuchos::rcpFromRef(A);
489 RCP<const Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bA =
490 Teuchos::rcp_dynamic_cast<const Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(rcpA);
491
492 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
493 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
494 TEUCHOS_TEST_FOR_EXCEPTION(
495 bB.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: B must be a BlockedMultiVector.");
496
497 // RCP<Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > me = Teuchos::rcp(new
498 // Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(bmap));
499
500 for (size_t m = 0; m < bmap->getNumMaps(); m++) {
501 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partA = bA->getMultiVector(m, bmap->getThyraMode())->getVector(0);
502 RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partB = bB->getMultiVector(m, bmap->getThyraMode());
503 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thisPart = this->getMultiVector(m, bmap->getThyraMode());
504
505 thisPart->elementWiseMultiply(scalarAB, *partA, *partB, scalarThis);
506 }
507}
508
509template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
510size_t
515
516template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
517size_t
518BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
519 getLocalLength() const {
520 XPETRA_MONITOR("BlockedMultiVector::getLocalLength()");
521 return map_->getFullMap()->getLocalNumElements();
522}
523
524template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
525global_size_t
526BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
527 getGlobalLength() const {
528 XPETRA_MONITOR("BlockedMultiVector::getGlobalLength()");
529 return map_->getFullMap()->getGlobalNumElements();
530}
531
532template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
533bool BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
534 isSameSize(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) const {
535 const BlockedMultiVector* Vb = dynamic_cast<const BlockedMultiVector*>(&vec);
536 if (!Vb)
537 return false;
538 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
539 RCP<const MultiVector> a = getMultiVector(r);
540 RCP<const MultiVector> b = Vb->getMultiVector(r);
541 if ((a == Teuchos::null && b != Teuchos::null) || (a != Teuchos::null && b == Teuchos::null))
542 return false;
543 if (a != Teuchos::null && b != Teuchos::null && !a->isSameSize(*b))
544 return false;
545 }
546 return true;
547}
548
549template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
550std::string
551BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
552 description() const {
553 return std::string("BlockedMultiVector");
554}
555
556template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
557void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
558 describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const {
559 out << description() << std::endl;
560 for (size_t r = 0; r < map_->getNumMaps(); r++)
561 getMultiVector(r)->describe(out, verbLevel);
562}
563
564template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
565void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
566 replaceMap(const RCP<const Map>& map) {
567 XPETRA_MONITOR("BlockedMultiVector::replaceMap");
568 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
569 if (bmap.is_null() == true) {
570 // if this has more than 1 sub blocks but "map" is not a blocked map, they are very likely not compatible
571 if (this->getBlockedMap()->getNumMaps() > 1) {
572 throw Xpetra::Exceptions::RuntimeError(
573 "BlockedMultiVector::replaceMap: map is not of type BlockedMap. General implementation not available, yet.");
574 TEUCHOS_UNREACHABLE_RETURN();
575 }
576 // special case: this is a blocked map with only one block
577 // TODO add more debug check (especially for Thyra mode)
578 std::vector<Teuchos::RCP<const Map>> subMaps(1, map);
579 map_ = Teuchos::rcp(new BlockedMap(map, subMaps, this->getBlockedMap()->getThyraMode()));
580 this->getMultiVector(0)->replaceMap(map);
581 return;
582 }
583 RCP<const BlockedMap> mybmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map_);
584
585 XPETRA_TEST_FOR_EXCEPTION(
586 mybmap->getThyraMode() != bmap->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::replaceMap: inconsistent Thyra mode");
587 map_ = bmap;
588 for (size_t r = 0; r < map_->getNumMaps(); r++)
589 getMultiVector(r)->replaceMap(bmap->getMap(r, map_->getThyraMode()));
590}
591
592template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
593void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
594 doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* source */,
595 const Import& /* importer */,
596 CombineMode /* CM */) {
597 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
598}
599
600template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
601void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
602 doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* dest */,
603 const Import& /* importer */,
604 CombineMode /* CM */) {
605 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
606}
607
608template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
609void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
610 doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* source */,
611 const Export& /* exporter */,
612 CombineMode /* CM */) {
613 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
614}
615
616template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
617void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
618 doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* dest */,
619 const Export& /* exporter */,
620 CombineMode /* CM */) {
621 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
622}
623
624template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
625void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
626 setSeed(unsigned int seed) {
627 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
628 getMultiVector(r)->setSeed(seed);
629 }
630}
631
632template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
633void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
634 randomize(bool bUseXpetraImplementation) {
635 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
636 getMultiVector(r)->randomize(bUseXpetraImplementation);
637 }
638}
639
640template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
641void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
642 randomize(const Scalar& minVal, const Scalar& maxVal, bool bUseXpetraImplementation) {
643 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
644 getMultiVector(r)->randomize(minVal, maxVal, bUseXpetraImplementation);
645 }
646}
647
648template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
653
654template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
659
660template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
661Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
662BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
663 getMap() const {
664 XPETRA_MONITOR("BlockedMultiVector::getMap");
665 return map_;
666}
667
668template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
669Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>>
670BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
671 getBlockedMap() const {
672 XPETRA_MONITOR("BlockedMultiVector::getBlockedMap");
673 return map_;
674}
675
676template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
677Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
678BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
679 getMultiVector(size_t r) const {
680 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r)");
681 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
682 std::out_of_range,
683 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
684 << " partial blocks.");
685 return vv_[r];
686}
687
688template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
689Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
690BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
691 getMultiVector(size_t r, bool bThyraMode) const {
692 XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r,bThyraMode)");
693 TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
694 std::out_of_range,
695 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
696 << " partial blocks.");
697 XPETRA_TEST_FOR_EXCEPTION(
698 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::getMultiVector: inconsistent Thyra mode");
699 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
700 return vv_[r];
701}
702
703template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
704void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
705 setMultiVector(size_t r,
706 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> v,
707 bool bThyraMode) {
708 // The map of the MultiVector should be the same as the stored submap
709 // In thyra mode the vectors should live on the thyra maps
710 // in xpetra mode the should live in the xpetra maps
711 // that should be also ok in the nested case for thyra (if the vectors are distributed accordingly)
712 XPETRA_MONITOR("BlockedMultiVector::setMultiVector");
713 XPETRA_TEST_FOR_EXCEPTION(r >= map_->getNumMaps(),
714 std::out_of_range,
715 "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
716 XPETRA_TEST_FOR_EXCEPTION(numVectors_ != v->getNumVectors(),
717 Xpetra::Exceptions::RuntimeError,
718 "The BlockedMultiVectors expects " << getNumVectors() << " vectors. The provided partial multivector has "
719 << v->getNumVectors() << " vectors.");
720 XPETRA_TEST_FOR_EXCEPTION(
721 map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::setMultiVector: inconsistent Thyra mode");
722 (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
723 Teuchos::RCP<MultiVector> vv = Teuchos::rcp_const_cast<MultiVector>(v);
724 TEUCHOS_TEST_FOR_EXCEPTION(vv == Teuchos::null, Xpetra::Exceptions::RuntimeError, "Partial vector must not be Teuchos::null");
725 vv_[r] = vv;
726}
727
728template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
729Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
730BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
731 Merge() const {
732 XPETRA_MONITOR("BlockedMultiVector::Merge");
733 using Teuchos::RCP;
734
735 RCP<MultiVector> v = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getFullMap(), getNumVectors());
736 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
737 RCP<MultiVector> vi = getMultiVector(r);
738 RCP<BlockedMultiVector> bvi = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vi);
739 if (bvi.is_null() == true) {
740 this->InsertVector(vi, r, v, map_->getThyraMode());
741 } else {
742 RCP<MultiVector> mvi = bvi->Merge();
743 this->InsertVector(mvi, r, v, map_->getThyraMode());
744 }
745 }
746
747 // TODO plausibility checks
748
749 return v;
750}
751
752template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
753void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
754 assign(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& rhs) {
755 Teuchos::RCP<const MultiVector> rcpRhs = Teuchos::rcpFromRef(rhs);
756 Teuchos::RCP<const BlockedMultiVector> bRhs = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpRhs);
757 if (bRhs == Teuchos::null)
758 throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: argument is not a blocked multi vector.");
759
760 map_ = Teuchos::rcp(new BlockedMap(*(bRhs->getBlockedMap())));
761 vv_ = std::vector<Teuchos::RCP<MultiVector>>(map_->getNumMaps());
762 for (size_t r = 0; r < map_->getNumMaps(); ++r) {
763 // extract source vector (is of type MultiVector or BlockedMultiVector)
764 RCP<MultiVector> src = bRhs->getMultiVector(r, map_->getThyraMode());
765
766 // create new (empty) multivector (is of type MultiVector or BlockedMultiVector)
767 RCP<MultiVector> vv = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(
768 map_->getMap(r, bRhs->getBlockedMap()->getThyraMode()), rcpRhs->getNumVectors(), true);
769
770 // check type of source and target vector
771 RCP<BlockedMultiVector> bsrc = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(src);
772 RCP<BlockedMultiVector> bvv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vv);
773 if (bsrc.is_null() == true && bvv.is_null() == true) {
774 *vv = *src; // deep copy
775 } else if (bsrc.is_null() == true && bvv.is_null() == false) {
776 // special case: source vector is a merged MultiVector but target vector is blocked
777 *vv = *src; // deep copy (is this a problem???)
778 } else if (bsrc.is_null() == false && bvv.is_null() == true) {
779 // special case: source vector is blocked but target vector is a merged MultiVector
780 // This is a problem and only works if bsrc has only one block
781 if (bsrc->getBlockedMap()->getNumMaps() > 1) {
782 throw Xpetra::Exceptions::RuntimeError(
783 "BlockedMultiVector::assign: source vector is of type BlockedMultiVector (with more than "
784 "1 blocks) and target is a MultiVector.");
785 TEUCHOS_UNREACHABLE_RETURN();
786 }
787 RCP<MultiVector> ssrc = bsrc->getMultiVector(0, map_->getThyraMode());
788 XPETRA_TEST_FOR_EXCEPTION(ssrc.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: cannot extract vector");
789 XPETRA_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<BlockedMultiVector>(ssrc) != Teuchos::null,
790 Xpetra::Exceptions::RuntimeError,
791 "BlockedMultiVector::assign: sub block must not be of type BlockedMultiVector.");
792 *vv = *ssrc;
793 } else {
794 // this should work (no exception necessary, i guess)
795 XPETRA_TEST_FOR_EXCEPTION(bsrc->getBlockedMap()->getNumMaps() != bvv->getBlockedMap()->getNumMaps(),
796 Xpetra::Exceptions::RuntimeError,
797 "BlockedMultiVector::assign: source and target are BlockedMultiVectors with a different number of submaps.");
798 *vv = *src; // deep copy
799 }
800 vv_[r] = vv;
801 }
802 numVectors_ = rcpRhs->getNumVectors();
803}
804
805template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
812
813template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
820
821template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
822RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
823BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
824 ExtractVector(RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& full,
825 size_t block,
826 bool bThyraMode) const {
827 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
828 std::out_of_range,
829 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
830 << " partial blocks.");
831 XPETRA_TEST_FOR_EXCEPTION(
832 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
833 RCP<const BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(full);
834 if (bfull.is_null() == true) {
835 // standard case: full is not of type BlockedMultiVector
836 // first extract partial vector from full vector (using xpetra style GIDs)
837 const RCP<MultiVector> vv =
838 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
839 // if(bThyraMode == false) {
840 // ExtractVector(*full, block, *vv);
841 // return vv;
842 //} else {
843 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
844 RCP<MultiVector> rcpNonConstFull = Teuchos::rcp_const_cast<MultiVector>(full);
845 rcpNonConstFull->replaceMap(map_->getImporter(block)->getSourceMap());
846 ExtractVector(*rcpNonConstFull, block, *vv);
847 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
848 Xpetra::Exceptions::RuntimeError,
849 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
850 "created using Thyra-style numbered submaps.");
851 if (bThyraMode == true)
852 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
853 rcpNonConstFull->replaceMap(oldThyMapFull);
854 return vv;
855 //}
856 } else {
857 // special case: full is of type BlockedMultiVector
858 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
859 Xpetra::Exceptions::RuntimeError,
860 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
861 << bfull->getBlockedMap()->getNumMaps()
862 << " (number of blocks in BlockedMultiVector)");
863 return bfull->getMultiVector(block, bThyraMode);
864 }
865}
866
867template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
868RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
869BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
870 ExtractVector(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& full,
871 size_t block,
872 bool bThyraMode) const {
873 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
874 std::out_of_range,
875 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
876 << " partial blocks.");
877 XPETRA_TEST_FOR_EXCEPTION(
878 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
879 RCP<BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(full);
880 if (bfull.is_null() == true) {
881 // standard case: full is not of type BlockedMultiVector
882 // first extract partial vector from full vector (using xpetra style GIDs)
883 const RCP<MultiVector> vv =
884 Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
885 // if(bThyraMode == false) {
886 // ExtractVector(*full, block, *vv);
887 // return vv;
888 //} else {
889 RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
890 full->replaceMap(map_->getImporter(block)->getSourceMap());
891 ExtractVector(*full, block, *vv);
892 TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
893 Xpetra::Exceptions::RuntimeError,
894 "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
895 "created using Thyra-style numbered submaps.");
896 if (bThyraMode == true)
897 vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
898 full->replaceMap(oldThyMapFull);
899 return vv;
900 //}
901 } else {
902 // special case: full is of type BlockedMultiVector
903 XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
904 Xpetra::Exceptions::RuntimeError,
905 "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
906 << bfull->getBlockedMap()->getNumMaps()
907 << " (number of blocks in BlockedMultiVector)");
908 return bfull->getMultiVector(block, bThyraMode);
909 }
910}
911
912template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
913void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
914 ExtractVector(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& full,
915 size_t block,
916 Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& partial) const {
917 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
918 std::out_of_range,
919 "ExtractVector: Error, block = " << block << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
920 << " partial blocks.");
921 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: maps_[" << block << "] is null");
922 partial.doImport(full, *(map_->getImporter(block)), Xpetra::INSERT);
923}
924
925template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
926void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
927 InsertVector(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& partial,
928 size_t block,
929 Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& full,
930 bool bThyraMode) const {
931 XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
932 std::out_of_range,
933 "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
934 << " partial blocks.");
935 XPETRA_TEST_FOR_EXCEPTION(
936 map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
937 XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
938 Xpetra::Exceptions::RuntimeError,
939 "MapExtractor::InsertVector: InsertVector in Thyra-style numbering only possible if MapExtractor has been created "
940 "using Thyra-style numbered submaps.");
941 if (bThyraMode) {
942 // NOTE: the importer objects in the BlockedMap are always using Xpetra GIDs (or Thyra style Xpetra GIDs)
943 // The source map corresponds to the full map (in Xpetra GIDs) starting with GIDs from zero. The GIDs are consecutive in Thyra mode
944 // The target map is the partial map (in the corresponding Xpetra GIDs)
945
946 // TODO can we skip the Export call in special cases (i.e. Src = Target map, same length, etc...)
947
948 // store original GIDs (could be Thyra GIDs)
949 RCP<const MultiVector> rcpPartial = Teuchos::rcpFromRef(partial);
950 RCP<MultiVector> rcpNonConstPartial = Teuchos::rcp_const_cast<MultiVector>(rcpPartial);
951 RCP<const Map> oldThyMapPartial = rcpNonConstPartial->getMap(); // temporarely store map of partial
952 RCP<const Map> oldThyMapFull = full.getMap(); // temporarely store map of full
953
954 // check whether getMap(block,false) is identical to target map of importer
955 // XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false)->isSameAs(*(map_->getImporter(block)->getTargetMap()))==false,
956 // Xpetra::Exceptions::RuntimeError,
957 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of partial vector are not identical to target Map
958 // of Importer. This should not be.");
959
960 // XPETRA_TEST_FOR_EXCEPTION(full.getMap()->isSameAs(*(map_->getImporter(block)->getSourceMap()))==false,
961 // Xpetra::Exceptions::RuntimeError,
962 // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of full vector are not identical to source Map of
963 // Importer. This should not be.");
964
965 rcpNonConstPartial->replaceMap(map_->getMap(block, false)); // temporarely switch to xpetra-style map
966 full.replaceMap(map_->getImporter(block)->getSourceMap()); // temporarely switch to Xpetra GIDs
967
968 // do the Export
969 full.doExport(*rcpNonConstPartial, *(map_->getImporter(block)), Xpetra::INSERT);
970
971 // switch back to original maps
972 full.replaceMap(oldThyMapFull); // reset original map (Thyra GIDs)
973 rcpNonConstPartial->replaceMap(oldThyMapPartial); // change map back to original map
974 } else {
975 // Xpetra style numbering
976 full.doExport(partial, *(map_->getImporter(block)), Xpetra::INSERT);
977 }
978}
979
980template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
981void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
982 InsertVector(RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partial,
983 size_t block,
984 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> full,
985 bool bThyraMode) const {
986 RCP<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bfull =
987 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
988 if (bfull.is_null() == true)
989 InsertVector(*partial, block, *full, bThyraMode);
990 else {
991 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
992
993#if 0
994 // WCMCLEN - ETI: MultiVector doesn't have a setMultiVector method,
995 // WCMCLEN - ETI: but BlockedMultiVector does... should this be bfull->...?
996 full->setMultiVector(block, partial, bThyraMode);
997#else
998 throw Xpetra::Exceptions::RuntimeError("MultiVector::setMultiVector() is not implemented.");
999#endif
1000 }
1001}
1002
1003template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1004void BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1005 InsertVector(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> partial,
1006 size_t block,
1007 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> full,
1008 bool bThyraMode) const {
1009 RCP<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bfull =
1010 Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(full);
1011 if (bfull.is_null() == true)
1012 InsertVector(*partial, block, *full, bThyraMode);
1013 else {
1014 XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1015 bfull->setMultiVector(block, partial, bThyraMode);
1016 }
1017}
1018
1019} // namespace Xpetra
1020
1021#endif // XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
XPETRA_MONITOR
#define XPETRA_MONITOR(funcName)
XPETRA_TEST_FOR_EXCEPTION
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra_BlockedMultiVector_decl.hpp
Teuchos::Array::begin
iterator begin()
Teuchos::Array::end
iterator end()
Teuchos::RCP::is_null
bool is_null() const
Xpetra::BlockedMultiVector::map_
Teuchos::RCP< const BlockedMap > map_
blocked map containing the sub block maps (either thyra or xpetra mode)
Definition Xpetra_BlockedMultiVector_decl.hpp:295
Xpetra::BlockedMultiVector::BlockedMultiVector
BlockedMultiVector(const Teuchos::RCP< const BlockedMap > &map, size_t NumVectors, bool zeroOut=true)
Constructor.
Definition Xpetra_BlockedMultiVector_def.hpp:23
Xpetra::BlockedMultiVector::getVector
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.
Definition Xpetra_BlockedMultiVector_def.hpp:179
Xpetra::BlockedMultiVector::getNumVectors
virtual size_t getNumVectors() const
Number of columns in the multivector.
Definition Xpetra_BlockedMultiVector_def.hpp:512
Xpetra::BlockedMultiVector::ExtractVector
void ExtractVector(RCP< const MultiVector > &full, size_t block, RCP< MultiVector > &partial) const
Definition Xpetra_BlockedMultiVector_def.hpp:807
Xpetra::BlockedMultiVector::InsertVector
void InsertVector(const MultiVector &partial, size_t block, MultiVector &full, bool bThyraMode=false) const
Definition Xpetra_BlockedMultiVector_def.hpp:927
Xpetra::BlockedMultiVector::numVectors_
size_t numVectors_
number of vectors (columns in multi vector)
Definition Xpetra_BlockedMultiVector_decl.hpp:297
Xpetra::BlockedMultiVector::description
virtual std::string description() const
A simple one-line description of this object.
Definition Xpetra_BlockedMultiVector_def.hpp:552
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Tpetra::KokkosClassic::DefaultNode::DefaultNodeType >::getBlockedMap
Teuchos::RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Tpetra::KokkosClassic::DefaultNode::DefaultNodeType > > getBlockedMap() const
Definition Xpetra_BlockedMultiVector_def.hpp:671
Xpetra::BlockedMultiVector::vv_
std::vector< Teuchos::RCP< MultiVector > > vv_
array containing RCPs of the partial vectors
Definition Xpetra_BlockedMultiVector_decl.hpp:296
Xpetra::DistObject::doExport
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").
Xpetra::DistObject::getMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Xpetra::DistObject::doImport
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").
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition Xpetra_Exceptions.hpp:65
Xpetra::MultiVectorFactory::Build
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.
Definition Xpetra_MultiVectorFactory_def.hpp:24
Xpetra::MultiVector::replaceMap
virtual void replaceMap(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map)=0
Xpetra::MultiVector::Xpetra_randomize
virtual void Xpetra_randomize()
Set multi-vector values to random numbers. XPetra implementation.
Definition Xpetra_MultiVector_def.hpp:32
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::FancyOStream
basic_FancyOStream< char > FancyOStream
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra::global_size_t
size_t global_size_t
Global size_t object.
Definition Xpetra_ConfigDefs.hpp:138
Xpetra::CombineMode
CombineMode
Xpetra::Combine Mode enumerable type.
Definition Xpetra_ConfigDefs.hpp:176
Teuchos::ScalarTraits::squareroot
static T squareroot(T x)