Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
BlockMultiVector_def.hpp
1#ifndef BlockMultiVector_DEF_hpp
2#define BlockMultiVector_DEF_hpp
3
4#include "MultiVector.hpp"
5#include "BlockMap.hpp"
6
15
16
17namespace FEDD {
18template <class SC, class LO, class GO, class NO>
19BlockMultiVector<SC,LO,GO,NO>::BlockMultiVector():
20blockMultiVector_(0),
21blockMap_(),
22mergedMultiVector_(),
23mergedMap_()
24{
25 blockMap_ = Teuchos::rcp( new BlockMap_Type( 0 ) );
26}
27
28template <class SC, class LO, class GO, class NO>
29BlockMultiVector<SC,LO,GO,NO>::BlockMultiVector(UN size):
30blockMultiVector_(size),
31blockMap_(),
32mergedMultiVector_(),
33mergedMap_()
34{
35 blockMap_ = Teuchos::rcp( new BlockMap_Type( size ) );
36}
37
38template <class SC, class LO, class GO, class NO>
39BlockMultiVector<SC,LO,GO,NO>::BlockMultiVector( BlockMultiVectorPtr_Type bMVIn):
40blockMultiVector_(bMVIn->size()),
41blockMap_(),
42mergedMultiVector_(),
43mergedMap_()
44{
45 blockMap_ = Teuchos::rcp( new BlockMap_Type( bMVIn->size() ) );
46 for (UN i=0; i<bMVIn->size(); i++) {
47 //copy MV and then insert
48 MultiVectorConstPtr_Type mvTmp = bMVIn->getBlock(i);
49 MultiVectorPtr_Type mvCopy = Teuchos::rcp(new MultiVector_Type( mvTmp ) );
50 this->addBlock( mvCopy, i );
51 }
52}
53
54template <class SC, class LO, class GO, class NO>
55BlockMultiVector<SC,LO,GO,NO>::BlockMultiVector( std::vector<MapConstPtr_Type>& maps, int numMV ):
56blockMultiVector_( maps.size() ),
57blockMap_(),
58mergedMultiVector_(),
59mergedMap_()
60{
61 blockMap_ = Teuchos::rcp( new BlockMap_Type( maps.size() ) );
62 buildFromMaps( maps, numMV );
63
64
65}
66
67template <class SC, class LO, class GO, class NO>
68BlockMultiVector<SC,LO,GO,NO>::BlockMultiVector( BlockMapConstPtr_Type blockMap, int numMV ):
69blockMultiVector_( blockMap->size() ),
70blockMap_( ),
71mergedMultiVector_(),
72mergedMap_()
73{
74
75 blockMap_ = Teuchos::rcp_const_cast<BlockMap_Type>( blockMap );
76 buildFromBlockMap( blockMap, numMV );
77}
78
79template <class SC, class LO, class GO, class NO>
80BlockMultiVector<SC,LO,GO,NO>::~BlockMultiVector(){
81
82}
83
84template <class SC, class LO, class GO, class NO>
85void BlockMultiVector<SC,LO,GO,NO>::resize(UN size){
86
87 blockMultiVector_.resize( size );
88 blockMap_->resize(size);
89 mergedMultiVector_.reset();
90 mergedMap_.reset();
91 localBlockOffsets_.resize( size );
92 globalBlockOffsets_.resize( size );
93
94}
95
96
97template <class SC, class LO, class GO, class NO>
98void BlockMultiVector<SC,LO,GO,NO>::buildFromMaps( std::vector<MapConstPtr_Type>& maps, int numMV ){
99
100 for (int i=0; i<maps.size(); i++){
101 blockMap_->addBlock( maps[i], i );
102 MultiVectorPtr_Type mv = Teuchos::rcp( new MultiVector_Type( maps[i], numMV ) );
103 this->addBlock( mv, i );
104 }
105}
106
107template <class SC, class LO, class GO, class NO>
108void BlockMultiVector<SC,LO,GO,NO>::buildFromBlockMap( BlockMapConstPtr_Type blockMap, int numMV ){
109
110 for (int i=0; i<blockMap->size(); i++){
111 MultiVectorPtr_Type mv = Teuchos::rcp( new MultiVector_Type( blockMap->getBlock( i ), numMV ) );
112 this->addBlock( mv, i );
113 }
114}
115
116template <class SC, class LO, class GO, class NO>
117UN BlockMultiVector<SC,LO,GO,NO>::getNumVectors() const{
118 TEUCHOS_TEST_FOR_EXCEPTION(blockMultiVector_.size()==0, std::logic_error,"No MultiVector in BlockMultiVector.")
119 TEUCHOS_TEST_FOR_EXCEPTION(blockMultiVector_[0].is_null(), std::runtime_error,"MultiVector in BlockMultiVector is null.")
120 return blockMultiVector_[0]->getNumVectors();
121}
122
123template <class SC, class LO, class GO, class NO>
124int BlockMultiVector<SC,LO,GO,NO>::size() const{
125 return blockMultiVector_.size();
126}
127
128template <class SC, class LO, class GO, class NO>
129typename BlockMultiVector<SC,LO,GO,NO>::MultiVectorConstPtr_Type BlockMultiVector<SC,LO,GO,NO>::getBlock(int i) const{
130
131 TEUCHOS_TEST_FOR_EXCEPTION( (blockMultiVector_.size()-1) < i, std::logic_error,"Block in BlockMultiVector does not exist.")
132 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_[0].is_null(), std::runtime_error, "Block in BlockMultiVector is null.")
133 return blockMultiVector_[i];
134}
135
136template <class SC, class LO, class GO, class NO>
137typename BlockMultiVector<SC,LO,GO,NO>::MultiVectorPtr_Type BlockMultiVector<SC,LO,GO,NO>::getBlockNonConst(int i){
138
139 TEUCHOS_TEST_FOR_EXCEPTION( (blockMultiVector_.size()-1) < i, std::logic_error,"Block in BlockMultiVector does not exist.")
140 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_[0].is_null(), std::runtime_error, "Block in BlockMultiVector is null.")
141 return blockMultiVector_[i];
142}
143
144template <class SC, class LO, class GO, class NO>
145void BlockMultiVector<SC,LO,GO,NO>::addBlock(const MultiVectorPtr_Type& multiVector, int i){
146
147 TEUCHOS_TEST_FOR_EXCEPTION( multiVector.is_null(), std::runtime_error,"MultiVector which you want to add to BlockMultiVector is null.")
148 if ( blockMultiVector_.size()>0 && !blockMultiVector_[0].is_null() )
149 TEUCHOS_TEST_FOR_EXCEPTION( multiVector->getNumVectors()!=blockMultiVector_[0]->getNumVectors(), std::logic_error,"MultiVectors for BlockMultiVector have different numbers of vectors.");
150
151 if (i>blockMultiVector_.size()-1)
152 blockMultiVector_.resize( blockMultiVector_.size()+1 );
153 //Do we need a warning here?
154 blockMultiVector_[i] = multiVector;
155
156 blockMap_->addBlock( multiVector->getMap(), i );
157
158}
159
160template <class SC, class LO, class GO, class NO>
161void BlockMultiVector<SC,LO,GO,NO>::merge(){
162 if ( mergedMap_.is_null() ) {
163 blockMap_->merge();
164 mergedMap_ = Teuchos::rcp_const_cast<Map_Type>(blockMap_->getMergedMap());
165 }
166
167
168 mergedMultiVector_ = Teuchos::rcp( new MultiVector_Type( mergedMap_, blockMultiVector_[0]->getNumVectors() ) );
169 this->determineLocalOffsets();
170 this->determineGlobalOffsets();
171
172 for (UN i=0; i<blockMultiVector_.size(); i++) {
173 if ( !blockMultiVector_[i].is_null() )
174 this->mergeBlock( i );
175 else
176 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, "MultiVector in BlockMultiVector is null.")
177 }
178}
179
180template <class SC, class LO, class GO, class NO>
181void BlockMultiVector<SC,LO,GO,NO>::mergeBlock(UN block){
182
183 MultiVectorPtr_Type mv = blockMultiVector_[block];
184 Teuchos::ArrayRCP<const SC> values;
185 GO offset = globalBlockOffsets_[block];
186 UN numVectors = mv->getNumVectors();
187 for (UN j=0; j<numVectors; j++) {
188 values = mv->getData(j);
189 for (UN i=0; i<values.size(); i++) {
190 mergedMultiVector_->replaceGlobalValue( mv->getMap()->getGlobalElement(i) + offset, j, values[i]);
191 }
192 }
193}
194
195template <class SC, class LO, class GO, class NO>
196void BlockMultiVector<SC,LO,GO,NO>::split(){
197 TEUCHOS_TEST_FOR_EXCEPTION( mergedMultiVector_.is_null(), std::runtime_error, "MergedMultiVector is null and we cannot use split on it. Use merge() first to generate the MergedMultiVector.")
198
199 for (UN i=0; i<blockMultiVector_.size(); i++) {
200 if ( !blockMultiVector_[i].is_null() )
201 this->splitBlock( i );
202 else
203 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, "MultiVector in BlockMultiVector is null.")
204 }
205}
206
207template <class SC, class LO, class GO, class NO>
208void BlockMultiVector<SC,LO,GO,NO>::splitBlock(UN block){
209
210 MultiVectorPtr_Type mv = blockMultiVector_[block];
211 MapConstPtr_Type map = mv->getMap();
212 Teuchos::ArrayRCP<const SC> values;
213 GO offset = globalBlockOffsets_[block];
214 UN numVectors = mv->getNumVectors();
215 UN numElements = (UN) map->getNodeNumElements();
216 for (UN j=0; j<numVectors; j++) {
217 values = mergedMultiVector_->getData(j);
218 for (UN i=0; i<numElements; i++) {
219 GO globalIndex = map->getGlobalElement(i) + offset;
220 LO index = mergedMultiVector_->getMap()->getLocalElement( globalIndex );
221 mv->replaceGlobalValue( mv->getMap()->getGlobalElement(i), j, values[index]);
222 }
223 }
224}
225
226template <class SC, class LO, class GO, class NO>
227void BlockMultiVector<SC,LO,GO,NO>::determineLocalOffsets(){
228
229 typedef Teuchos::ScalarTraits<LO> LOST;
230 localBlockOffsets_ = Teuchos::ArrayRCP<LO>( blockMultiVector_.size(), LOST::zero() );
231
232 for (UN i=1; i<blockMultiVector_.size() ; i++)
233 localBlockOffsets_[i] = localBlockOffsets_[i-1] + blockMultiVector_[i-1]->getMap()->getMaxLocalIndex() + 1;
234
235}
236
237template <class SC, class LO, class GO, class NO>
238void BlockMultiVector<SC,LO,GO,NO>::determineGlobalOffsets(){
239
240 typedef Teuchos::ScalarTraits<GO> GOST;
241 globalBlockOffsets_ = Teuchos::ArrayRCP<GO>( blockMultiVector_.size(), GOST::zero() );
242
243 for (UN i=1; i<blockMultiVector_.size() ; i++)
244 globalBlockOffsets_[i] = globalBlockOffsets_[i-1] + blockMultiVector_[i-1]->getMap()->getMaxAllGlobalIndex() + 1;
245
246}
247
248template <class SC, class LO, class GO, class NO>
249void BlockMultiVector<SC,LO,GO,NO>::setMergedVector( MultiVectorPtr_Type& mv ){
250
251 mergedMultiVector_ = mv;
252 mergedMap_ = mv->getMapNonConst();
253 this->determineLocalOffsets();
254 this->determineGlobalOffsets();
255
256}
257
258template <class SC, class LO, class GO, class NO>
259Teuchos::RCP< Thyra::MultiVectorBase<SC> > BlockMultiVector<SC,LO,GO,NO>::getThyraMultiVector( ) {
260 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_.size() == 0, std::logic_error,"BlockMultiVector size is 0.")
261 Teuchos::RCP<Thyra::MultiVectorBase<SC> > thyraMV;
262 this->merge();
263 thyraMV = mergedMultiVector_->getThyraMultiVector( );
264 return thyraMV;
265}
266
267template <class SC, class LO, class GO, class NO>
268Teuchos::RCP<const Thyra::MultiVectorBase<SC> > BlockMultiVector<SC,LO,GO,NO>::getThyraMultiVectorConst( ) {
269 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_.size() == 0, std::logic_error,"BlockMultiVector size is 0.")
270 Teuchos::RCP<const Thyra::MultiVectorBase<SC> > thyraMV;
271 this->merge();
272 thyraMV = mergedMultiVector_->getThyraMultiVector( );
273 return thyraMV;
274}
275
276template <class SC, class LO, class GO, class NO>
277Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > BlockMultiVector<SC,LO,GO,NO>::getProdThyraMultiVector( ) {
278
279 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_.size() == 0, std::logic_error,"BlockMultiVector size is 0.")
280
281 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpaces( blockMultiVector_.size() );
282 Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > multiVecs( blockMultiVector_.size() );
283 for (int i=0; i<blockMultiVector_.size(); i++){
284 multiVecs[i] = blockMultiVector_[i]->getThyraMultiVector();
285 vectorSpaces[i] = multiVecs[i]->range();
286 }
287
288 const Teuchos::RCP< Thyra::DefaultProductVectorSpace<SC> > productSpace = Thyra::productVectorSpace<SC>( vectorSpaces );
289
290 return Thyra::defaultProductMultiVector<SC>( productSpace, multiVecs );
291}
292
293template <class SC, class LO, class GO, class NO>
294Teuchos::RCP<const Thyra::ProductMultiVectorBase<SC> > BlockMultiVector<SC,LO,GO,NO>::getThyraProdMultiVectorConst( ) const{
295 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,"We need to implement this.")
296 Teuchos::RCP< const Thyra::ProductMultiVectorBase<SC> > thyraProdMV;
297 return thyraProdMV;
298}
299
300template <class SC, class LO, class GO, class NO>
301void BlockMultiVector<SC,LO,GO,NO>::fromThyraMultiVector( Teuchos::RCP< Thyra::MultiVectorBase<SC> > thyraMV) {
302 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_.size() == 0, std::logic_error,"BlockMultiVector size is 0.")
303 if (blockMultiVector_.size() == 1){
304 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_[0].is_null(), std::runtime_error, "Block in BlockMultiVector is null.")
305 blockMultiVector_[0]->fromThyraMultiVector( thyraMV );
306 }
307 else {
308 if (mergedMultiVector_.is_null())
309 this->merge();
310 mergedMultiVector_->fromThyraMultiVector( thyraMV );
311 this->split();
312 }
313}
314
315
316template <class SC, class LO, class GO, class NO>
317void BlockMultiVector<SC,LO,GO,NO>::fromThyraProdMultiVector( Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > thyraMV) {
318 TEUCHOS_TEST_FOR_EXCEPTION( blockMultiVector_.size() == 0, std::logic_error,"BlockMultiVector size is 0.")
319 for (int i=0; i<this->size(); i++) {
320 this->getBlockNonConst(i)->fromThyraMultiVector( thyraMV->getNonconstMultiVectorBlock( i ) );
321 }
322}
323
324
325template <class SC, class LO, class GO, class NO>
326void BlockMultiVector<SC,LO,GO,NO>::norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &norms) const {
327 typedef Teuchos::ScalarTraits<SC> ST;
328 for (int j=0; j<norms.size(); j++)
329 norms[j] = ST::zero();
330
331 for (int i=0; i<size(); i++) {
332 Teuchos::Array<SC> partialNorm(norms.size());
333 blockMultiVector_[i]->norm2(partialNorm());
334 for (int j=0; j<partialNorm.size(); j++)
335 norms[j] += partialNorm[j]*partialNorm[j];
336 }
337 for (int j=0; j<norms.size(); j++)
338 norms[j] = ST::squareroot( norms[j] );
339}
340
341template <class SC, class LO, class GO, class NO>
342void BlockMultiVector<SC,LO,GO,NO>::normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &norms) const {
343 typedef Teuchos::ScalarTraits<SC> ST;
344 for (int j=0; j<norms.size(); j++)
345 norms[j] = ST::zero();
346
347 for (int i=0; i<size(); i++) {
348 Teuchos::Array<SC> partialNorm(norms.size());
349 blockMultiVector_[i]->normInf(partialNorm());
350 for (int j=0; j<partialNorm.size(); j++)
351 if(partialNorm[j]>norms[j])
352 norms[j] = partialNorm[j];
353 }
354
355}
356
357template <class SC, class LO, class GO, class NO>
358void BlockMultiVector<SC,LO,GO,NO>::dot(BlockMultiVectorConstPtr_Type a, const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &dots) {
359
360 for (int j=0; j<dots.size(); j++)
361 dots[j] = Teuchos::ScalarTraits<SC>::zero();
362
363 Teuchos::Array<SC> partialDots(dots.size());
364 for (int i=0; i<size(); i++) {
365 blockMultiVector_[i]->dot( a->getBlock(i), partialDots() );
366 for (int j=0; j<partialDots.size(); j++)
367 dots[j] += partialDots[j];
368 }
369}
370
371template <class SC, class LO, class GO, class NO>
372typename BlockMultiVector<SC,LO,GO,NO>::BlockMultiVectorPtr_Type BlockMultiVector<SC,LO,GO,NO>::sumColumns() const{
373
374 BlockMultiVectorPtr_Type sumBlockMV = Teuchos::rcp( new BlockMultiVector_Type ( blockMap_, 1 ) );
375
376 for (int i=0; i<this->size(); i++) {
377 MultiVectorConstPtr_Type tmpMV = this->getBlock(i);
378 MultiVectorPtr_Type sumMV = tmpMV->sumColumns();
379 sumBlockMV->addBlock( sumMV, i );
380 }
381 return sumBlockMV;
382}
383
384template <class SC, class LO, class GO, class NO>
385void BlockMultiVector<SC,LO,GO,NO>::update( const SC& alpha, const BlockMultiVector_Type& A, const SC& beta) {
386 //this = alpha*A + beta*this
387 TEUCHOS_TEST_FOR_EXCEPTION( size() != A.size(), std::logic_error,"BlockMultiVector sizes are not equal for update.")
388
389 for (int i=0; i<size(); i++) {
390 MultiVectorConstPtr_Type tmpA = A[i];
391 blockMultiVector_[i]->update( alpha, *tmpA, beta );
392 }
393}
394
395template <class SC, class LO, class GO, class NO>
396void BlockMultiVector<SC,LO,GO,NO>::multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const SC &alpha, BlockMultiVectorConstPtr_Type &A, BlockMultiVectorConstPtr_Type &B, const SC &beta){
397// if (this->getMap()->getCommNonConst()->getRank()==0)
398// std::cout << "### For testing purposes only." << std::endl;
399
400 for (int i=0; i<this->size(); i++){
401 MultiVectorConstPtr_Type a = A->getBlock(i);
402 MultiVectorConstPtr_Type b = B->getBlock(i);
403 this->getBlockNonConst(i)->multiply( transA, transB, alpha, a, b, beta );
404 }
405}
406
407template <class SC, class LO, class GO, class NO>
408void BlockMultiVector<SC,LO,GO,NO>::update( const SC& alpha, const BlockMultiVector_Type& A, const SC& beta , const BlockMultiVector_Type& B, const SC& gamma) {
409 //this = alpha*A + beta*B + gamma*this
410 TEUCHOS_TEST_FOR_EXCEPTION( size() != A.size(), std::logic_error,"BlockMultiVector sizes are not equal for update.")
411
412 for (int i=0; i<size(); i++) {
413 MultiVectorConstPtr_Type tmpA = A[i];
414 MultiVectorConstPtr_Type tmpB = B[i];
415 blockMultiVector_[i]->update( alpha, *tmpA, beta, *tmpB, gamma);
416 }
417}
418
419template <class SC, class LO, class GO, class NO>
420void BlockMultiVector<SC,LO,GO,NO>::putScalar( const SC& alpha ) {
421 for (int i=0; i<size(); i++) {
422 blockMultiVector_[i]->putScalar( alpha );
423 }
424}
425
426template <class SC, class LO, class GO, class NO>
427void BlockMultiVector<SC,LO,GO,NO>::scale( const SC& alpha ) {
428 for (int i=0; i<size(); i++) {
429 blockMultiVector_[i]->scale( alpha );
430 }
431}
432
433template <class SC, class LO, class GO, class NO>
434void BlockMultiVector<SC,LO,GO,NO>::writeMM(std::string fN) const{
435 for (int i=0; i<size(); i++) {
436 if (!blockMultiVector_[i].is_null() ){
437 std::string fileName = fN + std::to_string(i) + ".mm" ;
438 blockMultiVector_[i]->writeMM(fileName);
439 }
440 }
441}
442
443template <class SC, class LO, class GO, class NO>
444void BlockMultiVector<SC,LO,GO,NO>::print(Teuchos::EVerbosityLevel verbLevel){
445
446 for (UN i=0; i<blockMultiVector_.size(); i++) {
447 if ( !blockMultiVector_[i].is_null() ) {
448 blockMultiVector_[i]->print( verbLevel );
449 }
450 }
451}
452
453template <class SC, class LO, class GO, class NO>
454typename BlockMultiVector<SC,LO,GO,NO>::BlockMapConstPtr_Type BlockMultiVector<SC,LO,GO,NO>::getMap() const{
455
456 return blockMap_;
457}
458
459template <class SC, class LO, class GO, class NO>
460typename BlockMultiVector<SC,LO,GO,NO>::MultiVectorConstPtr_Type BlockMultiVector<SC,LO,GO,NO>::getMergedVector(){
461 if (this->size()>1) {
462 this->merge();
463 return mergedMultiVector_;
464 }
465 else
466 return this->getBlockNonConst(0);
467}
468
469template <class SC, class LO, class GO, class NO>
470typename BlockMultiVector<SC,LO,GO,NO>::MultiVectorPtr_Type BlockMultiVector<SC,LO,GO,NO>::getMergedVectorNonConst(){
471 if (this->size()>1) {
472 this->merge();
473 return mergedMultiVector_;
474 }
475 else
476 return this->getBlockNonConst(0);
477}
478
479}
480#endif
FEDD::BlockMultiVector::getMergedVectorNonConst
MultiVectorPtr_Type getMergedVectorNonConst()
Return the merged vector as non-const vector. This finds application in pressure projection.
Definition BlockMultiVector_def.hpp:470
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36