Teko Version of the Day
Teko_BlockedReordering.cpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #include <iostream>
00048 #include <string>
00049 #include <vector>
00050 #include <stack>
00051 
00052 #include "Teko_BlockedReordering.hpp"
00053 #include "Teko_Utilities.hpp"
00054 
00055 #include "Teuchos_VerboseObject.hpp"
00056 #include "Teuchos_RCP.hpp"
00057 #include "Teuchos_StrUtils.hpp"
00058 
00059 #include "Thyra_DefaultProductMultiVector.hpp"
00060 #include "Thyra_DefaultProductVectorSpace.hpp"
00061 
00062 using Teuchos::RCP;
00063 using Teuchos::rcp;
00064 using Teuchos::rcp_dynamic_cast;
00065 using Teuchos::Array;
00066 
00067 namespace Teko {
00068 
00069 void BlockReorderManager::SetBlock(int blockIndex,int reorder)
00070 { 
00071    TEUCHOS_ASSERT(blockIndex<(int) children_.size());
00072 
00073    RCP<BlockReorderManager> child = rcp(new BlockReorderLeaf(reorder));
00074 
00075    children_[blockIndex] = child;
00076 }
00077 
00090 void BlockReorderManager::SetBlock(int blockIndex,const RCP<BlockReorderManager> & reorder)
00091 {
00092    TEUCHOS_ASSERT(blockIndex<(int) children_.size());
00093 
00094    children_[blockIndex] = reorder;
00095 }
00096 
00097 const Teuchos::RCP<BlockReorderManager> BlockReorderManager::GetBlock(int blockIndex)
00098 {
00099    TEUCHOS_ASSERT(blockIndex<(int) children_.size());
00100 
00101    if(children_[blockIndex]==Teuchos::null)
00102       children_[blockIndex] = rcp(new BlockReorderManager());
00103 
00104    return children_[blockIndex];
00105 }
00106 
00107 const Teuchos::RCP<const BlockReorderManager> BlockReorderManager::GetBlock(int blockIndex) const
00108 {
00109    TEUCHOS_ASSERT(blockIndex<(int) children_.size());
00110 
00111    return children_[blockIndex];
00112 }
00113 
00114 std::string BlockReorderManager::toString() const
00115 {
00116    // build the string by recursively calling each child
00117    std::stringstream ss;
00118    ss << "[";
00119    for(unsigned int i=0;i<children_.size();i++) {
00120       if(children_[i]==Teuchos::null) 
00121          ss << " <NULL> ";
00122       else
00123          ss << " " << children_[i]->toString() << " "; 
00124    }
00125    ss << "]";
00126 
00127    return ss.str();
00128 }
00129 
00130 int BlockReorderManager::LargestIndex() const
00131 {
00132    int max = 0;
00133    for(unsigned int i=0;i<children_.size();i++) {
00134       // see if current child is larger 
00135       if(children_[i]!=Teuchos::null) {
00136          int subMax = children_[i]->LargestIndex();
00137          max = max > subMax ? max : subMax;
00138       }
00139    }
00140 
00141    return max;
00142 } 
00143 
00144 Teuchos::RCP<const Thyra::LinearOpBase<double> >
00145 buildReorderedLinearOp(const BlockReorderManager & bmm,
00146                        const Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > & blkOp)
00147 {
00148    return buildReorderedLinearOp(bmm,bmm,blkOp);
00149 }
00150 
00151 Teuchos::RCP<const Thyra::LinearOpBase<double> >
00152 buildReorderedLinearOp(const BlockReorderManager & rowMgr,const BlockReorderManager & colMgr,
00153                        const Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > & blkOp)
00154 {
00155    typedef RCP<const BlockReorderManager> BRMptr;
00156 
00157    int rowSz = rowMgr.GetNumBlocks();
00158    int colSz = colMgr.GetNumBlocks();
00159 
00160    if(rowSz==0 && colSz==0) {
00161       // both are leaf nodes
00162       const BlockReorderLeaf & rowLeaf = dynamic_cast<const BlockReorderLeaf &>(rowMgr);
00163       const BlockReorderLeaf & colLeaf = dynamic_cast<const BlockReorderLeaf &>(colMgr);
00164 
00165       // simply return entry in matrix
00166       Teko::LinearOp linOp = blkOp->getBlock(rowLeaf.GetIndex(),colLeaf.GetIndex());
00167 
00168       // somehow we need to set this operator up
00169       if(linOp==Teuchos::null) {
00170          linOp = Thyra::zero(blkOp->productRange()->getBlock(rowLeaf.GetIndex()),
00171                              blkOp->productDomain()->getBlock(colLeaf.GetIndex()));
00172       }
00173 
00174       return linOp;
00175    }
00176    else if(rowSz==0) {
00177       Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
00178 
00179       // operator will be rowSz by colSz
00180       reBlkOp->beginBlockFill(1,colSz);   
00181 
00182       // fill the column entries
00183       for(int col=0;col<colSz;col++) {
00184          BRMptr colPtr = colMgr.GetBlock(col);
00185 
00186          reBlkOp->setBlock(0,col,buildReorderedLinearOp(rowMgr,*colPtr,blkOp));
00187       } 
00188 
00189       // done building
00190       reBlkOp->endBlockFill();   
00191 
00192       return reBlkOp;
00193    }
00194    else if(colSz==0) {
00195       Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
00196 
00197       // operator will be rowSz by colSz
00198       reBlkOp->beginBlockFill(rowSz,1);   
00199 
00200       // fill the row entries
00201       for(int row=0;row<rowSz;row++) {
00202          BRMptr rowPtr = rowMgr.GetBlock(row);
00203 
00204          reBlkOp->setBlock(row,0,buildReorderedLinearOp(*rowPtr,colMgr,blkOp));
00205       } 
00206 
00207       // done building
00208       reBlkOp->endBlockFill();   
00209 
00210       return reBlkOp;
00211    }
00212    else {
00213       Teko::BlockedLinearOp reBlkOp = Teko::createBlockedOp();
00214   
00215       // this is the general case
00216       TEUCHOS_ASSERT(rowSz>0);
00217       TEUCHOS_ASSERT(colSz>0);
00218 
00219       // operator will be rowSz by colSz
00220       reBlkOp->beginBlockFill(rowSz,colSz);   
00221    
00222       for(int row=0;row<rowSz;row++) {
00223          BRMptr rowPtr = rowMgr.GetBlock(row);
00224    
00225          for(int col=0;col<colSz;col++) {
00226             BRMptr colPtr = colMgr.GetBlock(col);
00227    
00228             reBlkOp->setBlock(row,col,buildReorderedLinearOp(*rowPtr,*colPtr,blkOp));
00229          } 
00230       }
00231 
00232       // done building
00233       reBlkOp->endBlockFill();   
00234 
00235       return reBlkOp;
00236    }
00237 }
00238 
00260 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
00261 buildReorderedVectorSpace(const BlockReorderManager & mgr,
00262                           const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<double> > & blkSpc)
00263 {
00264    typedef RCP<const BlockReorderManager> BRMptr;
00265 
00266    int sz = mgr.GetNumBlocks();
00267 
00268    if(sz==0) {
00269       // its a  leaf nodes
00270       const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
00271 
00272       // simply return entry in matrix
00273       return blkSpc->getBlock(leaf.GetIndex());
00274    } 
00275    else {
00276       Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
00277 
00278       // loop over each row
00279       for(int i=0;i<sz;i++) {
00280          BRMptr blkMgr = mgr.GetBlock(i);
00281 
00282          const RCP<const Thyra::VectorSpaceBase<double> > lvs = buildReorderedVectorSpace(*blkMgr,blkSpc);
00283 
00284          vecSpaces.push_back(lvs);
00285       }
00286 
00287       // build a vector space
00288       const RCP<const Thyra::DefaultProductVectorSpace<double> > vs 
00289             = Thyra::productVectorSpace<double>(vecSpaces);
00290 
00291       // build the vector
00292       return vs;
00293    }
00294 }
00295 
00300 Teuchos::RCP<Thyra::MultiVectorBase<double> >
00301 buildReorderedMultiVector(const BlockReorderManager & mgr,
00302                           const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > & blkVec)
00303 {
00304    using Teuchos::rcp_const_cast;
00305 
00306    // give vector const so that the right function is called
00307    const Teuchos::RCP<const Thyra::ProductMultiVectorBase<double> > blkVecConst
00308       = rcp_const_cast<const Thyra::ProductMultiVectorBase<double> >(blkVec);
00309    
00310    // its not really const, so take it away
00311    const Teuchos::RCP<Thyra::MultiVectorBase<double> > result
00312          = rcp_const_cast<Thyra::MultiVectorBase<double> >(buildReorderedMultiVector(mgr,blkVecConst));
00313 
00314    return result; 
00315 }
00316 
00321 Teuchos::RCP<const Thyra::MultiVectorBase<double> >
00322 buildReorderedMultiVector(const BlockReorderManager & mgr,
00323                           const Teuchos::RCP<const Thyra::ProductMultiVectorBase<double> > & blkVec)
00324 {
00325    typedef RCP<const BlockReorderManager> BRMptr;
00326 
00327    int sz = mgr.GetNumBlocks();
00328 
00329    if(sz==0) {
00330       // its a  leaf nodes
00331       const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
00332 
00333       // simply return entry in matrix
00334       return blkVec->getMultiVectorBlock(leaf.GetIndex());
00335    } 
00336    else {
00337       Array<RCP<const Thyra::MultiVectorBase<double> > > multiVecs;
00338       Array<RCP<const Thyra::VectorSpaceBase<double> > > vecSpaces;
00339 
00340       // loop over each row
00341       for(int i=0;i<sz;i++) {
00342          BRMptr blkMgr = mgr.GetBlock(i);
00343 
00344          const RCP<const Thyra::MultiVectorBase<double> > lmv = buildReorderedMultiVector(*blkMgr,blkVec);
00345          const RCP<const Thyra::VectorSpaceBase<double> > lvs = lmv->range();
00346 
00347          multiVecs.push_back(lmv);
00348          vecSpaces.push_back(lvs);
00349       }
00350 
00351       // build a vector space
00352       const RCP<const Thyra::DefaultProductVectorSpace<double> > vs 
00353             = Thyra::productVectorSpace<double>(vecSpaces);
00354 
00355       // build the vector
00356       return Thyra::defaultProductMultiVector<double>(vs,multiVecs);
00357    }
00358 }
00359 
00363 void buildNonconstFlatMultiVector(const BlockReorderManager & mgr,
00364                                   const RCP<Thyra::MultiVectorBase<double> > & blkVec,
00365                                   Array<RCP<Thyra::MultiVectorBase<double> > > & multivecs,
00366                                   Array<RCP<const Thyra::VectorSpaceBase<double> > > & vecspaces)
00367 {
00368    typedef RCP<const BlockReorderManager> BRMptr;
00369 
00370    int sz = mgr.GetNumBlocks();
00371 
00372    if(sz==0) {
00373       // its a  leaf nodes
00374       const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
00375       int index = leaf.GetIndex();
00376 
00377       // simply return entry in matrix
00378       multivecs[index] = blkVec;
00379       vecspaces[index] = blkVec->range();
00380    } 
00381    else {
00382       const RCP<Thyra::ProductMultiVectorBase<double> > prodMV
00383             = rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(blkVec);
00384  
00385       // get flattened elements from each child
00386       for(int i=0;i<sz;i++) {
00387          const RCP<Thyra::MultiVectorBase<double> > mv = prodMV->getNonconstMultiVectorBlock(i);
00388          buildNonconstFlatMultiVector(*mgr.GetBlock(i),mv,multivecs,vecspaces);
00389       }
00390    }
00391    
00392 }
00393 
00397 void buildFlatMultiVector(const BlockReorderManager & mgr,
00398                           const RCP<const Thyra::MultiVectorBase<double> > & blkVec,
00399                           Array<RCP<const Thyra::MultiVectorBase<double> > > & multivecs,
00400                           Array<RCP<const Thyra::VectorSpaceBase<double> > > & vecspaces)
00401 {
00402    typedef RCP<const BlockReorderManager> BRMptr;
00403 
00404    int sz = mgr.GetNumBlocks();
00405 
00406    if(sz==0) {
00407       // its a  leaf nodes
00408       const BlockReorderLeaf & leaf = dynamic_cast<const BlockReorderLeaf &>(mgr);
00409       int index = leaf.GetIndex();
00410 
00411       // simply return entry in matrix
00412       multivecs[index] = blkVec;
00413       vecspaces[index] = blkVec->range();
00414    } 
00415    else {
00416       const RCP<const Thyra::ProductMultiVectorBase<double> > prodMV
00417             = rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<double> >(blkVec);
00418  
00419       // get flattened elements from each child
00420       for(int i=0;i<sz;i++) {
00421          RCP<const Thyra::MultiVectorBase<double> > mv = prodMV->getMultiVectorBlock(i);
00422          buildFlatMultiVector(*mgr.GetBlock(i),mv,multivecs,vecspaces);
00423       }
00424    }
00425    
00426 }
00427 
00432 Teuchos::RCP<Thyra::MultiVectorBase<double> >
00433 buildFlatMultiVector(const BlockReorderManager & mgr,
00434                      const Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > & blkVec)
00435 {
00436    int numBlocks = mgr.LargestIndex()+1;
00437  
00438    Array<RCP<Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
00439    Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
00440 
00441    // flatten everything into a vector first
00442    buildNonconstFlatMultiVector(mgr,blkVec,multivecs,vecspaces);
00443 
00444    // build a vector space
00445    const RCP<Thyra::DefaultProductVectorSpace<double> > vs 
00446          = Thyra::productVectorSpace<double>(vecspaces);
00447 
00448    // build the vector
00449    return Thyra::defaultProductMultiVector<double>(vs,multivecs);
00450 }
00451 
00456 Teuchos::RCP<const Thyra::MultiVectorBase<double> >
00457 buildFlatMultiVector(const BlockReorderManager & mgr,
00458                      const Teuchos::RCP<const Thyra::ProductMultiVectorBase<double> > & blkVec)
00459 {
00460    int numBlocks = mgr.LargestIndex()+1;
00461  
00462    Array<RCP<const Thyra::MultiVectorBase<double> > > multivecs(numBlocks);
00463    Array<RCP<const Thyra::VectorSpaceBase<double> > > vecspaces(numBlocks);
00464 
00465    // flatten everything into a vector first
00466    buildFlatMultiVector(mgr,blkVec,multivecs,vecspaces);
00467 
00468    // build a vector space
00469    const RCP<const Thyra::DefaultProductVectorSpace<double> > vs 
00470          = Thyra::productVectorSpace<double>(vecspaces);
00471 
00472    // build the vector
00473    return Thyra::defaultProductMultiVector<double>(vs,multivecs);
00474 }
00475 
00477 // The next three functions are useful for parsing the string
00478 // description of a BlockReorderManager.
00480 
00481 // this function tokenizes a string, breaking out whitespace but saving the
00482 // brackets [,] as special tokens.
00483 static void tokenize(std::string srcInput,std::string whitespace,std::string prefer,
00484               std::vector<std::string> & tokens)
00485 {
00486    std::string input = srcInput;
00487    std::vector<std::string> wsTokens;
00488    std::size_t endPos   = input.length()-1;
00489    while(endPos<input.length()) {
00490       std::size_t next = input.find_first_of(whitespace);
00491 
00492 
00493       // get the sub string
00494       std::string s;
00495       if(next!=std::string::npos) {
00496          s = input.substr(0,next);
00497 
00498          // break out the old substring
00499          input = input.substr(next+1,endPos);
00500       }
00501       else {
00502          s = input;
00503          input = "";
00504       }
00505 
00506       endPos   = input.length()-1;
00507 
00508       // add it to the WS tokens list
00509       if(s=="") continue;
00510       wsTokens.push_back(s);
00511    }
00512 
00513    for(unsigned int i=0;i<wsTokens.size();i++) {
00514       // get string to break up
00515       input = wsTokens[i];
00516 
00517       std::size_t endPos   = input.length()-1;
00518       while(endPos<input.length()) {
00519          std::size_t next = input.find_first_of(prefer);
00520 
00521          std::string s = input;
00522          if(next>0 && next<input.length()) {
00523 
00524             // get the sub string
00525             s = input.substr(0,next);
00526  
00527             input = input.substr(next,endPos);
00528          }
00529          else if(next==0) {
00530             // get the sub string
00531             s = input.substr(0,next+1);
00532  
00533             input = input.substr(next+1,endPos);
00534          }
00535          else input = "";
00536 
00537          // break out the old substring
00538          endPos   = input.length()-1;
00539 
00540          // add it to the tokens list
00541          tokens.push_back(s);
00542       }
00543    }
00544 }
00545 
00546 // this function takes a set of tokens and returns the first "block", i.e. those
00547 // values (including) brackets that correspond to the first block
00548 static std::vector<std::string>::const_iterator 
00549 buildSubBlock(std::vector<std::string>::const_iterator begin,
00550               std::vector<std::string>::const_iterator end, std::vector<std::string> & subBlock)
00551 {
00552    std::stack<std::string> matched;
00553    std::vector<std::string>::const_iterator itr;
00554    for(itr=begin;itr!=end;++itr) { 
00555 
00556       subBlock.push_back(*itr);
00557 
00558       // push/pop brackets as they are discovered
00559       if(*itr=="[") matched.push("["); 
00560       else if(*itr=="]") matched.pop(); 
00561 
00562       // found all matching brackets
00563       if(matched.empty()) 
00564          return itr;
00565    }
00566 
00567    TEUCHOS_ASSERT(matched.empty());
00568 
00569    return itr-1;
00570 }
00571 
00572 // This function takes a tokenized vector and converts it to a block reorder manager
00573 static RCP<BlockReorderManager> blockedReorderFromTokens(const std::vector<std::string> & tokens)
00574 {
00575    // base case
00576    if(tokens.size()==1)
00577       return rcp(new Teko::BlockReorderLeaf(Teuchos::StrUtils::atoi(tokens[0])));
00578 
00579    // check first and last character
00580    TEUCHOS_ASSERT(*(tokens.begin())=="[")
00581    TEUCHOS_ASSERT(*(tokens.end()-1)=="]");
00582 
00583    std::vector<RCP<Teko::BlockReorderManager> > vecRMgr;
00584    std::vector<std::string>::const_iterator itr = tokens.begin()+1; 
00585    while(itr!=tokens.end()-1) {
00586       // figure out which tokens are relevant for this block
00587       std::vector<std::string> subBlock;
00588       itr = buildSubBlock(itr,tokens.end()-1,subBlock); 
00589 
00590       // build the child block reorder manager
00591       vecRMgr.push_back(blockedReorderFromTokens(subBlock));
00592 
00593       // move the iterator one more
00594       itr++;
00595    }
00596 
00597    // build the parent reorder manager
00598    RCP<Teko::BlockReorderManager> rMgr = rcp(new Teko::BlockReorderManager(vecRMgr.size()));
00599    for(unsigned int i=0;i<vecRMgr.size();i++)
00600       rMgr->SetBlock(i,vecRMgr[i]);
00601 
00602    return rMgr;
00603 }
00604 
00606 
00618 Teuchos::RCP<const BlockReorderManager> blockedReorderFromString(std::string & reorder)
00619 {
00620    // vector of tokens to use
00621    std::vector<std::string> tokens;
00622    
00623    // manager to be returned
00624 
00625    // build tokens vector
00626    tokenize(reorder," \t\n","[]",tokens);
00627 
00628    // parse recursively and build reorder manager
00629    Teuchos::RCP<BlockReorderManager> mgr = blockedReorderFromTokens(tokens);
00630 
00631    return mgr;
00632 }
00633 
00634 } // end namespace Teko
 All Classes Files Functions Variables