Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
57 namespace Xpetra {
58 
59 
60 template<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 
89 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
90 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
91 BlockedMultiVector(Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>> bmap,
92  Teuchos::RCP<const MultiVector> v)
93 {
94  XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(),
96  "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has "
97  << bmap->getMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements()
98  << ".");
99  XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(),
101  "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has "
102  << bmap->getMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements()
103  << ".");
104  // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError,
105  // "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " <<
106  // bmap->getFullMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements() << ".");
107  // TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError,
108  // "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has " <<
109  // bmap->getFullMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements() << ".");
110 
111  numVectors_ = v->getNumVectors();
112 
113  map_ = bmap;
114 
115  // Extract vector
116  vv_.reserve(bmap->getNumMaps());
117 
118  // add CrsMatrix objects in row,column order
119  for(size_t r = 0; r < bmap->getNumMaps(); ++r)
120  vv_.push_back(this->ExtractVector(v, r, bmap->getThyraMode()));
121 }
122 
123 
135 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
137 BlockedMultiVector(Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>> mapExtractor,
138  Teuchos::RCP<const MultiVector> v)
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 
158 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
160 BlockedMultiVector(const Teuchos::RCP<const BlockedMap>& map,
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 ];
168 }
169 
170 
171 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
173 ~BlockedMultiVector()
174 {
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 
185 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
189 {
190  assign(rhs); // dispatch to protected virtual method
191  return *this;
192 }
193 
194 
195 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
196 void
198 replaceGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */)
199 {
200  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceGlobalValue: Not (yet) supported by BlockedMultiVector.");
201 }
202 
203 
204 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205 void
207 sumIntoGlobalValue(GlobalOrdinal /* globalRow */, size_t /* vectorIndex */, const Scalar& /* value */)
208 {
209  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoGlobalValue: Not (yet) supported by BlockedMultiVector.");
210 }
211 
212 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213 void
215 replaceLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */)
216 {
217  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceLocalValue: Not supported by BlockedMultiVector.");
218 }
219 
220 
221 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
222 void
224 sumIntoLocalValue(LocalOrdinal /* myRow */, size_t /* vectorIndex */, const Scalar& /* value */)
225 {
226  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoLocalValue:Not (yet) supported by BlockedMultiVector.");
227 }
228 
229 
230 
231 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
232 void
234 putScalar(const Scalar& value)
235 {
236  XPETRA_MONITOR("BlockedMultiVector::putScalar");
237  for(size_t r = 0; r < map_->getNumMaps(); r++)
238  {
239  getMultiVector(r)->putScalar(value);
240  }
241 }
242 
243 
244 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
245 Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
247 getVector(size_t j) const
248 {
249  RCP<Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
250  Teuchos::rcp(new Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(this->getBlockedMap(), false));
251 
252  for(size_t r = 0; r < getBlockedMap()->getNumMaps(); r++)
253  {
254  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> subvec =
255  this->getMultiVector(r, this->getBlockedMap()->getThyraMode())->getVector(j);
256  ret->setMultiVector(r, subvec, this->getBlockedMap()->getThyraMode());
257  }
258  return ret;
259 }
260 
261 
262 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
263 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
266 {
267  RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
268  Teuchos::rcp_const_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(getVector(j));
269  return ret;
270 }
271 
272 
273 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
274 Teuchos::ArrayRCP<const Scalar>
276 getData(size_t j) const
277 {
278  if(map_->getNumMaps() == 1)
279  {
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 
287 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
288 Teuchos::ArrayRCP<Scalar>
291 {
292  if(map_->getNumMaps() == 1)
293  {
294  return vv_[ 0 ]->getDataNonConst(j);
295  }
296  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getDataNonConst: Not (yet) supported by BlockedMultiVector.");
297  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
298 }
299 
300 
301 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
302 void
304 dot(const MultiVector& /* A */, const Teuchos::ArrayView<Scalar>& /* dots */) const
305 {
306  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::dot: Not (yet) supported by BlockedMultiVector.");
307 }
308 
309 
310 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
311 void
313 abs(const MultiVector& /* A */)
314 {
315  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::abs: Not (yet) supported by BlockedMultiVector.");
316 }
317 
318 
319 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
320 void
322 reciprocal(const MultiVector& /* A */)
323 {
324  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::reciprocal: Not (yet) supported by BlockedMultiVector.");
325 }
326 
327 
328 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329 void
331 scale(const Scalar& alpha)
332 {
333  XPETRA_MONITOR("BlockedMultiVector::scale (Scalar)");
334  for(size_t r = 0; r < map_->getNumMaps(); ++r)
335  {
336  if(getMultiVector(r) != Teuchos::null)
337  {
338  getMultiVector(r)->scale(alpha);
339  }
340  }
341 }
342 
343 
344 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
345 void
347 scale(Teuchos::ArrayView<const Scalar> alpha)
348 {
349  XPETRA_MONITOR("BlockedMultiVector::scale (Array)");
350  for(size_t r = 0; r < map_->getNumMaps(); ++r)
351  {
352  if(getMultiVector(r) != Teuchos::null)
353  {
354  getMultiVector(r)->scale(alpha);
355  }
356  }
357 }
358 
359 
360 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
361 void
363 update(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++)
382  {
383 
384  XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != bA->getMultiVector(r)->getMap()->getNodeNumElements(),
386  "BlockedMultiVector::update: in subvector "
387  << r << ": Cannot add a vector of (local) length " << bA->getMultiVector(r)->getMap()->getNodeNumElements()
388  << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
389  XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != bA->getMultiVector(r)->getMap()->getGlobalNumElements(),
391  "BlockedMultiVector::update: in subvector "
392  << r << ": Cannot add a vector of length " << bA->getMultiVector(r)->getMap()->getGlobalNumElements()
393  << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
394 
395  // TAW: 12/6 We basically want to do something like:
396  // getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta);
397  // But if the left hand side is a standard MultiVector and bA->getMultiVector(r) is
398  // a blocked MultiVector this will not work, as the update implementation of the standard
399  // multivector cannot deal with blocked multivectors.
400  // The only use case where this could happen is, if the rhs vector is a single block multivector
401  Teuchos::RCP<MultiVector> lmv = getMultiVector(r);
402  Teuchos::RCP<MultiVector> rmv = bA->getMultiVector(r);
403 
404  // check whether lmv/rmv are blocked or not
405  Teuchos::RCP<BlockedMultiVector> blmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(lmv);
406  Teuchos::RCP<BlockedMultiVector> brmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rmv);
407 
408  if(blmv.is_null() == true && brmv.is_null() == false)
409  {
410  // special case: lmv is standard MultiVector but rmv is BlockedMultiVector
411  TEUCHOS_TEST_FOR_EXCEPTION(brmv->getBlockedMap()->getNumMaps() > 1,
413  "BlockedMultiVector::update: Standard MultiVector object does not accept BlockedMultVector object as "
414  "parameter in update call.");
415  lmv->update(alpha, *(brmv->getMultiVector(0)), beta);
416  }
417  else
418  lmv->update(alpha, *rmv, beta);
419  }
420  }
421  else
422  {
423  // A is a MultiVector
424  // If this is a BlockedMultiVector with only one sub-vector of same length we can just update
425  // Otherwise, A is not compatible with this as BlockedMultiVector and we have to extract the vector data from A
426 
427  if(getBlockedMap()->getNumMaps() == 1)
428  {
429  // Actually call "update" on the underlying MultiVector (Epetra or Tpetra).
430  // The maps have to be compatible.
432  getMultiVector(0)->getMap()->isSameAs(*(rcpA->getMap())) == false,
434  "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor).");
435  getMultiVector(0)->update(alpha, *rcpA, beta);
436  }
437  else
438  {
439  // general case: A has to be splitted and subparts have to be extracted and stored in this BlockedMultiVector
440  // XPETRA_TEST_FOR_EXCEPTION(map_->getFullMap()->isSameAs(*(rcpA->getMap()))==false, Xpetra::Exceptions::RuntimeError,
441  // "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor). - Note:
442  // This test might be too strict and can probably be relaxed!");
443  for(size_t r = 0; r < map_->getNumMaps(); r++)
444  {
445  // Call "update" on the subvector. Note, that getMultiVector(r) could return another BlockedMultiVector.
446  // That is, in Thyra mode the maps could differ (local Xpetra versus Thyra style gids)
447  Teuchos::RCP<const MultiVector> part = this->ExtractVector(rcpA, r, map_->getThyraMode());
448  XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != part->getMap()->getNodeNumElements(),
450  "BlockedMultiVector::update: in subvector "
451  << r << ": Cannot add a vector of (local) length " << part->getMap()->getNodeNumElements()
452  << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
453  XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != part->getMap()->getGlobalNumElements(),
455  "BlockedMultiVector::update: in subvector "
456  << r << ": Cannot add a vector of length " << part->getMap()->getGlobalNumElements()
457  << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
458  getMultiVector(r)->update(alpha, *part, beta);
459  }
460  }
461  }
462 }
463 
464 
465 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
466 void
468 update(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).");
483  TEUCHOS_TEST_FOR_EXCEPTION(
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).");
493  TEUCHOS_TEST_FOR_EXCEPTION(
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 
511 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
512 void
514 norm1(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const
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 
534 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
535 void
537 norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const
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 
560 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
561 void
563 normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& norms) const
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 
582 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
583 void
585 meanValue(const Teuchos::ArrayView<Scalar>& /* means */) const
586 {
587  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::meanValue: Not (yet) supported by BlockedMultiVector.");
588 }
589 
590 
591 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
592 void
594 multiply(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 
605 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
606 void
608 elementWiseMultiply(Scalar scalarAB,
610  const MultiVector& B,
611  Scalar scalarThis)
612 {
613  XPETRA_MONITOR("BlockedMultiVector::elementWiseMultiply");
614  XPETRA_TEST_FOR_EXCEPTION(B.getMap()->isSameAs(*(this->getMap())) == false,
616  "BlockedMultiVector::elementWiseMultipy: B must have same blocked map than this.");
617  // XPETRA_TEST_FOR_EXCEPTION(A.getMap()->isSameAs(*(this->getMap()))==false, Xpetra::Exceptions::RuntimeError,
618  // "BlockedMultiVector::elementWiseMultipy: A must have same blocked map than this.");
619  TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getNodeNumElements() != B.getMap()->getNodeNumElements(),
621  "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getNodeNumElements() << " elements, B has "
622  << B.getMap()->getNodeNumElements() << ".");
623  TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getGlobalNumElements() != B.getMap()->getGlobalNumElements(),
625  "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getGlobalNumElements() << " elements, B has "
626  << B.getMap()->getGlobalNumElements() << ".");
627 
628  RCP<const BlockedMap> bmap = getBlockedMap();
629  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rcpA = Teuchos::rcpFromRef(A);
630  RCP<const Xpetra::BlockedVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bA =
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);
635  TEUCHOS_TEST_FOR_EXCEPTION(
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 
653 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
654 size_t
657 {
658  return numVectors_;
659 }
660 
661 
662 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
663 size_t
666 {
667  XPETRA_MONITOR("BlockedMultiVector::getLocalLength()");
668  return map_->getFullMap()->getNodeNumElements();
669 }
670 
671 
672 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
676 {
677  XPETRA_MONITOR("BlockedMultiVector::getGlobalLength()");
678  return map_->getFullMap()->getGlobalNumElements();
679 }
680 
681 
682 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
683 bool
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);
693  RCP<const MultiVector> b = Vb->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 
703 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
704 std::string
706 description() const
707 {
708  return std::string("BlockedMultiVector");
709 }
710 
711 
712 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
713 void
715 describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const
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 
723 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
724 void
726 replaceMap(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.");
737  TEUCHOS_UNREACHABLE_RETURN();
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 
756 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
757 void
760  const Import& /* importer */,
761  CombineMode /* CM */)
762 {
763  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
764 }
765 
766 
767 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
768 void
771  const Import& /* importer */,
772  CombineMode /* CM */)
773 {
774  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
775 }
776 
777 
778 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
779 void
782  const Export& /* exporter */,
783  CombineMode /* CM */)
784 {
785  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
786 }
787 
788 
789 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
790 void
793  const Export& /* exporter */,
794  CombineMode /* CM */)
795 {
796  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
797 }
798 
799 
800 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
801 void
803 setSeed(unsigned int seed)
804 {
805  for(size_t r = 0; r < map_->getNumMaps(); ++r)
806  {
807  getMultiVector(r)->setSeed(seed);
808  }
809 }
810 
811 
812 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
813 void
815 randomize(bool bUseXpetraImplementation)
816 {
817  for(size_t r = 0; r < map_->getNumMaps(); ++r)
818  {
819  getMultiVector(r)->randomize(bUseXpetraImplementation);
820  }
821 }
822 
823 
824 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
825 void
828 {
830 }
831 
832 
833 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
834 Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
836 getMap() const
837 {
838  XPETRA_MONITOR("BlockedMultiVector::getMap");
839  return map_;
840 }
841 
842 
843 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
844 Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>>
847 {
848  XPETRA_MONITOR("BlockedMultiVector::getBlockedMap");
849  return map_;
850 }
851 
852 
853 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
854 Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
856 getMultiVector(size_t r) const
857 {
858  XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r)");
859  TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
860  std::out_of_range,
861  "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
862  << " partial blocks.");
863  return vv_[ r ];
864 }
865 
866 
867 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
868 Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
870 getMultiVector(size_t r, bool bThyraMode) const
871 {
872  XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r,bThyraMode)");
873  TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(),
874  std::out_of_range,
875  "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
876  << " partial blocks.");
878  map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::getMultiVector: inconsistent Thyra mode");
879  (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
880  return vv_[ r ];
881 }
882 
883 
884 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
885 void
887 setMultiVector(size_t r,
889  bool bThyraMode)
890 {
891  // The map of the MultiVector should be the same as the stored submap
892  // In thyra mode the vectors should live on the thyra maps
893  // in xpetra mode the should live in the xpetra maps
894  // that should be also ok in the nested case for thyra (if the vectors are distributed accordingly)
895  XPETRA_MONITOR("BlockedMultiVector::setMultiVector");
896  XPETRA_TEST_FOR_EXCEPTION(r >= map_->getNumMaps(),
897  std::out_of_range,
898  "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
899  XPETRA_TEST_FOR_EXCEPTION(numVectors_ != v->getNumVectors(),
901  "The BlockedMultiVectors expects " << getNumVectors() << " vectors. The provided partial multivector has "
902  << v->getNumVectors() << " vectors.");
904  map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::setMultiVector: inconsistent Thyra mode");
905  (void)bThyraMode; // avoid unused parameter warning when HAVE_XPETRA_DEBUG isn't defined
906  Teuchos::RCP<MultiVector> vv = Teuchos::rcp_const_cast<MultiVector>(v);
907  TEUCHOS_TEST_FOR_EXCEPTION(vv == Teuchos::null, Xpetra::Exceptions::RuntimeError, "Partial vector must not be Teuchos::null");
908  vv_[ r ] = vv;
909 }
910 
911 
912 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
913 Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
915 Merge() const
916 {
917  XPETRA_MONITOR("BlockedMultiVector::Merge");
918  using Teuchos::RCP;
919 
920  RCP<MultiVector> v = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getFullMap(), getNumVectors());
921  for(size_t r = 0; r < map_->getNumMaps(); ++r)
922  {
923  RCP<MultiVector> vi = getMultiVector(r);
924  RCP<BlockedMultiVector> bvi = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vi);
925  if(bvi.is_null() == true)
926  {
927  this->InsertVector(vi, r, v, map_->getThyraMode());
928  }
929  else
930  {
931  RCP<MultiVector> mvi = bvi->Merge();
932  this->InsertVector(mvi, r, v, map_->getThyraMode());
933  }
934  }
935 
936  // TODO plausibility checks
937 
938  return v;
939 }
940 
941 
942 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
943 void
946 {
947  Teuchos::RCP<const MultiVector> rcpRhs = Teuchos::rcpFromRef(rhs);
948  Teuchos::RCP<const BlockedMultiVector> bRhs = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpRhs);
949  if(bRhs == Teuchos::null)
950  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: argument is not a blocked multi vector.");
951 
952  map_ = Teuchos::rcp(new BlockedMap(*(bRhs->getBlockedMap())));
953  vv_ = std::vector<Teuchos::RCP<MultiVector>>(map_->getNumMaps());
954  for(size_t r = 0; r < map_->getNumMaps(); ++r)
955  {
956  // extract source vector (is of type MultiVector or BlockedMultiVector)
957  RCP<MultiVector> src = bRhs->getMultiVector(r, map_->getThyraMode());
958 
959  // create new (empty) multivector (is of type MultiVector or BlockedMultiVector)
961  map_->getMap(r, bRhs->getBlockedMap()->getThyraMode()), rcpRhs->getNumVectors(), true);
962 
963  // check type of source and target vector
964  RCP<BlockedMultiVector> bsrc = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(src);
965  RCP<BlockedMultiVector> bvv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vv);
966  if(bsrc.is_null() == true && bvv.is_null() == true)
967  {
968  *vv = *src; // deep copy
969  }
970  else if(bsrc.is_null() == true && bvv.is_null() == false)
971  {
972  // special case: source vector is a merged MultiVector but target vector is blocked
973  *vv = *src; // deep copy (is this a problem???)
974  }
975  else if(bsrc.is_null() == false && bvv.is_null() == true)
976  {
977  // special case: source vector is blocked but target vector is a merged MultiVector
978  // This is a problem and only works if bsrc has only one block
979  if(bsrc->getBlockedMap()->getNumMaps() > 1)
980  {
981  throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: source vector is of type BlockedMultiVector (with more than "
982  "1 blocks) and target is a MultiVector.");
983  TEUCHOS_UNREACHABLE_RETURN();
984  }
985  RCP<MultiVector> ssrc = bsrc->getMultiVector(0, map_->getThyraMode());
986  XPETRA_TEST_FOR_EXCEPTION(ssrc.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: cannot extract vector");
987  XPETRA_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<BlockedMultiVector>(ssrc) != Teuchos::null,
989  "BlockedMultiVector::assign: sub block must not be of type BlockedMultiVector.");
990  *vv = *ssrc;
991  }
992  else
993  {
994  // this should work (no exception necessary, i guess)
995  XPETRA_TEST_FOR_EXCEPTION(bsrc->getBlockedMap()->getNumMaps() != bvv->getBlockedMap()->getNumMaps(),
997  "BlockedMultiVector::assign: source and target are BlockedMultiVectors with a different number of submaps.");
998  *vv = *src; // deep copy
999  }
1000  vv_[ r ] = vv;
1001  }
1002  numVectors_ = rcpRhs->getNumVectors();
1003 }
1004 
1005 
1006 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1007 void
1010  size_t block,
1012 {
1013  ExtractVector(*full, block, *partial);
1014 }
1015 
1016 
1017 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1018 void
1019 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1021  size_t block,
1023 {
1024  ExtractVector(*full, block, *partial);
1025 }
1026 
1027 
1028 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1029 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1030 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1032  size_t block,
1033  bool bThyraMode) const
1034 {
1035  XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1036  std::out_of_range,
1037  "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1038  << " partial blocks.");
1040  map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1041  RCP<const BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(full);
1042  if(bfull.is_null() == true)
1043  {
1044  // standard case: full is not of type BlockedMultiVector
1045  // first extract partial vector from full vector (using xpetra style GIDs)
1046  const RCP<MultiVector> vv =
1047  Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1048  // if(bThyraMode == false) {
1049  // ExtractVector(*full, block, *vv);
1050  // return vv;
1051  //} else {
1052  RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1053  RCP<MultiVector> rcpNonConstFull = Teuchos::rcp_const_cast<MultiVector>(full);
1054  rcpNonConstFull->replaceMap(map_->getImporter(block)->getSourceMap());
1055  ExtractVector(*rcpNonConstFull, block, *vv);
1056  TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1057  Xpetra::Exceptions::RuntimeError,
1058  "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1059  "created using Thyra-style numbered submaps.");
1060  if(bThyraMode == true)
1061  vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1062  rcpNonConstFull->replaceMap(oldThyMapFull);
1063  return vv;
1064  //}
1065  }
1066  else
1067  {
1068  // special case: full is of type BlockedMultiVector
1069  XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1070  Xpetra::Exceptions::RuntimeError,
1071  "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1072  << bfull->getBlockedMap()->getNumMaps()
1073  << " (number of blocks in BlockedMultiVector)");
1074  return bfull->getMultiVector(block, bThyraMode);
1075  }
1076 }
1077 
1078 
1079 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1080 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1081 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1083  size_t block,
1084  bool bThyraMode) const
1085 {
1086  XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1087  std::out_of_range,
1088  "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1089  << " partial blocks.");
1091  map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1092  RCP<BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(full);
1093  if(bfull.is_null() == true)
1094  {
1095  // standard case: full is not of type BlockedMultiVector
1096  // first extract partial vector from full vector (using xpetra style GIDs)
1097  const RCP<MultiVector> vv =
1098  Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map_->getMap(block, false), full->getNumVectors(), false);
1099  // if(bThyraMode == false) {
1100  // ExtractVector(*full, block, *vv);
1101  // return vv;
1102  //} else {
1103  RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
1104  full->replaceMap(map_->getImporter(block)->getSourceMap());
1105  ExtractVector(*full, block, *vv);
1106  TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1107  Xpetra::Exceptions::RuntimeError,
1108  "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been "
1109  "created using Thyra-style numbered submaps.");
1110  if(bThyraMode == true)
1111  vv->replaceMap(map_->getMap(block, true)); // switch to Thyra-style map
1112  full->replaceMap(oldThyMapFull);
1113  return vv;
1114  //}
1115  }
1116  else
1117  {
1118  // special case: full is of type BlockedMultiVector
1119  XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(),
1120  Xpetra::Exceptions::RuntimeError,
1121  "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be "
1122  << bfull->getBlockedMap()->getNumMaps()
1123  << " (number of blocks in BlockedMultiVector)");
1124  return bfull->getMultiVector(block, bThyraMode);
1125  }
1126 }
1127 
1128 
1129 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1130 void
1131 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1133  size_t block,
1135 {
1136  XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1137  std::out_of_range,
1138  "ExtractVector: Error, block = " << block << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps()
1139  << " partial blocks.");
1140  XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: maps_[" << block << "] is null");
1141  partial.doImport(full, *(map_->getImporter(block)), Xpetra::INSERT);
1142 }
1143 
1144 
1145 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1146 void
1149  size_t block,
1151  bool bThyraMode) const
1152 {
1153  XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(),
1154  std::out_of_range,
1155  "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps()
1156  << " partial blocks.");
1158  map_->getMap(block, false) == null, Xpetra::Exceptions::RuntimeError, "ExtractVector: map_->getmap(" << block << ",false) is null");
1159  XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true,
1160  Xpetra::Exceptions::RuntimeError,
1161  "MapExtractor::InsertVector: InsertVector in Thyra-style numbering only possible if MapExtractor has been created "
1162  "using Thyra-style numbered submaps.");
1163  if(bThyraMode)
1164  {
1165  // NOTE: the importer objects in the BlockedMap are always using Xpetra GIDs (or Thyra style Xpetra GIDs)
1166  // The source map corresponds to the full map (in Xpetra GIDs) starting with GIDs from zero. The GIDs are consecutive in Thyra mode
1167  // The target map is the partial map (in the corresponding Xpetra GIDs)
1168 
1169  // TODO can we skip the Export call in special cases (i.e. Src = Target map, same length, etc...)
1170 
1171  // store original GIDs (could be Thyra GIDs)
1172  RCP<const MultiVector> rcpPartial = Teuchos::rcpFromRef(partial);
1173  RCP<MultiVector> rcpNonConstPartial = Teuchos::rcp_const_cast<MultiVector>(rcpPartial);
1174  RCP<const Map> oldThyMapPartial = rcpNonConstPartial->getMap(); // temporarely store map of partial
1175  RCP<const Map> oldThyMapFull = full.getMap(); // temporarely store map of full
1176 
1177  // check whether getMap(block,false) is identical to target map of importer
1178  // XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false)->isSameAs(*(map_->getImporter(block)->getTargetMap()))==false,
1179  // Xpetra::Exceptions::RuntimeError,
1180  // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of partial vector are not identical to target Map
1181  // of Importer. This should not be.");
1182 
1183  // XPETRA_TEST_FOR_EXCEPTION(full.getMap()->isSameAs(*(map_->getImporter(block)->getSourceMap()))==false,
1184  // Xpetra::Exceptions::RuntimeError,
1185  // "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of full vector are not identical to source Map of
1186  // Importer. This should not be.");
1187 
1188  rcpNonConstPartial->replaceMap(map_->getMap(block, false)); // temporarely switch to xpetra-style map
1189  full.replaceMap(map_->getImporter(block)->getSourceMap()); // temporarely switch to Xpetra GIDs
1190 
1191  // do the Export
1192  full.doExport(*rcpNonConstPartial, *(map_->getImporter(block)), Xpetra::INSERT);
1193 
1194  // switch back to original maps
1195  full.replaceMap(oldThyMapFull); // reset original map (Thyra GIDs)
1196  rcpNonConstPartial->replaceMap(oldThyMapPartial); // change map back to original map
1197  }
1198  else
1199  {
1200  // Xpetra style numbering
1201  full.doExport(partial, *(map_->getImporter(block)), Xpetra::INSERT);
1202  }
1203 }
1204 
1205 
1206 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1207 void
1210  size_t block,
1212  bool bThyraMode) const
1213 {
1214  RCP<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bfull =
1216  if(bfull.is_null() == true)
1217  InsertVector(*partial, block, *full, bThyraMode);
1218  else
1219  {
1220  XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1221 
1222 #if 0
1223  // WCMCLEN - ETI: MultiVector doesn't have a setMultiVector method,
1224  // WCMCLEN - ETI: but BlockedMultiVector does... should this be bfull->...?
1225  full->setMultiVector(block, partial, bThyraMode);
1226 #else
1227  throw Xpetra::Exceptions::RuntimeError("MultiVector::setMultiVector() is not implemented.");
1228 #endif
1229  }
1230 }
1231 
1232 
1233 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1234 void
1235 BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1237  size_t block,
1239  bool bThyraMode) const
1240 {
1241  RCP<Xpetra::BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bfull =
1243  if(bfull.is_null() == true)
1244  InsertVector(*partial, block, *full, bThyraMode);
1245  else
1246  {
1247  XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError, "InsertVector: maps_[" << block << "] is null");
1248  bfull->setMultiVector(block, partial, bThyraMode);
1249  }
1250 }
1251 
1252 
1253 } // namespace Xpetra
1254 
1255 
1256 #endif // XPETRA_BLOCKEDMULTIVECTOR_DEF_HPP
virtual void replaceMap(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map)=0
virtual void doExport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)=0
Export data into this object using an Export object (&quot;forward mode&quot;).
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.
BlockedMultiVector(const Teuchos::RCP< const BlockedMap > &map, size_t NumVectors, bool zeroOut=true)
Constructor.
Exception throws to report errors in the internal logical of the program.
virtual void Xpetra_randomize()
Set multi-vector values to random numbers. XPetra implementation.
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 doImport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM)=0
Import data into this object using an Import object (&quot;forward mode&quot;).
Teuchos::RCP< MultiVector > Merge() const
merge BlockedMultiVector blocks to a single MultiVector
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 Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
virtual void setSeed(unsigned int seed)
Set seed for Random function.
size_t global_size_t
Global size_t object.
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const
Const view of the local values in a particular vector of this multivector.
Teuchos::RCP< MultiVector > getMultiVector(size_t r) const
return partial multivector associated with block row r
virtual Teuchos::ArrayRCP< Scalar > getDataNonConst(size_t j)
View of the local values in a particular vector of this multivector.
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual size_t getNumVectors() const
Number of columns in the multivector.
CombineMode
Xpetra::Combine Mode enumerable type.
#define XPETRA_MONITOR(funcName)
virtual void randomize(bool bUseXpetraImplementation=false)
Set multi-vector values to random numbers.