Teko Version of the Day
Loading...
Searching...
No Matches
Teko_BlockedReordering.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include <iostream>
48#include <string>
49#include <vector>
50#include <stack>
51
52#include "Teko_BlockedReordering.hpp"
53#include "Teko_Utilities.hpp"
54
55#include "Teuchos_VerboseObject.hpp"
56#include "Teuchos_RCP.hpp"
57#include "Teuchos_StrUtils.hpp"
58
59#include "Thyra_DefaultProductMultiVector.hpp"
60#include "Thyra_DefaultProductVectorSpace.hpp"
61
62using Teuchos::RCP;
63using Teuchos::rcp;
64using Teuchos::rcp_dynamic_cast;
65using Teuchos::Array;
66
67namespace Teko {
68
69void BlockReorderManager::SetBlock(int blockIndex,int reorder)
70{
71 TEUCHOS_ASSERT(blockIndex<(int) children_.size());
72
73 RCP<BlockReorderManager> child = rcp(new BlockReorderLeaf(reorder));
74
75 children_[blockIndex] = child;
76}
77
90void BlockReorderManager::SetBlock(int blockIndex,const RCP<BlockReorderManager> & reorder)
91{
92 TEUCHOS_ASSERT(blockIndex<(int) children_.size());
93
94 children_[blockIndex] = reorder;
95}
96
97const Teuchos::RCP<BlockReorderManager> BlockReorderManager::GetBlock(int blockIndex)
98{
99 TEUCHOS_ASSERT(blockIndex<(int) children_.size());
100
101 if(children_[blockIndex]==Teuchos::null)
102 children_[blockIndex] = rcp(new BlockReorderManager());
103
104 return children_[blockIndex];
105}
106
107const Teuchos::RCP<const BlockReorderManager> BlockReorderManager::GetBlock(int blockIndex) const
108{
109 TEUCHOS_ASSERT(blockIndex<(int) children_.size());
110
111 return children_[blockIndex];
112}
113
115{
116 // build the string by recursively calling each child
117 std::stringstream ss;
118 ss << "[";
119 for(unsigned int i=0;i<children_.size();i++) {
120 if(children_[i]==Teuchos::null)
121 ss << " <NULL> ";
122 else
123 ss << " " << children_[i]->toString() << " ";
124 }
125 ss << "]";
126
127 return ss.str();
128}
129
131{
132 int max = 0;
133 for(unsigned int i=0;i<children_.size();i++) {
134 // see if current child is larger
135 if(children_[i]!=Teuchos::null) {
136 int subMax = children_[i]->LargestIndex();
137 max = max > subMax ? max : subMax;
138 }
139 }
140
141 return max;
142}
143
144Teuchos::RCP<const Thyra::LinearOpBase<double> >
145buildReorderedLinearOp(const BlockReorderManager & bmm,
146 const Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > & blkOp)
147{
148 return buildReorderedLinearOp(bmm,bmm,blkOp);
149}
151Teuchos::RCP<const Thyra::LinearOpBase<double> >
152buildReorderedLinearOp(const BlockReorderManager & rowMgr,const BlockReorderManager & colMgr,
153 const Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > & blkOp)
154{
155 typedef RCP<const BlockReorderManager> BRMptr;
156
157 int rowSz = rowMgr.GetNumBlocks();
158 int colSz = colMgr.GetNumBlocks();
159
160 if(rowSz==0 && colSz==0) {
161 // both are leaf nodes
162 const BlockReorderLeaf & rowLeaf = dynamic_cast<const BlockReorderLeaf &>(rowMgr);
163 const BlockReorderLeaf & colLeaf = dynamic_cast<const BlockReorderLeaf &>(colMgr);
164
165 // simply return entry in matrix
166 Teko::LinearOp linOp = blkOp->getBlock(rowLeaf.GetIndex(),colLeaf.GetIndex());
167
168 // somehow we need to set this operator up
169 if(linOp==Teuchos::null) {
170 linOp = Thyra::zero(blkOp->productRange()->getBlock(rowLeaf.GetIndex()),
171 blkOp->productDomain()->getBlock(colLeaf.GetIndex()));
172 }
173
174 return linOp;
175 }
176 else if(rowSz==0) {
177 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
178
179 // operator will be rowSz by colSz
180 reBlkOp->beginBlockFill(1,colSz);
181
182 // fill the column entries
183 for(int col=0;col<colSz;col++) {
184 BRMptr colPtr = colMgr.GetBlock(col);
185
186 reBlkOp->setBlock(0,col,buildReorderedLinearOp(rowMgr,*colPtr,blkOp));
187 }
188
189 // done building
190 reBlkOp->endBlockFill();
191
192 return reBlkOp;
193 }
194 else if(colSz==0) {
195 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
196
197 // operator will be rowSz by colSz
198 reBlkOp->beginBlockFill(rowSz,1);
199
200 // fill the row entries
201 for(int row=0;row<rowSz;row++) {
202 BRMptr rowPtr = rowMgr.GetBlock(row);
203
204 reBlkOp->setBlock(row,0,buildReorderedLinearOp(*rowPtr,colMgr,blkOp));
205 }
206
207 // done building
208 reBlkOp->endBlockFill();
209
210 return reBlkOp;
211 }
212 else {
213 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
214
215 // this is the general case
216 TEUCHOS_ASSERT(rowSz>0);
217 TEUCHOS_ASSERT(colSz>0);
218
219 // operator will be rowSz by colSz
220 reBlkOp->beginBlockFill(rowSz,colSz);
221
222 for(int row=0;row<rowSz;row++) {
223 BRMptr rowPtr = rowMgr.GetBlock(row);
224
225 for(int col=0;col<colSz;col++) {
226 BRMptr colPtr = colMgr.GetBlock(col);
227
228 reBlkOp->setBlock(row,col,buildReorderedLinearOp(*rowPtr,*colPtr,blkOp));
229 }
230 }
231
232 // done building
233 reBlkOp->endBlockFill();
234
235 return reBlkOp;
236 }
237}
238
260Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
262 const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<double> > & blkSpc)
263{
264 typedef RCP<const BlockReorderManager> BRMptr;
265
266 int sz = mgr.GetNumBlocks();
267
268 if(sz==0) {
269 // its a leaf nodes
270 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
271
272 // simply return entry in matrix
273 return blkSpc->getBlock(leaf.GetIndex());
274 }
275 else {
276 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
277
278 // loop over each row
279 for(int i=0;i<sz;i++) {
280 BRMptr blkMgr = mgr.GetBlock(i);
281
282 const RCP<const Thyra::VectorSpaceBase<double> > lvs = buildReorderedVectorSpace(*blkMgr,blkSpc);
283
284 vecSpaces.push_back(lvs);
285 }
286
287 // build a vector space
288 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
289 = Thyra::productVectorSpace<double>(vecSpaces);
290
291 // build the vector
292 return vs;
293 }
294}
295
300Teuchos::RCP<Thyra::MultiVectorBase<double> >
302 const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > & blkVec)
303{
304 typedef RCP<const BlockReorderManager> BRMptr;
305
306 int sz = mgr.GetNumBlocks();
307
308 if(sz==0) {
309 // its a leaf nodes
310 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
311
312 // simply return entry in matrix
313 return blkVec->getNonconstMultiVectorBlock(leaf.GetIndex());
314 }
315 else {
316 Array<RCP<Thyra::MultiVectorBase<double> > > multiVecs;
317 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
318
319 // loop over each row
320 for(int i=0;i<sz;i++) {
321 BRMptr blkMgr = mgr.GetBlock(i);
322
323 const RCP<Thyra::MultiVectorBase<double> > lmv = buildReorderedMultiVector(*blkMgr,blkVec);
324 const RCP<const Thyra::VectorSpaceBase<double> > lvs = lmv->range();
325
326 multiVecs.push_back(lmv);
327 vecSpaces.push_back(lvs);
328 }
329
330 // build a vector space
331 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
332 = Thyra::productVectorSpace<double>(vecSpaces);
333
334 // build the vector
335 return Thyra::defaultProductMultiVector<double>(vs,multiVecs);
336 }
337}
338
343Teuchos::RCP<const Thyra::MultiVectorBase<double> >
345 const Teuchos::RCP<const Thyra::ProductMultiVectorBase<double> > & blkVec)
346{
347 typedef RCP<const BlockReorderManager> BRMptr;
348
349 int sz = mgr.GetNumBlocks();
350
351 if(sz==0) {
352 // its a leaf nodes
353 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
354
355 // simply return entry in matrix
356 return blkVec->getMultiVectorBlock(leaf.GetIndex());
357 }
358 else {
359 Array<RCP<const Thyra::MultiVectorBase<double> > > multiVecs;
360 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
361
362 // loop over each row
363 for(int i=0;i<sz;i++) {
364 BRMptr blkMgr = mgr.GetBlock(i);
365
366 const RCP<const Thyra::MultiVectorBase<double> > lmv = buildReorderedMultiVector(*blkMgr,blkVec);
367 const RCP<const Thyra::VectorSpaceBase<double> > lvs = lmv->range();
368
369 multiVecs.push_back(lmv);
370 vecSpaces.push_back(lvs);
371 }
372
373 // build a vector space
374 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
375 = Thyra::productVectorSpace<double>(vecSpaces);
376
377 // build the vector
378 return Thyra::defaultProductMultiVector<double>(vs,multiVecs);
379 }
380}
381
385void buildNonconstFlatMultiVector(const BlockReorderManager & mgr,
386 const RCP<Thyra::MultiVectorBase<double> > & blkVec,
387 Array<RCP<Thyra::MultiVectorBase<double> > > & multivecs,
388 Array<RCP<const Thyra::VectorSpaceBase<double> > > & vecspaces)
389{
390 int sz = mgr.GetNumBlocks();
391
392 if(sz==0) {
393 // its a leaf nodes
394 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
395 int index = leaf.GetIndex();
396
397 // simply return entry in matrix
398 multivecs[index] = blkVec;
399 vecspaces[index] = blkVec->range();
400 }
401 else {
402 const RCP<Thyra::ProductMultiVectorBase<double> > prodMV
403 = rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(blkVec);
404
405 // get flattened elements from each child
406 for(int i=0;i<sz;i++) {
407 const RCP<Thyra::MultiVectorBase<double> > mv = prodMV->getNonconstMultiVectorBlock(i);
408 buildNonconstFlatMultiVector(*mgr.GetBlock(i),mv,multivecs,vecspaces);
409 }
410 }
411
412}
413
417void buildFlatMultiVector(const BlockReorderManager & mgr,
418 const RCP<const Thyra::MultiVectorBase<double> > & blkVec,
419 Array<RCP<const Thyra::MultiVectorBase<double> > > & multivecs,
420 Array<RCP<const Thyra::VectorSpaceBase<double> > > & vecspaces)
421{
422 int sz = mgr.GetNumBlocks();
423
424 if(sz==0) {
425 // its a leaf nodes
426 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
427 int index = leaf.GetIndex();
428
429 // simply return entry in matrix
430 multivecs[index] = blkVec;
431 vecspaces[index] = blkVec->range();
432 }
433 else {
434 const RCP<const Thyra::ProductMultiVectorBase<double> > prodMV
435 = rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<double> >(blkVec);
436
437 // get flattened elements from each child
438 for(int i=0;i<sz;i++) {
439 RCP<const Thyra::MultiVectorBase<double> > mv = prodMV->getMultiVectorBlock(i);
440 buildFlatMultiVector(*mgr.GetBlock(i),mv,multivecs,vecspaces);
441 }
442 }
443
444}
445
450Teuchos::RCP<Thyra::MultiVectorBase<double> >
451buildFlatMultiVector(const BlockReorderManager & mgr,
452 const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > & blkVec)
453{
454 int numBlocks = mgr.LargestIndex()+1;
455
456 Array<RCP<Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
457 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
458
459 // flatten everything into a vector first
460 buildNonconstFlatMultiVector(mgr,blkVec,multivecs,vecspaces);
461
462 // build a vector space
463 const RCP<Thyra::DefaultProductVectorSpace<double> > vs
464 = Thyra::productVectorSpace<double>(vecspaces);
465
466 // build the vector
467 return Thyra::defaultProductMultiVector<double>(vs,multivecs);
468}
469
474Teuchos::RCP<const Thyra::MultiVectorBase<double> >
475buildFlatMultiVector(const BlockReorderManager & mgr,
476 const Teuchos::RCP<const Thyra::ProductMultiVectorBase<double> > & blkVec)
477{
478 int numBlocks = mgr.LargestIndex()+1;
479
480 Array<RCP<const Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
481 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
482
483 // flatten everything into a vector first
484 buildFlatMultiVector(mgr,blkVec,multivecs,vecspaces);
485
486 // build a vector space
487 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
488 = Thyra::productVectorSpace<double>(vecspaces);
489
490 // build the vector
491 return Thyra::defaultProductMultiVector<double>(vs,multivecs);
492}
493
497void buildFlatVectorSpace(const BlockReorderManager & mgr,
498 const RCP<const Thyra::VectorSpaceBase<double> > & blkSpc,
499 Array<RCP<const Thyra::VectorSpaceBase<double> > > & vecspaces)
500{
501 int sz = mgr.GetNumBlocks();
502
503 if(sz==0) {
504 // its a leaf nodes
505 const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
506 int index = leaf.GetIndex();
507
508 // simply return entry in matrix
509 vecspaces[index] = blkSpc;
510 }
511 else {
512 const RCP<const Thyra::ProductVectorSpaceBase<double> > prodSpc
513 = rcp_dynamic_cast<const Thyra::ProductVectorSpaceBase<double> >(blkSpc);
514
515 // get flattened elements from each child
516 for(int i=0;i<sz;i++) {
517 RCP<const Thyra::VectorSpaceBase<double> > space = prodSpc->getBlock(i);
518 buildFlatVectorSpace(*mgr.GetBlock(i),space,vecspaces);
519 }
520 }
521
522}
523
526Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
527buildFlatVectorSpace(const BlockReorderManager & mgr,
528 const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & blkSpc)
529{
530 int numBlocks = mgr.LargestIndex()+1;
531
532 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
533
534 // flatten everything into a vector first
535 buildFlatVectorSpace(mgr,blkSpc,vecspaces);
536
537 // build a vector space
538 return Thyra::productVectorSpace<double>(vecspaces);
539}
540
542// The next three functions are useful for parsing the string
543// description of a BlockReorderManager.
545
546// this function tokenizes a string, breaking out whitespace but saving the
547// brackets [,] as special tokens.
548static void tokenize(std::string srcInput,std::string whitespace,std::string prefer,
549 std::vector<std::string> & tokens)
550{
551 std::string input = srcInput;
552 std::vector<std::string> wsTokens;
553 std::size_t endPos = input.length()-1;
554 while(endPos<input.length()) {
555 std::size_t next = input.find_first_of(whitespace);
556
557
558 // get the sub string
559 std::string s;
560 if(next!=std::string::npos) {
561 s = input.substr(0,next);
562
563 // break out the old substring
564 input = input.substr(next+1,endPos);
565 }
566 else {
567 s = input;
568 input = "";
569 }
570
571 endPos = input.length()-1;
572
573 // add it to the WS tokens list
574 if(s=="") continue;
575 wsTokens.push_back(s);
576 }
577
578 for(unsigned int i=0;i<wsTokens.size();i++) {
579 // get string to break up
580 input = wsTokens[i];
581
582 endPos = input.length()-1;
583 while(endPos<input.length()) {
584 std::size_t next = input.find_first_of(prefer);
585
586 std::string s = input;
587 if(next>0 && next<input.length()) {
588
589 // get the sub string
590 s = input.substr(0,next);
591
592 input = input.substr(next,endPos);
593 }
594 else if(next==0) {
595 // get the sub string
596 s = input.substr(0,next+1);
597
598 input = input.substr(next+1,endPos);
599 }
600 else input = "";
601
602 // break out the old substring
603 endPos = input.length()-1;
604
605 // add it to the tokens list
606 tokens.push_back(s);
607 }
608 }
609}
610
611// this function takes a set of tokens and returns the first "block", i.e. those
612// values (including) brackets that correspond to the first block
613static std::vector<std::string>::const_iterator
614buildSubBlock(std::vector<std::string>::const_iterator begin,
615 std::vector<std::string>::const_iterator end, std::vector<std::string> & subBlock)
616{
617 std::stack<std::string> matched;
618 std::vector<std::string>::const_iterator itr;
619 for(itr=begin;itr!=end;++itr) {
620
621 subBlock.push_back(*itr);
622
623 // push/pop brackets as they are discovered
624 if(*itr=="[") matched.push("[");
625 else if(*itr=="]") matched.pop();
626
627 // found all matching brackets
628 if(matched.empty())
629 return itr;
630 }
631
632 TEUCHOS_ASSERT(matched.empty());
633
634 return itr-1;
635}
636
637// This function takes a tokenized vector and converts it to a block reorder manager
638static RCP<BlockReorderManager> blockedReorderFromTokens(const std::vector<std::string> & tokens)
639{
640 // base case
641 if(tokens.size()==1)
642 return rcp(new Teko::BlockReorderLeaf(Teuchos::StrUtils::atoi(tokens[0])));
643
644 // check first and last character
645 TEUCHOS_ASSERT(*(tokens.begin())=="[")
646 TEUCHOS_ASSERT(*(tokens.end()-1)=="]");
647
648 std::vector<RCP<Teko::BlockReorderManager> > vecRMgr;
649 std::vector<std::string>::const_iterator itr = tokens.begin()+1;
650 while(itr!=tokens.end()-1) {
651 // figure out which tokens are relevant for this block
652 std::vector<std::string> subBlock;
653 itr = buildSubBlock(itr,tokens.end()-1,subBlock);
654
655 // build the child block reorder manager
656 vecRMgr.push_back(blockedReorderFromTokens(subBlock));
657
658 // move the iterator one more
659 itr++;
660 }
661
662 // build the parent reorder manager
663 RCP<Teko::BlockReorderManager> rMgr = rcp(new Teko::BlockReorderManager(vecRMgr.size()));
664 for(unsigned int i=0;i<vecRMgr.size();i++)
665 rMgr->SetBlock(i,vecRMgr[i]);
666
667 return rMgr;
668}
669
671
683Teuchos::RCP<const BlockReorderManager> blockedReorderFromString(std::string & reorder)
684{
685 // vector of tokens to use
686 std::vector<std::string> tokens;
687
688 // manager to be returned
689
690 // build tokens vector
691 tokenize(reorder," \t\n","[]",tokens);
692
693 // parse recursively and build reorder manager
694 Teuchos::RCP<BlockReorderManager> mgr = blockedReorderFromTokens(tokens);
695
696 return mgr;
697}
698
699} // end namespace Teko
Teuchos::RCP< Thyra::MultiVectorBase< double > > buildReorderedMultiVector(const BlockReorderManager &mgr, const Teuchos::RCP< Thyra::ProductMultiVectorBase< double > > &blkVec)
Convert a flat multi vector into a reordered multivector.
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
int GetIndex() const
Get the the index that is stored in this block.
Class that describes how a flat blocked operator should be reordered.
virtual std::string toString() const
For sanities sake, print a readable string.
virtual void SetBlock(int blockIndex, int reorder)
Sets the sublock to a specific index value.
virtual const Teuchos::RCP< BlockReorderManager > GetBlock(int blockIndex)
Get a particular block. If there is no block at this index location return a new one.
virtual int LargestIndex() const
Largest index in this manager.
BlockReorderManager()
Basic empty constructor.
std::vector< Teuchos::RCP< BlockReorderManager > > children_
Definitions of the subblocks.
virtual int GetNumBlocks() const
Gets the number of subblocks.
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > buildReorderedVectorSpace(const BlockReorderManager &mgr, const Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > &blkSpc)
Use the BlockReorderManager to change a flat vector space into a composite vector space.