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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #ifndef TPETRA_VBRUTILS_HPP
00030 #define TPETRA_VBRUTILS_HPP
00031 
00032 #include "Teuchos_Array.hpp"
00033 #include "Teuchos_ArrayRCP.hpp"
00034 #include "Tpetra_BlockMap.hpp"
00035 #include <map>
00036 
00041 namespace Tpetra {
00045 namespace VbrUtils {
00046 
00047 template<typename LocalOrdinal, typename Scalar>
00048 struct BlkInfo {
00049   LocalOrdinal numPtRows;
00050   LocalOrdinal numPtCols;
00051   Teuchos::ArrayRCP<Scalar> blkEntry;
00052 };
00053 
00054 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00055 struct VbrData {
00056 
00057   typedef typename std::map<GlobalOrdinal,BlkInfo<LocalOrdinal,Scalar> > RowGlobalCols;
00058 
00059   VbrData()
00060    : row_map(),
00061      data(Teuchos::rcp(new Teuchos::Array<RowGlobalCols>))
00062   {}
00063 
00064   ~VbrData(){}
00065 
00066   std::map<GlobalOrdinal,LocalOrdinal> row_map;
00067   Teuchos::RCP<Teuchos::Array<RowGlobalCols> > data;
00068 };
00069 
00070 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00071 inline
00072 void zeroEntries(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata)
00073 {
00074   typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00075 
00076   typedef typename Teuchos::Array<RowGlobalCols>::size_type Tsize_t;
00077 
00078   for(Tsize_t i=0; i<vbrdata.data->size(); ++i) {
00079     RowGlobalCols& rgc = (*vbrdata.data)[i];
00080     typename RowGlobalCols::iterator
00081       rgc_it = rgc.begin(), rgc_end = rgc.end();
00082     for(; rgc_it != rgc_end; ++rgc_it) {
00083       BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00084       std::fill(blk.blkEntry.begin(), blk.blkEntry.end(), 0.0);
00085     }    
00086   }
00087 }
00088 
00089 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00090 inline
00091 void getGlobalBlockEntryViewNonConst(
00092                    VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00093                    GlobalOrdinal globalBlockRow,
00094                    GlobalOrdinal globalBlockCol,
00095                    LocalOrdinal& numPtRows,
00096                    LocalOrdinal& numPtCols,
00097                    Teuchos::ArrayRCP<Scalar>& blockEntry)
00098 {
00099   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator miter =
00100       vbrdata.row_map.find(globalBlockRow);
00101 
00102   LocalOrdinal localBlockRow = 0;
00103 
00104   if (miter == vbrdata.row_map.end()) {
00105     localBlockRow = vbrdata.data->size();
00106     vbrdata.row_map.insert(std::make_pair(globalBlockRow,localBlockRow));
00107     vbrdata.data->resize(localBlockRow+1);
00108   }
00109   else {
00110     localBlockRow = miter->second;
00111   }
00112 
00113   typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00114   RowGlobalCols& blkrow = (*vbrdata.data)[localBlockRow];
00115   typename RowGlobalCols::iterator iter = blkrow.find(globalBlockCol);
00116   if (iter == blkrow.end()) {
00117     BlkInfo<LocalOrdinal,Scalar> blk;
00118     blk.numPtRows = numPtRows;
00119     blk.numPtCols = numPtCols;
00120     size_t blockSize = numPtRows*numPtCols;
00121     blk.blkEntry = Teuchos::arcp(new Scalar[blockSize], 0, blockSize);
00122     std::fill(blk.blkEntry.begin(), blk.blkEntry.end(), 0);
00123     blkrow.insert(iter, std::make_pair(globalBlockCol, blk));
00124     blockEntry = blk.blkEntry;
00125   }
00126   else {
00127     BlkInfo<LocalOrdinal,Scalar>& blk = iter->second;
00128     numPtRows = blk.numPtRows;
00129     numPtCols = blk.numPtCols;
00130     blockEntry = blk.blkEntry;
00131   }
00132 }
00133 
00134 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar, class Node>
00135 Teuchos::RCP<const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node> >
00136 createOverlapMap(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00137                 const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& rowmap)
00138 {
00139   Teuchos::Array<GlobalOrdinal> importBlocks;
00140   Teuchos::Array<int> src_procs;
00141 
00142   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator
00143     iter = vbrdata.row_map.begin(),
00144     iter_end = vbrdata.row_map.end();
00145   for(; iter != iter_end; ++iter) {
00146     importBlocks.push_back(iter->first);
00147   }
00148 
00149   Teuchos::Array<GlobalOrdinal> importFirstPointInBlocks(importBlocks.size());
00150   Teuchos::Array<LocalOrdinal>  importBlockSizes(importBlocks.size());
00151   rowmap.getRemoteBlockInfo(importBlocks(), importFirstPointInBlocks(), importBlockSizes());
00152 
00153   return Teuchos::rcp(
00154      new Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>(
00155            Teuchos::OrdinalTraits<global_size_t>::invalid(),
00156            importBlocks(),
00157            importFirstPointInBlocks(),
00158            importBlockSizes(),
00159            rowmap.getPointMap()->getIndexBase(),
00160            rowmap.getPointMap()->getComm(),
00161            rowmap.getPointMap()->getNode()
00162          )
00163   );
00164 }
00165 
00166 template<typename LocalOrdinal,typename GlobalOrdinal,typename Scalar,typename Node>
00167 struct VbrDataDist : public Tpetra::DistObject<char,LocalOrdinal,GlobalOrdinal,Node> {
00168   VbrDataDist(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbr_data,
00169               const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& importmap,
00170               const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& rowmap)
00171    : Tpetra::DistObject<char,LocalOrdinal,GlobalOrdinal,Node>(
00172         convertBlockMapToPointMap(importmap)),
00173      vbrdata(vbr_data)
00174   {}
00175 
00176   ~VbrDataDist(){}
00177 
00178   VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata;
00179   Teuchos::RCP<Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node> > importMap;
00180 
00181   bool checkSizes(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source)
00182   {
00183     bool ok = this->getMap()->getMinAllGlobalIndex() <= source.getMap()->getMinAllGlobalIndex();
00184     ok = ok && this->getMap()->getMaxAllGlobalIndex() >= source.getMap()->getMaxAllGlobalIndex();
00185     return ok;
00186   }
00187 
00188   void copyAndPermute(
00189    const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source,
00190    size_t numSameIDs,
00191    const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs,
00192    const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs)
00193   {
00194     throw std::runtime_error("VbrDataDist hasn't implemented copyAndPermute!!");
00195   }
00196 
00197   void packAndPrepare(
00198      const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source,
00199      const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
00200      Teuchos::Array<char>& exports,
00201      const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00202      size_t& constantNumPackets,
00203      Distributor& distor)
00204   {
00205     const VbrDataDist<LocalOrdinal,GlobalOrdinal,Scalar,Node>* vdd = dynamic_cast<const VbrDataDist<LocalOrdinal,GlobalOrdinal,Scalar,Node>*>(&source);
00206     if (vdd == NULL) {
00207       throw std::runtime_error("VbrDataDist::packAndPrepare ERROR, dynamic_cast failed.");
00208     }
00209     typedef typename Teuchos::ArrayView<const LocalOrdinal>::size_type Tsize_t;
00210     typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00211   
00212     //We will pack each row's data into the exports buffer as follows:
00213     //[num-block-cols,numPtRows,{list-of-blk-cols},{list-of-ptColsPerBlkCol},{all vals}]
00214     //so the length of the char exports buffer for a row is:
00215     //sizeof(LocalOrdinal)*(2+2*(num-block-cols)) + sizeof(Scalar)*numPtRows*sum(numPtCols_i)
00216   
00217     //For each row corresponding to exportLIDs, accumulate the size that it will
00218     //occupy in the exports buffer:
00219     size_t total_exports_size = 0;
00220     for(Tsize_t i=0; i<exportLIDs.size(); ++i) {
00221       LocalOrdinal numBlkCols = 0;
00222       LocalOrdinal numScalars = 0;
00223   
00224       LocalOrdinal localIndex = exportLIDs[i];
00225       const RowGlobalCols& rgc = (*vdd->vbrdata.data)[localIndex];
00226       numBlkCols = rgc.size();
00227       typename RowGlobalCols::const_iterator rgc_it = rgc.begin(), rgc_end = rgc.end();
00228       for(; rgc_it!=rgc_end; ++rgc_it) {
00229         const BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00230         numScalars += blk.numPtRows*blk.numPtCols;
00231       }
00232   
00233       size_t size_for_this_row = sizeof(GlobalOrdinal)*(2+2*numBlkCols)
00234                    + sizeof(Scalar)*numScalars;
00235       numPacketsPerLID[i] = size_for_this_row;
00236       total_exports_size += size_for_this_row;
00237     }
00238   
00239     exports.resize(total_exports_size);
00240   
00241     ArrayView<char> avIndsC, avValsC;
00242     ArrayView<Scalar> avVals;
00243   
00244     Teuchos::Array<GlobalOrdinal> blkCols;
00245     Teuchos::Array<LocalOrdinal> ptColsPerBlkCol;
00246     Teuchos::Array<Scalar> blkEntries;
00247   
00248     size_t offset = 0;
00249     for(Tsize_t i=0; i<exportLIDs.size(); ++i) {
00250       blkCols.clear();
00251       ptColsPerBlkCol.clear();
00252       blkEntries.clear();
00253   
00254       LocalOrdinal numPtRows = 0;
00255       LocalOrdinal localIndex = exportLIDs[i];
00256       const RowGlobalCols& rgc = (*vdd->vbrdata.data)[localIndex];
00257       typename RowGlobalCols::const_iterator rgc_it = rgc.begin(), rgc_end = rgc.end();
00258       for(; rgc_it!=rgc_end; ++rgc_it) {
00259         blkCols.push_back(rgc_it->first);
00260         const BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00261         numPtRows = blk.numPtRows;//should be the same for all cols
00262         ptColsPerBlkCol.push_back(blk.numPtCols);
00263         for(LocalOrdinal j=0; j<blk.numPtRows*blk.numPtCols; ++j) {
00264           blkEntries.push_back(blk.blkEntry[j]);
00265         }
00266       }
00267   
00268       LocalOrdinal numBlkCols = blkCols.size();
00269       LocalOrdinal numScalars = blkEntries.size();
00270   
00271       LocalOrdinal num_chars_for_ordinals = (2*numBlkCols+2)*sizeof(GlobalOrdinal);
00272       //get export views
00273       avIndsC = exports(offset, num_chars_for_ordinals);
00274       avValsC = exports(offset+ num_chars_for_ordinals, numScalars*sizeof(Scalar));
00275       ArrayView<GlobalOrdinal> avInds = av_reinterpret_cast<GlobalOrdinal>(avIndsC);
00276       typename ArrayView<GlobalOrdinal>::iterator ind_it = avInds.begin();
00277 
00278       *ind_it++ = numBlkCols;
00279       *ind_it++ = numPtRows;
00280       for(Tsize_t j=0; j<blkCols.size(); ++j) {
00281         *ind_it++ = blkCols[j];
00282       }
00283       for(Tsize_t j=0; j<ptColsPerBlkCol.size(); ++j) {
00284         *ind_it++ = ptColsPerBlkCol[j];
00285       }
00286 
00287       avVals = av_reinterpret_cast<Scalar>(avValsC);
00288       std::copy(blkEntries.begin(), blkEntries.end(), avVals.begin());
00289 
00290       size_t size_for_this_row = sizeof(GlobalOrdinal)*(2+2*numBlkCols)
00291                      + sizeof(Scalar)*numScalars;
00292       offset += size_for_this_row;
00293     }
00294     constantNumPackets = 0;
00295   }
00296 
00297   void unpackAndCombine(
00298       const Teuchos::ArrayView<const LocalOrdinal>& importLIDs,
00299       const Teuchos::ArrayView<const char>& imports,
00300       const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00301       size_t constantNumPackets,
00302       Distributor& distor,
00303       CombineMode CM)
00304   {
00305     throw std::runtime_error("VbrDataDist hasn't implemented unpackAndCombine!!");
00306   }
00307 
00308 };
00309 
00310 }//namespace VbrUtils
00311 }//namespace Tpetra
00312 
00313 #endif
00314 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines