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