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 "Tpetra_SrcDistObject.hpp"
00049 #include "Tpetra_Packable.hpp"
00050 #include <map>
00051 
00056 namespace Tpetra {
00057 
00058   // Forward declaration
00059   class Distributor;
00060   
00061 
00065 namespace VbrUtils {
00066 
00067 template<typename LocalOrdinal, typename Scalar>
00068 struct BlkInfo {
00069   LocalOrdinal numPtRows;
00070   LocalOrdinal numPtCols;
00071   Teuchos::ArrayRCP<Scalar> blkEntry;
00072 };
00073 
00074 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00075 struct VbrData {
00076 
00077   typedef typename std::map<GlobalOrdinal,BlkInfo<LocalOrdinal,Scalar> > RowGlobalCols;
00078 
00079   VbrData() {}
00080   ~VbrData() {}
00081 
00082   void zeroEntries () {
00083     typedef typename Teuchos::Array<RowGlobalCols>::size_type size_type;
00084     typedef typename Teuchos::ScalarTraits<Scalar> STS;
00085 
00086     const size_type numEntries = this->data.size ();
00087     for (size_type i = 0; i < numEntries; ++i) {
00088       RowGlobalCols& rgc = (this->data)[i];
00089       typename RowGlobalCols::iterator rgc_it = rgc.begin(), rgc_end = rgc.end();
00090       for ( ; rgc_it != rgc_end; ++rgc_it) {
00091   BlkInfo<LocalOrdinal,Scalar>& blk = rgc_it->second;
00092   std::fill (blk.blkEntry.begin (), blk.blkEntry.end (), STS::zero ());
00093       }
00094     }
00095   }
00096 
00097   std::map<GlobalOrdinal,LocalOrdinal> row_map;
00098   Teuchos::Array<RowGlobalCols> data;
00099 };
00100 
00101 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar>
00102 void
00103 getGlobalBlockEntryViewNonConst (VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00104          GlobalOrdinal globalBlockRow,
00105          GlobalOrdinal globalBlockCol,
00106          LocalOrdinal& numPtRows,
00107          LocalOrdinal& numPtCols,
00108          Teuchos::ArrayRCP<Scalar>& blockEntry)
00109 {
00110   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator miter =
00111       vbrdata.row_map.find(globalBlockRow);
00112 
00113   LocalOrdinal localBlockRow = 0;
00114 
00115   if (miter == vbrdata.row_map.end()) {
00116     localBlockRow = vbrdata.data.size ();
00117     vbrdata.row_map.insert(std::make_pair(globalBlockRow,localBlockRow));
00118     vbrdata.data.resize (localBlockRow+1);
00119   }
00120   else {
00121     localBlockRow = miter->second;
00122   }
00123 
00124   typedef typename VbrData<LocalOrdinal,GlobalOrdinal,Scalar>::RowGlobalCols RowGlobalCols;
00125   RowGlobalCols& blkrow = (vbrdata.data)[localBlockRow];
00126   typename RowGlobalCols::iterator iter = blkrow.find(globalBlockCol);
00127   if (iter == blkrow.end()) {
00128     BlkInfo<LocalOrdinal,Scalar> blk;
00129     blk.numPtRows = numPtRows;
00130     blk.numPtCols = numPtCols;
00131     size_t blockSize = numPtRows*numPtCols;
00132     blk.blkEntry = Teuchos::arcp(new Scalar[blockSize], 0, blockSize);
00133     std::fill(blk.blkEntry.begin(), blk.blkEntry.end(), 0);
00134     blkrow.insert(iter, std::make_pair(globalBlockCol, blk));
00135     blockEntry = blk.blkEntry;
00136   }
00137   else {
00138     BlkInfo<LocalOrdinal,Scalar>& blk = iter->second;
00139     numPtRows = blk.numPtRows;
00140     numPtCols = blk.numPtCols;
00141     blockEntry = blk.blkEntry;
00142   }
00143 }
00144 
00145 template<typename LocalOrdinal, typename GlobalOrdinal, typename Scalar, class Node>
00146 Teuchos::RCP<const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node> >
00147 createOverlapMap(VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrdata,
00148                 const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& rowmap)
00149 {
00150   Teuchos::Array<GlobalOrdinal> importBlocks;
00151   Teuchos::Array<int> src_procs;
00152 
00153   typename std::map<GlobalOrdinal,LocalOrdinal>::iterator
00154     iter = vbrdata.row_map.begin(),
00155     iter_end = vbrdata.row_map.end();
00156   for(; iter != iter_end; ++iter) {
00157     importBlocks.push_back(iter->first);
00158   }
00159 
00160   Teuchos::Array<GlobalOrdinal> importFirstPointInBlocks(importBlocks.size());
00161   Teuchos::Array<LocalOrdinal>  importBlockSizes(importBlocks.size());
00162   rowmap.getRemoteBlockInfo(importBlocks(), importFirstPointInBlocks(), importBlockSizes());
00163 
00164   return Teuchos::rcp(
00165      new Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>(
00166            Teuchos::OrdinalTraits<global_size_t>::invalid(),
00167            importBlocks(),
00168            importFirstPointInBlocks(),
00169            importBlockSizes(),
00170            rowmap.getPointMap()->getIndexBase(),
00171            rowmap.getPointMap()->getComm(),
00172            rowmap.getPointMap()->getNode()
00173          )
00174   );
00175 }
00176 
00177 template<class LocalOrdinal, class GlobalOrdinal, class Scalar, class Node>
00178 class VbrDataDist : public Tpetra::SrcDistObject, 
00179         public Tpetra::Packable<char, LocalOrdinal> {
00180 private:
00181   Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > pointMap_;
00182   VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrData_;
00183 
00184 public:
00185   VbrDataDist (VbrData<LocalOrdinal,GlobalOrdinal,Scalar>& vbrData,
00186          const Tpetra::BlockMap<LocalOrdinal,GlobalOrdinal,Node>& importMap)
00187     : pointMap_ (convertBlockMapToPointMap (importMap)),
00188       vbrData_ (vbrData)
00189   {}
00190 
00191   ~VbrDataDist() {}
00192 
00193   Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > 
00194   getMap () const { return pointMap_; }
00195 
00196   virtual void 
00197   pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
00198   Teuchos::Array<char>& exports,
00199   const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00200   size_t& constantNumPackets,
00201   Distributor& distor) const
00202   {
00203     (void) distor; // forestall "unused argument" compiler warning
00204     using Teuchos::Array;
00205     using Teuchos::ArrayView;
00206     using Teuchos::as;
00207     using Teuchos::av_reinterpret_cast;
00208     typedef LocalOrdinal LO;
00209     typedef GlobalOrdinal GO;
00210     typedef typename ArrayView<const LO>::size_type Tsize_t;
00211     typedef typename VbrData<LO,GO,Scalar>::RowGlobalCols RowGlobalCols;
00212 
00213     //We will pack each row's data into the exports buffer as follows:
00214     //[num-block-cols,numPtRows,{list-of-blk-cols},{list-of-ptColsPerBlkCol},{all vals}]
00215     //so the length of the char exports buffer for a row is:
00216     //sizeof(LO)*(2+2*(num-block-cols)) + sizeof(Scalar)*numPtRows*sum(numPtCols_i)
00217 
00218     //For each row corresponding to exportLIDs, accumulate the size that it will
00219     //occupy in the exports buffer:
00220     size_t total_exports_size = 0;
00221     for (Tsize_t i = 0; i < exportLIDs.size (); ++i) {
00222       const RowGlobalCols& rgc = (this->vbrData_.data)[exportLIDs[i]];
00223       typename RowGlobalCols::const_iterator rgc_it = rgc.begin();
00224       typename RowGlobalCols::const_iterator rgc_end = rgc.end();
00225       size_t numScalars = 0;
00226       for ( ; rgc_it != rgc_end; ++rgc_it) {
00227         const BlkInfo<LO,Scalar>& blk = rgc_it->second;
00228         numScalars += as<size_t> (blk.numPtRows) * as<size_t> (blk.numPtCols);
00229       }
00230       const size_t numBlkCols = as<size_t> (rgc.size ());
00231       const size_t size_for_this_row = 
00232   sizeof (GO) * (2 + 2*numBlkCols)
00233   + sizeof (Scalar) * numScalars;
00234       numPacketsPerLID[i] = size_for_this_row;
00235       total_exports_size += size_for_this_row;
00236     }
00237 
00238     exports.resize (total_exports_size);
00239 
00240     ArrayView<char> avIndsC, avValsC;
00241     ArrayView<Scalar> avVals;
00242 
00243     Array<GO> blkCols;
00244     Array<LO> ptColsPerBlkCol;
00245     Array<Scalar> blkEntries;
00246 
00247     size_t offset = 0;
00248     for (Tsize_t i = 0; i < exportLIDs.size (); ++i) {
00249       blkCols.clear();
00250       ptColsPerBlkCol.clear();
00251       blkEntries.clear();
00252 
00253       const RowGlobalCols& rgc = (this->vbrData_.data)[exportLIDs[i]];
00254       typename RowGlobalCols::const_iterator rgc_it = rgc.begin ();
00255       typename RowGlobalCols::const_iterator rgc_end = rgc.end ();
00256       LO numPtRows = 0;
00257       for ( ; rgc_it != rgc_end; ++rgc_it) {
00258         blkCols.push_back (rgc_it->first);
00259         const BlkInfo<LO,Scalar>& blk = rgc_it->second;
00260         numPtRows = blk.numPtRows; // should be the same for all cols
00261         ptColsPerBlkCol.push_back (blk.numPtCols);
00262         for (LO j = 0; j < blk.numPtRows * blk.numPtCols; ++j) {
00263           blkEntries.push_back (blk.blkEntry[j]);
00264         }
00265       }
00266 
00267       LO numBlkCols = blkCols.size();
00268       LO numScalars = blkEntries.size();
00269 
00270       LO num_chars_for_ordinals = (2*numBlkCols + 2) * sizeof (GO);
00271       //get export views
00272       avIndsC = exports(offset, num_chars_for_ordinals);
00273       avValsC = exports(offset+ num_chars_for_ordinals, numScalars*sizeof(Scalar));
00274       ArrayView<GO> avInds = av_reinterpret_cast<GO>(avIndsC);
00275       typename ArrayView<GO>::iterator ind_it = avInds.begin ();
00276 
00277       *ind_it++ = numBlkCols;
00278       *ind_it++ = numPtRows;
00279       for (Tsize_t j = 0; j < blkCols.size (); ++j) {
00280         *ind_it++ = blkCols[j];
00281       }
00282       for (Tsize_t j = 0; j < ptColsPerBlkCol.size (); ++j) {
00283         *ind_it++ = ptColsPerBlkCol[j];
00284       }
00285 
00286       avVals = av_reinterpret_cast<Scalar> (avValsC);
00287       std::copy (blkEntries.begin (), blkEntries.end (), avVals.begin ());
00288 
00289       const size_t size_for_this_row = 
00290   sizeof (GO) * (2 + 2*numBlkCols) + 
00291   sizeof (Scalar) * numScalars;
00292       offset += size_for_this_row;
00293     }
00294     constantNumPackets = 0;
00295   }
00296 
00297 };
00298 
00299 }//namespace VbrUtils
00300 }//namespace Tpetra
00301 
00302 #endif
00303 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines