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