52 #include "Teko_BlockedReordering.hpp"
55 #include "Teuchos_VerboseObject.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_StrUtils.hpp"
59 #include "Thyra_DefaultProductMultiVector.hpp"
60 #include "Thyra_DefaultProductVectorSpace.hpp"
65 using Teuchos::rcp_dynamic_cast;
70 TEUCHOS_ASSERT(blockIndex < (
int)
children_.size());
90 TEUCHOS_ASSERT(blockIndex < (
int)
children_.size());
96 TEUCHOS_ASSERT(blockIndex < (
int)
children_.size());
98 if (
children_[blockIndex] == Teuchos::null)
105 TEUCHOS_ASSERT(blockIndex < (
int)
children_.size());
112 std::stringstream ss;
114 for (
unsigned int i = 0; i <
children_.size(); i++) {
118 ss <<
" " <<
children_[i]->toString() <<
" ";
127 for (
unsigned int i = 0; i <
children_.size(); i++) {
130 int subMax =
children_[i]->LargestIndex();
131 max = max > subMax ? max : subMax;
138 Teuchos::RCP<const Thyra::LinearOpBase<double> > buildReorderedLinearOp(
140 const Teuchos::RCP<
const Thyra::BlockedLinearOpBase<double> > &blkOp) {
141 return buildReorderedLinearOp(bmm, bmm, blkOp);
144 Teuchos::RCP<const Thyra::LinearOpBase<double> > buildReorderedLinearOp(
146 const Teuchos::RCP<
const Thyra::BlockedLinearOpBase<double> > &blkOp) {
147 typedef RCP<const BlockReorderManager> BRMptr;
152 if (rowSz == 0 && colSz == 0) {
158 Teko::LinearOp linOp = blkOp->getBlock(rowLeaf.
GetIndex(), colLeaf.
GetIndex());
161 if (linOp == Teuchos::null) {
162 linOp = Thyra::zero(blkOp->productRange()->getBlock(rowLeaf.
GetIndex()),
163 blkOp->productDomain()->getBlock(colLeaf.
GetIndex()));
167 }
else if (rowSz == 0) {
168 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
171 reBlkOp->beginBlockFill(1, colSz);
174 for (
int col = 0; col < colSz; col++) {
175 BRMptr colPtr = colMgr.
GetBlock(col);
177 reBlkOp->setBlock(0, col, buildReorderedLinearOp(rowMgr, *colPtr, blkOp));
181 reBlkOp->endBlockFill();
184 }
else if (colSz == 0) {
185 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
188 reBlkOp->beginBlockFill(rowSz, 1);
191 for (
int row = 0; row < rowSz; row++) {
192 BRMptr rowPtr = rowMgr.
GetBlock(row);
194 reBlkOp->setBlock(row, 0, buildReorderedLinearOp(*rowPtr, colMgr, blkOp));
198 reBlkOp->endBlockFill();
202 Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
205 TEUCHOS_ASSERT(rowSz > 0);
206 TEUCHOS_ASSERT(colSz > 0);
209 reBlkOp->beginBlockFill(rowSz, colSz);
211 for (
int row = 0; row < rowSz; row++) {
212 BRMptr rowPtr = rowMgr.
GetBlock(row);
214 for (
int col = 0; col < colSz; col++) {
215 BRMptr colPtr = colMgr.
GetBlock(col);
217 reBlkOp->setBlock(row, col, buildReorderedLinearOp(*rowPtr, *colPtr, blkOp));
222 reBlkOp->endBlockFill();
251 const Teuchos::RCP<
const Thyra::ProductVectorSpaceBase<double> > &blkSpc) {
252 typedef RCP<const BlockReorderManager> BRMptr;
261 return blkSpc->getBlock(leaf.
GetIndex());
263 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
266 for (
int i = 0; i < sz; i++) {
269 const RCP<const Thyra::VectorSpaceBase<double> > lvs =
270 buildReorderedVectorSpace(*blkMgr, blkSpc);
272 vecSpaces.push_back(lvs);
276 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs =
277 Thyra::productVectorSpace<double>(vecSpaces);
290 const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > &blkVec) {
291 typedef RCP<const BlockReorderManager> BRMptr;
300 return blkVec->getNonconstMultiVectorBlock(leaf.
GetIndex());
302 Array<RCP<Thyra::MultiVectorBase<double> > > multiVecs;
303 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
306 for (
int i = 0; i < sz; i++) {
310 const RCP<const Thyra::VectorSpaceBase<double> > lvs = lmv->range();
312 multiVecs.push_back(lmv);
313 vecSpaces.push_back(lvs);
317 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs =
318 Thyra::productVectorSpace<double>(vecSpaces);
321 return Thyra::defaultProductMultiVector<double>(vs, multiVecs);
331 const Teuchos::RCP<
const Thyra::ProductMultiVectorBase<double> > &blkVec) {
332 typedef RCP<const BlockReorderManager> BRMptr;
341 return blkVec->getMultiVectorBlock(leaf.
GetIndex());
343 Array<RCP<const Thyra::MultiVectorBase<double> > > multiVecs;
344 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
347 for (
int i = 0; i < sz; i++) {
350 const RCP<const Thyra::MultiVectorBase<double> > lmv =
352 const RCP<const Thyra::VectorSpaceBase<double> > lvs = lmv->range();
354 multiVecs.push_back(lmv);
355 vecSpaces.push_back(lvs);
359 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs =
360 Thyra::productVectorSpace<double>(vecSpaces);
363 return Thyra::defaultProductMultiVector<double>(vs, multiVecs);
370 void buildNonconstFlatMultiVector(
const BlockReorderManager &mgr,
371 const RCP<Thyra::MultiVectorBase<double> > &blkVec,
372 Array<RCP<Thyra::MultiVectorBase<double> > > &multivecs,
373 Array<RCP<
const Thyra::VectorSpaceBase<double> > > &vecspaces) {
374 int sz = mgr.GetNumBlocks();
378 const BlockReorderLeaf &leaf =
dynamic_cast<const BlockReorderLeaf &
>(mgr);
379 int index = leaf.GetIndex();
382 multivecs[index] = blkVec;
383 vecspaces[index] = blkVec->range();
385 const RCP<Thyra::ProductMultiVectorBase<double> > prodMV =
386 rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(blkVec);
389 for (
int i = 0; i < sz; i++) {
390 const RCP<Thyra::MultiVectorBase<double> > mv = prodMV->getNonconstMultiVectorBlock(i);
391 buildNonconstFlatMultiVector(*mgr.GetBlock(i), mv, multivecs, vecspaces);
399 void buildFlatMultiVector(
const BlockReorderManager &mgr,
400 const RCP<
const Thyra::MultiVectorBase<double> > &blkVec,
401 Array<RCP<
const Thyra::MultiVectorBase<double> > > &multivecs,
402 Array<RCP<
const Thyra::VectorSpaceBase<double> > > &vecspaces) {
403 int sz = mgr.GetNumBlocks();
407 const BlockReorderLeaf &leaf =
dynamic_cast<const BlockReorderLeaf &
>(mgr);
408 int index = leaf.GetIndex();
411 multivecs[index] = blkVec;
412 vecspaces[index] = blkVec->range();
414 const RCP<const Thyra::ProductMultiVectorBase<double> > prodMV =
415 rcp_dynamic_cast<
const Thyra::ProductMultiVectorBase<double> >(blkVec);
418 for (
int i = 0; i < sz; i++) {
419 RCP<const Thyra::MultiVectorBase<double> > mv = prodMV->getMultiVectorBlock(i);
420 buildFlatMultiVector(*mgr.GetBlock(i), mv, multivecs, vecspaces);
429 Teuchos::RCP<Thyra::MultiVectorBase<double> > buildFlatMultiVector(
431 const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > &blkVec) {
434 Array<RCP<Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
435 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
438 buildNonconstFlatMultiVector(mgr, blkVec, multivecs, vecspaces);
441 const RCP<Thyra::DefaultProductVectorSpace<double> > vs =
442 Thyra::productVectorSpace<double>(vecspaces);
445 return Thyra::defaultProductMultiVector<double>(vs, multivecs);
452 Teuchos::RCP<const Thyra::MultiVectorBase<double> > buildFlatMultiVector(
454 const Teuchos::RCP<
const Thyra::ProductMultiVectorBase<double> > &blkVec) {
457 Array<RCP<const Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
458 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
461 buildFlatMultiVector(mgr, blkVec, multivecs, vecspaces);
464 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs =
465 Thyra::productVectorSpace<double>(vecspaces);
468 return Thyra::defaultProductMultiVector<double>(vs, multivecs);
474 void buildFlatVectorSpace(
const BlockReorderManager &mgr,
475 const RCP<
const Thyra::VectorSpaceBase<double> > &blkSpc,
476 Array<RCP<
const Thyra::VectorSpaceBase<double> > > &vecspaces) {
477 int sz = mgr.GetNumBlocks();
481 const BlockReorderLeaf &leaf =
dynamic_cast<const BlockReorderLeaf &
>(mgr);
482 int index = leaf.GetIndex();
485 vecspaces[index] = blkSpc;
487 const RCP<const Thyra::ProductVectorSpaceBase<double> > prodSpc =
488 rcp_dynamic_cast<
const Thyra::ProductVectorSpaceBase<double> >(blkSpc);
491 for (
int i = 0; i < sz; i++) {
492 RCP<const Thyra::VectorSpaceBase<double> > space = prodSpc->getBlock(i);
493 buildFlatVectorSpace(*mgr.GetBlock(i), space, vecspaces);
500 Teuchos::RCP<const Thyra::VectorSpaceBase<double> > buildFlatVectorSpace(
501 const BlockReorderManager &mgr,
502 const Teuchos::RCP<
const Thyra::VectorSpaceBase<double> > &blkSpc) {
503 int numBlocks = mgr.LargestIndex() + 1;
505 Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
508 buildFlatVectorSpace(mgr, blkSpc, vecspaces);
511 return Thyra::productVectorSpace<double>(vecspaces);
521 static void tokenize(std::string srcInput, std::string whitespace, std::string prefer,
522 std::vector<std::string> &tokens) {
523 std::string input = srcInput;
524 std::vector<std::string> wsTokens;
525 std::size_t endPos = input.length() - 1;
526 while (endPos < input.length()) {
527 std::size_t next = input.find_first_of(whitespace);
531 if (next != std::string::npos) {
532 s = input.substr(0, next);
535 input = input.substr(next + 1, endPos);
541 endPos = input.length() - 1;
544 if (s ==
"")
continue;
545 wsTokens.push_back(s);
548 for (
unsigned int i = 0; i < wsTokens.size(); i++) {
552 endPos = input.length() - 1;
553 while (endPos < input.length()) {
554 std::size_t next = input.find_first_of(prefer);
556 std::string s = input;
557 if (next > 0 && next < input.length()) {
559 s = input.substr(0, next);
561 input = input.substr(next, endPos);
562 }
else if (next == 0) {
564 s = input.substr(0, next + 1);
566 input = input.substr(next + 1, endPos);
571 endPos = input.length() - 1;
581 static std::vector<std::string>::const_iterator buildSubBlock(
582 std::vector<std::string>::const_iterator begin, std::vector<std::string>::const_iterator end,
583 std::vector<std::string> &subBlock) {
584 std::stack<std::string> matched;
585 std::vector<std::string>::const_iterator itr;
586 for (itr = begin; itr != end; ++itr) {
587 subBlock.push_back(*itr);
592 else if (*itr ==
"]")
596 if (matched.empty())
return itr;
599 TEUCHOS_ASSERT(matched.empty());
605 static RCP<BlockReorderManager> blockedReorderFromTokens(
const std::vector<std::string> &tokens) {
607 if (tokens.size() == 1)
611 TEUCHOS_ASSERT(*(tokens.begin()) ==
"[")
612 TEUCHOS_ASSERT(*(tokens.end() - 1) == "]");
614 std::vector<RCP<Teko::BlockReorderManager> > vecRMgr;
615 std::vector<std::
string>::const_iterator itr = tokens.begin() + 1;
616 while (itr != tokens.end() - 1) {
618 std::vector<std::string> subBlock;
619 itr = buildSubBlock(itr, tokens.end() - 1, subBlock);
622 vecRMgr.push_back(blockedReorderFromTokens(subBlock));
630 for (
unsigned int i = 0; i < vecRMgr.size(); i++) rMgr->SetBlock(i, vecRMgr[i]);
650 std::vector<std::string> tokens;
655 tokenize(reorder,
" \t\n",
"[]", tokens);
658 Teuchos::RCP<BlockReorderManager> mgr = blockedReorderFromTokens(tokens);
int GetIndex() const
Get the the index that is stored in this block.
std::vector< Teuchos::RCP< BlockReorderManager > > children_
Definitions of the subblocks.
Class that describes how a flat blocked operator should be reordered.
virtual int GetNumBlocks() const
Gets the number of subblocks.
virtual int LargestIndex() const
Largest index in this manager.
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.
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 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...
virtual std::string toString() const
For sanities sake, print a readable string.
BlockReorderManager()
Basic empty constructor.
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.
virtual void SetBlock(int blockIndex, int reorder)
Sets the sublock to a specific index value.