Tpetra Matrix/Vector Services Version of the Day
Tpetra_VbrUtils.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_VBRUTILS_HPP
00043 #define TPETRA_VBRUTILS_HPP
00044 
00045 #include "Teuchos_Array.hpp"
00046 #include "Teuchos_ArrayRCP.hpp"
00047 #include "Tpetra_BlockMap.hpp"
00048 #include <map>
00049 
00054 namespace Tpetra {
00058 namespace VbrUtils {
00059 
00060 template<typename LocalOrdinal, typename Scalar>
00061 struct BlkInfo {
00062   LocalOrdinal numPtRows;
00063   LocalOrdinal numPtCols;
00064   Teuchos::ArrayRCP<Scalar> blkEntry;
00065 };
00066 
00067 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00068 struct VbrData {
00069 
00070   typedef typename std::map<GlobalOrdinal,BlkInfo<LocalOrdinal,Scalar> > RowGlobalCols;
00071 
00072   VbrData()
00073    : row_map(),
00074      data(Teuchos::rcp(new Teuchos::Array<RowGlobalCols>))
00075   {}
00076 
00077   ~VbrData(){}
00078 
00079   std::map<GlobalOrdinal,LocalOrdinal> row_map;
00080   Teuchos::RCP<Teuchos::Array<RowGlobalCols> > data;
00081 };
00082 
00083 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00084 inline
00085 void zeroEntries(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata)
00086 {
00087   typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00088 
00089   typedef typename Teuchos::Array<RowGlobalCols>::size_type Tsize_t;
00090 
00091   for(Tsize_t i=0; i<vbrdata.data->size(); ++i) {
00092     RowGlobalCols& rgc = (*vbrdata.data)[i];
00093     typename RowGlobalCols::iterator
00094       rgc_it = rgc.begin(), rgc_end = rgc.end();
00095     for(; rgc_it != rgc_end; ++rgc_it) {
00096       BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00097       std::fill(blk.blkEntry.begin(), blk.blkEntry.end(), 0.0);
00098     }    
00099   }
00100 }
00101 
00102 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00103 inline
00104 void getGlobalBlockEntryViewNonConst(
00105                    VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00106                    GlobalOrdinal globalBlockRow,
00107                    GlobalOrdinal globalBlockCol,
00108                    LocalOrdinal& numPtRows,
00109                    LocalOrdinal& numPtCols,
00110                    Teuchos::ArrayRCP<Scalar>& blockEntry)
00111 {
00112   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator miter =
00113       vbrdata.row_map.find(globalBlockRow);
00114 
00115   LocalOrdinal localBlockRow = 0;
00116 
00117   if (miter == vbrdata.row_map.end()) {
00118     localBlockRow = vbrdata.data->size();
00119     vbrdata.row_map.insert(std::make_pair(globalBlockRow,localBlockRow));
00120     vbrdata.data->resize(localBlockRow+1);
00121   }
00122   else {
00123     localBlockRow = miter->second;
00124   }
00125 
00126   typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00127   RowGlobalCols& blkrow = (*vbrdata.data)[localBlockRow];
00128   typename RowGlobalCols::iterator iter = blkrow.find(globalBlockCol);
00129   if (iter == blkrow.end()) {
00130     BlkInfo<LocalOrdinal,Scalar> blk;
00131     blk.numPtRows = numPtRows;
00132     blk.numPtCols = numPtCols;
00133     size_t blockSize = numPtRows*numPtCols;
00134     blk.blkEntry = Teuchos::arcp(new Scalar[blockSize], 0, blockSize);
00135     std::fill(blk.blkEntry.begin(), blk.blkEntry.end(), 0);
00136     blkrow.insert(iter, std::make_pair(globalBlockCol, blk));
00137     blockEntry = blk.blkEntry;
00138   }
00139   else {
00140     BlkInfo<LocalOrdinal,Scalar>& blk = iter->second;
00141     numPtRows = blk.numPtRows;
00142     numPtCols = blk.numPtCols;
00143     blockEntry = blk.blkEntry;
00144   }
00145 }
00146 
00147 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar, class Node>
00148 Teuchos::RCP<const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node> >
00149 createOverlapMap(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00150                 const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& rowmap)
00151 {
00152   Teuchos::Array<GlobalOrdinal> importBlocks;
00153   Teuchos::Array<int> src_procs;
00154 
00155   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator
00156     iter = vbrdata.row_map.begin(),
00157     iter_end = vbrdata.row_map.end();
00158   for(; iter != iter_end; ++iter) {
00159     importBlocks.push_back(iter->first);
00160   }
00161 
00162   Teuchos::Array<GlobalOrdinal> importFirstPointInBlocks(importBlocks.size());
00163   Teuchos::Array<LocalOrdinal>  importBlockSizes(importBlocks.size());
00164   rowmap.getRemoteBlockInfo(importBlocks(), importFirstPointInBlocks(), importBlockSizes());
00165 
00166   return Teuchos::rcp(
00167      new Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>(
00168            Teuchos::OrdinalTraits<global_size_t>::invalid(),
00169            importBlocks(),
00170            importFirstPointInBlocks(),
00171            importBlockSizes(),
00172            rowmap.getPointMap()->getIndexBase(),
00173            rowmap.getPointMap()->getComm(),
00174            rowmap.getPointMap()->getNode()
00175          )
00176   );
00177 }
00178 
00179 template<typename LocalOrdinal,typename GlobalOrdinal,typename Scalar,typename Node>
00180 struct VbrDataDist : public Tpetra::DistObject<char,LocalOrdinal,GlobalOrdinal,Node> {
00181   VbrDataDist(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbr_data,
00182               const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& importmap,
00183               const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& rowmap)
00184    : Tpetra::DistObject<char,LocalOrdinal,GlobalOrdinal,Node>(
00185         convertBlockMapToPointMap(importmap)),
00186      vbrdata(vbr_data)
00187   {}
00188 
00189   ~VbrDataDist(){}
00190 
00191   VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata;
00192   Teuchos::RCP<Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node> > importMap;
00193 
00194   bool checkSizes(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source)
00195   {
00196     bool ok = this->getMap()->getMinAllGlobalIndex() <= source.getMap()->getMinAllGlobalIndex();
00197     ok = ok && this->getMap()->getMaxAllGlobalIndex() >= source.getMap()->getMaxAllGlobalIndex();
00198     return ok;
00199   }
00200 
00201   void copyAndPermute(
00202    const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source,
00203    size_t numSameIDs,
00204    const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs,
00205    const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs)
00206   {
00207     throw std::runtime_error("VbrDataDist hasn't implemented copyAndPermute!!");
00208   }
00209 
00210   void packAndPrepare(
00211      const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source,
00212      const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
00213      Teuchos::Array<char>& exports,
00214      const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00215      size_t& constantNumPackets,
00216      Distributor& distor)
00217   {
00218     const VbrDataDist<LocalOrdinal,GlobalOrdinal,Scalar,Node>* vdd = dynamic_cast<const VbrDataDist<LocalOrdinal,GlobalOrdinal,Scalar,Node>*>(&source);
00219     if (vdd == NULL) {
00220       throw std::runtime_error("VbrDataDist::packAndPrepare ERROR, dynamic_cast failed.");
00221     }
00222     typedef typename Teuchos::ArrayView<const LocalOrdinal>::size_type Tsize_t;
00223     typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00224   
00225     //We will pack each row's data into the exports buffer as follows:
00226     //[num-block-cols,numPtRows,{list-of-blk-cols},{list-of-ptColsPerBlkCol},{all vals}]
00227     //so the length of the char exports buffer for a row is:
00228     //sizeof(LocalOrdinal)*(2+2*(num-block-cols)) + sizeof(Scalar)*numPtRows*sum(numPtCols_i)
00229   
00230     //For each row corresponding to exportLIDs, accumulate the size that it will
00231     //occupy in the exports buffer:
00232     size_t total_exports_size = 0;
00233     for(Tsize_t i=0; i<exportLIDs.size(); ++i) {
00234       LocalOrdinal numBlkCols = 0;
00235       LocalOrdinal numScalars = 0;
00236   
00237       LocalOrdinal localIndex = exportLIDs[i];
00238       const RowGlobalCols& rgc = (*vdd->vbrdata.data)[localIndex];
00239       numBlkCols = rgc.size();
00240       typename RowGlobalCols::const_iterator rgc_it = rgc.begin(), rgc_end = rgc.end();
00241       for(; rgc_it!=rgc_end; ++rgc_it) {
00242         const BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00243         numScalars += blk.numPtRows*blk.numPtCols;
00244       }
00245   
00246       size_t size_for_this_row = sizeof(GlobalOrdinal)*(2+2*numBlkCols)
00247                    + sizeof(Scalar)*numScalars;
00248       numPacketsPerLID[i] = size_for_this_row;
00249       total_exports_size += size_for_this_row;
00250     }
00251   
00252     exports.resize(total_exports_size);
00253   
00254     ArrayView<char> avIndsC, avValsC;
00255     ArrayView<Scalar> avVals;
00256   
00257     Teuchos::Array<GlobalOrdinal> blkCols;
00258     Teuchos::Array<LocalOrdinal> ptColsPerBlkCol;
00259     Teuchos::Array<Scalar> blkEntries;
00260   
00261     size_t offset = 0;
00262     for(Tsize_t i=0; i<exportLIDs.size(); ++i) {
00263       blkCols.clear();
00264       ptColsPerBlkCol.clear();
00265       blkEntries.clear();
00266   
00267       LocalOrdinal numPtRows = 0;
00268       LocalOrdinal localIndex = exportLIDs[i];
00269       const RowGlobalCols& rgc = (*vdd->vbrdata.data)[localIndex];
00270       typename RowGlobalCols::const_iterator rgc_it = rgc.begin(), rgc_end = rgc.end();
00271       for(; rgc_it!=rgc_end; ++rgc_it) {
00272         blkCols.push_back(rgc_it->first);
00273         const BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00274         numPtRows = blk.numPtRows;//should be the same for all cols
00275         ptColsPerBlkCol.push_back(blk.numPtCols);
00276         for(LocalOrdinal j=0; j<blk.numPtRows*blk.numPtCols; ++j) {
00277           blkEntries.push_back(blk.blkEntry[j]);
00278         }
00279       }
00280   
00281       LocalOrdinal numBlkCols = blkCols.size();
00282       LocalOrdinal numScalars = blkEntries.size();
00283   
00284       LocalOrdinal num_chars_for_ordinals = (2*numBlkCols+2)*sizeof(GlobalOrdinal);
00285       //get export views
00286       avIndsC = exports(offset, num_chars_for_ordinals);
00287       avValsC = exports(offset+ num_chars_for_ordinals, numScalars*sizeof(Scalar));
00288       ArrayView<GlobalOrdinal> avInds = av_reinterpret_cast<GlobalOrdinal>(avIndsC);
00289       typename ArrayView<GlobalOrdinal>::iterator ind_it = avInds.begin();
00290 
00291       *ind_it++ = numBlkCols;
00292       *ind_it++ = numPtRows;
00293       for(Tsize_t j=0; j<blkCols.size(); ++j) {
00294         *ind_it++ = blkCols[j];
00295       }
00296       for(Tsize_t j=0; j<ptColsPerBlkCol.size(); ++j) {
00297         *ind_it++ = ptColsPerBlkCol[j];
00298       }
00299 
00300       avVals = av_reinterpret_cast<Scalar>(avValsC);
00301       std::copy(blkEntries.begin(), blkEntries.end(), avVals.begin());
00302 
00303       size_t size_for_this_row = sizeof(GlobalOrdinal)*(2+2*numBlkCols)
00304                      + sizeof(Scalar)*numScalars;
00305       offset += size_for_this_row;
00306     }
00307     constantNumPackets = 0;
00308   }
00309 
00310   void unpackAndCombine(
00311       const Teuchos::ArrayView<const LocalOrdinal>& importLIDs,
00312       const Teuchos::ArrayView<const char>& imports,
00313       const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00314       size_t constantNumPackets,
00315       Distributor& distor,
00316       CombineMode CM)
00317   {
00318     throw std::runtime_error("VbrDataDist hasn't implemented unpackAndCombine!!");
00319   }
00320 
00321 };
00322 
00323 }//namespace VbrUtils
00324 }//namespace Tpetra
00325 
00326 #endif
00327 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines