Teko Version of the Day
Teko_InterlacedEpetra.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 "Teko_InterlacedEpetra.hpp"
00048 
00049 #include <vector>
00050 
00051 using Teuchos::RCP;
00052 using Teuchos::rcp;
00053 
00054 namespace Teko {
00055 namespace Epetra {
00056 namespace Strided {
00057 
00058 // this assumes that there are numGlobals with numVars each interlaced
00059 // i.e. for numVars = 2 (u,v) then the vector is
00060 //    [u_0,v_0,u_1,v_1,u_2,v_2, ..., u_(numGlobals-1),v_(numGlobals-1)]
00061 void buildSubMaps(int numGlobals,int numVars,const Epetra_Comm & comm,std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps)
00062 {
00063    std::vector<int> vars;
00064    
00065    // build vector describing the sub maps
00066    for(int i=0;i<numVars;i++) vars.push_back(1);
00067 
00068    // build all the submaps
00069    buildSubMaps(numGlobals,vars,comm,subMaps);
00070 }
00071 
00072 // build maps to make other conversions
00073 void buildSubMaps(const Epetra_Map & globalMap,const std::vector<int> & vars,const Epetra_Comm & comm,
00074                   std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
00075 {
00076    buildSubMaps(globalMap.NumGlobalElements(),globalMap.NumMyElements(),globalMap.MinMyGID(),
00077                 vars,comm,subMaps);
00078 }
00079 
00080 // build maps to make other conversions
00081 void buildSubMaps(int numGlobals,const std::vector<int> & vars,const Epetra_Comm & comm,std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
00082 {
00083    std::vector<int>::const_iterator varItr;
00084 
00085    // compute total number of variables
00086    int numGlobalVars = 0;
00087    for(varItr=vars.begin();varItr!=vars.end();++varItr)
00088       numGlobalVars += *varItr;
00089 
00090    // must be an even number of globals
00091    TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
00092 
00093    Epetra_Map sampleMap(numGlobals/numGlobalVars,0,comm);
00094 
00095    buildSubMaps(numGlobals,numGlobalVars*sampleMap.NumMyElements(),numGlobalVars*sampleMap.MinMyGID(),vars,comm,subMaps);
00096 }
00097 
00098 // build maps to make other conversions
00099 void buildSubMaps(int numGlobals,int numMyElements,int minMyGID,const std::vector<int> & vars,const Epetra_Comm & comm,
00100                   std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
00101 {
00102    std::vector<int>::const_iterator varItr;
00103 
00104    // compute total number of variables
00105    int numGlobalVars = 0;
00106    for(varItr=vars.begin();varItr!=vars.end();++varItr)
00107       numGlobalVars += *varItr;
00108 
00109    // must be an even number of globals
00110    TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
00111    TEUCHOS_ASSERT((numMyElements%numGlobalVars)==0);
00112    TEUCHOS_ASSERT((minMyGID%numGlobalVars)==0);
00113 
00114    int numBlocks  = numMyElements/numGlobalVars;
00115    int minBlockID = minMyGID/numGlobalVars;
00116 
00117    subMaps.clear();
00118 
00119    // index into local block in strided map
00120    int blockOffset = 0;
00121    for(varItr=vars.begin();varItr!=vars.end();++varItr) {
00122       int numLocalVars = *varItr;
00123       int numAllElmts = numLocalVars*numGlobals/numGlobalVars;
00124       int numMyElmts = numLocalVars * numBlocks;
00125 
00126       // create global arrays describing the as of yet uncreated maps
00127       std::vector<int> subGlobals;
00128       std::vector<int> contigGlobals; // the contiguous globals
00129 
00130       // loop over each block of variables
00131       int count = 0;
00132       for(int blockNum=0;blockNum<numBlocks;blockNum++) {
00133 
00134          // loop over each local variable in the block
00135          for(int local=0;local<numLocalVars;++local) {
00136             // global block number = minGID+blockNum 
00137             // block begin global id = numGlobalVars*(minGID+blockNum)
00138             // global id block offset = blockOffset+local
00139             subGlobals.push_back((minBlockID+blockNum)*numGlobalVars+blockOffset+local);
00140 
00141             // also build the contiguous IDs
00142             contigGlobals.push_back(numLocalVars*minBlockID+count);
00143             count++;
00144          }
00145       }
00146 
00147       // sanity check
00148       assert((unsigned int) numMyElmts==subGlobals.size());
00149 
00150       // create the map with contiguous elements and the map with global elements
00151       RCP<Epetra_Map> subMap = rcp(new Epetra_Map(numAllElmts,numMyElmts,&subGlobals[0],0,comm));
00152       RCP<Epetra_Map> contigMap = rcp(new Epetra_Map(numAllElmts,numMyElmts,&contigGlobals[0],0,comm));
00153 
00154       Teuchos::set_extra_data(contigMap,"contigMap",Teuchos::inOutArg(subMap));
00155       subMaps.push_back(std::make_pair(numLocalVars,subMap));
00156 
00157       // update the block offset
00158       blockOffset += numLocalVars;
00159    }
00160 }
00161 
00162 void buildExportImport(const Epetra_Map & baseMap, const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,
00163                        std::vector<RCP<Epetra_Export> > & subExport,
00164                        std::vector<RCP<Epetra_Import> > & subImport)
00165 {
00166    std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
00167 
00168    // build importers and exporters
00169    for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
00170       // exctract basic map
00171       const Epetra_Map & map = *(mapItr->second);
00172 
00173       // add new elements to vectors
00174       subImport.push_back(rcp(new Epetra_Import(map,baseMap)));
00175       subExport.push_back(rcp(new Epetra_Export(map,baseMap)));
00176    }
00177 }
00178 
00179 void buildSubVectors(const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,std::vector<RCP<Epetra_MultiVector> > & subVectors,int count)
00180 {
00181    std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
00182 
00183    // build vectors
00184    for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
00185       // exctract basic map
00186       const Epetra_Map & map = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(mapItr->second,"contigMap"));
00187 
00188       // add new elements to vectors
00189       RCP<Epetra_MultiVector> mv = rcp(new Epetra_MultiVector(map,count));
00190       Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(mv)); 
00191       subVectors.push_back(mv);
00192    }
00193 }
00194 
00195 void associateSubVectors(const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,std::vector<RCP<const Epetra_MultiVector> > & subVectors)
00196 {
00197    std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
00198    std::vector<RCP<const Epetra_MultiVector> >::iterator vecItr;
00199 
00200    TEUCHOS_ASSERT(subMaps.size()==subVectors.size());
00201 
00202    // associate the sub vectors with the subMaps
00203    for(mapItr=subMaps.begin(),vecItr=subVectors.begin();mapItr!=subMaps.end();++mapItr,++vecItr)
00204       Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(*vecItr)); 
00205 }
00206 
00207 // build a single subblock Epetra_CrsMatrix
00208 RCP<Epetra_CrsMatrix> buildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps)
00209 {
00210    // get the number of variables families
00211    int numVarFamily = subMaps.size();
00212 
00213    TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
00214    TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
00215 
00216    const Epetra_Map & gRowMap = *subMaps[i].second;
00217    const Epetra_Map & rowMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[i].second,"contigMap");
00218    const Epetra_Map & colMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[j].second,"contigMap");
00219    int colFamilyCnt = subMaps[j].first;
00220 
00221    // compute the number of global variables
00222    // and the row and column block offset
00223    int numGlobalVars = 0;
00224    int rowBlockOffset = 0;
00225    int colBlockOffset = 0;
00226    for(int k=0;k<numVarFamily;k++) {
00227       numGlobalVars += subMaps[k].first;
00228  
00229       // compute block offsets
00230       if(k<i) rowBlockOffset += subMaps[k].first;
00231       if(k<j) colBlockOffset += subMaps[k].first;
00232    }
00233 
00234    // copy all global rows to here
00235    Epetra_Import import(gRowMap,A.RowMap());
00236    Epetra_CrsMatrix localA(Copy,gRowMap,0);
00237    localA.Import(A,import,Insert);
00238 
00239    RCP<Epetra_CrsMatrix> mat = rcp(new Epetra_CrsMatrix(Copy,rowMap,0));
00240 
00241    // get entry information
00242    int numMyRows = rowMap.NumMyElements();
00243    int maxNumEntries = A.GlobalMaxNumEntries();
00244 
00245    // for extraction
00246    std::vector<int> indices(maxNumEntries);
00247    std::vector<double> values(maxNumEntries);
00248 
00249    // for insertion
00250    std::vector<int> colIndices(maxNumEntries);
00251    std::vector<double> colValues(maxNumEntries);
00252 
00253    // insert each row into subblock
00254    // let FillComplete handle column distribution
00255    for(int localRow=0;localRow<numMyRows;localRow++) {
00256       int numEntries = -1; 
00257       int globalRow = gRowMap.GID(localRow);
00258       int contigRow = rowMap.GID(localRow);
00259 
00260       TEUCHOS_ASSERT(globalRow>=0);
00261       TEUCHOS_ASSERT(contigRow>=0);
00262 
00263       // extract a global row copy
00264       int err = localA.ExtractGlobalRowCopy(globalRow, maxNumEntries, numEntries, &values[0], &indices[0]);
00265       TEUCHOS_ASSERT(err==0);
00266 
00267       int numOwnedCols = 0;
00268       for(int localCol=0;localCol<numEntries;localCol++) {
00269          int globalCol = indices[localCol];
00270 
00271          // determinate which block this column ID is in
00272          int block = globalCol / numGlobalVars;
00273          
00274          bool inFamily = true; 
00275  
00276          // test the beginning of the block
00277          inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
00278          inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
00279 
00280          // is this column in the variable family
00281          if(inFamily) {
00282             int familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
00283 
00284             colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
00285             colValues[numOwnedCols] = values[localCol];
00286 
00287             numOwnedCols++;
00288          }
00289       }
00290 
00291       // insert it into the new matrix
00292       mat->InsertGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
00293    }
00294 
00295    // fill it and automagically optimize the storage
00296    mat->FillComplete(colMap,rowMap);
00297 
00298    return mat;
00299 }
00300 
00301 // rebuild a single subblock Epetra_CrsMatrix
00302 void rebuildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,Epetra_CrsMatrix & mat)
00303 {
00304    // get the number of variables families
00305    int numVarFamily = subMaps.size();
00306 
00307    TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
00308    TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
00309    TEUCHOS_ASSERT(mat.Filled());
00310 
00311    const Epetra_Map & gRowMap = *subMaps[i].second;
00312    const Epetra_Map & rowMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[i].second,"contigMap");
00313    int colFamilyCnt = subMaps[j].first;
00314 
00315    // compute the number of global variables
00316    // and the row and column block offset
00317    int numGlobalVars = 0;
00318    int rowBlockOffset = 0;
00319    int colBlockOffset = 0;
00320    for(int k=0;k<numVarFamily;k++) {
00321       numGlobalVars += subMaps[k].first;
00322  
00323       // compute block offsets
00324       if(k<i) rowBlockOffset += subMaps[k].first;
00325       if(k<j) colBlockOffset += subMaps[k].first;
00326    }
00327 
00328    // copy all global rows to here
00329    Epetra_Import import(gRowMap,A.RowMap());
00330    Epetra_CrsMatrix localA(Copy,gRowMap,0);
00331    localA.Import(A,import,Insert);
00332 
00333    // clear out the old matrix
00334    mat.PutScalar(0.0);
00335 
00336    // get entry information
00337    int numMyRows = rowMap.NumMyElements();
00338    int maxNumEntries = A.GlobalMaxNumEntries();
00339 
00340    // for extraction
00341    std::vector<int> indices(maxNumEntries);
00342    std::vector<double> values(maxNumEntries);
00343 
00344    // for insertion
00345    std::vector<int> colIndices(maxNumEntries);
00346    std::vector<double> colValues(maxNumEntries);
00347 
00348    // insert each row into subblock
00349    // let FillComplete handle column distribution
00350    for(int localRow=0;localRow<numMyRows;localRow++) {
00351       int numEntries = -1; 
00352       int globalRow = gRowMap.GID(localRow);
00353       int contigRow = rowMap.GID(localRow);
00354 
00355       TEUCHOS_ASSERT(globalRow>=0);
00356       TEUCHOS_ASSERT(contigRow>=0);
00357 
00358       // extract a global row copy
00359       int err = localA.ExtractGlobalRowCopy(globalRow, maxNumEntries, numEntries, &values[0], &indices[0]);
00360       TEUCHOS_ASSERT(err==0);
00361 
00362       int numOwnedCols = 0;
00363       for(int localCol=0;localCol<numEntries;localCol++) {
00364          int globalCol = indices[localCol];
00365 
00366          // determinate which block this column ID is in
00367          int block = globalCol / numGlobalVars;
00368          
00369          bool inFamily = true; 
00370  
00371          // test the beginning of the block
00372          inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
00373          inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
00374 
00375          // is this column in the variable family
00376          if(inFamily) {
00377             int familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
00378 
00379             colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
00380             colValues[numOwnedCols] = values[localCol];
00381 
00382             numOwnedCols++;
00383          }
00384       }
00385 
00386       // insert it into the new matrix
00387       mat.SumIntoGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
00388    }
00389 }
00390 
00391 
00392 // collect subvectors into a single global vector
00393 void many2one(Epetra_MultiVector & one, const std::vector<RCP<const Epetra_MultiVector> > & many,
00394                                    const std::vector<RCP<Epetra_Export> > & subExport)
00395 {
00396    // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
00397    std::vector<RCP<const Epetra_MultiVector> >::const_iterator vecItr;
00398    std::vector<RCP<Epetra_Export> >::const_iterator expItr;
00399 
00400    // using Exporters fill the empty vector from the sub-vectors
00401    for(vecItr=many.begin(),expItr=subExport.begin();
00402        vecItr!=many.end();++vecItr,++expItr) {
00403 
00404       // for ease of access to the source
00405       RCP<const Epetra_MultiVector> srcVec = *vecItr;
00406 
00407       // extract the map with global indicies from the current vector
00408       const Epetra_Map & globalMap = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(srcVec,"globalMap"));
00409 
00410       // build the export vector as a view of the destination
00411       Epetra_MultiVector exportVector(View,globalMap,srcVec->Values(),srcVec->Stride(),srcVec->NumVectors());
00412       one.Export(exportVector,**expItr,Insert);
00413    }
00414 }
00415 
00416 // distribute one global vector into a many subvectors
00417 void one2many(std::vector<RCP<Epetra_MultiVector> > & many,const Epetra_MultiVector & single,
00418                                    const std::vector<RCP<Epetra_Import> > & subImport)
00419 {
00420    // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
00421    std::vector<RCP<Epetra_MultiVector> >::const_iterator vecItr;
00422    std::vector<RCP<Epetra_Import> >::const_iterator impItr;
00423 
00424    // using Importers fill the sub vectors from the mama vector
00425    for(vecItr=many.begin(),impItr=subImport.begin();
00426        vecItr!=many.end();++vecItr,++impItr) {
00427       // for ease of access to the destination
00428       RCP<Epetra_MultiVector> destVec = *vecItr;
00429 
00430       // extract the map with global indicies from the current vector
00431       const Epetra_Map & globalMap = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(destVec,"globalMap"));
00432 
00433       // build the import vector as a view on the destination
00434       Epetra_MultiVector importVector(View,globalMap,destVec->Values(),destVec->Stride(),destVec->NumVectors());
00435 
00436       // perform the import
00437       importVector.Import(single,**impItr,Insert);
00438    }
00439 }
00440 
00441 }
00442 } // end namespace Epetra 
00443 } // end namespace Teko
 All Classes Files Functions Variables