Tpetra Matrix/Vector Services Version of the Day
TpetraExt_MMHelpers_def.hpp
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_MMHELPERS_DEF_HPP
00043 #define TPETRA_MMHELPERS_DEF_HPP
00044 
00045 #include "Teuchos_VerboseObject.hpp"
00046 #include "Tpetra_ConfigDefs.hpp"
00047 #include "Tpetra_CrsMatrix.hpp"
00048 #include "TpetraExt_MMHelpers_decl.hpp"
00049 
00054 namespace Tpetra {
00055 
00056 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00057 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::CrsMatrixStruct()
00058 {
00059 }
00060 
00061 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00062 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::~CrsMatrixStruct()
00063 {
00064   deleteContents();
00065 }
00066 
00067 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00068 void CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::deleteContents()
00069 {
00070   numRows = 0;
00071   numEntriesPerRow.clear();
00072   indices.clear();
00073   values.clear();
00074   remote.clear();
00075   numRemote = 0;
00076   importMatrix.reset();
00077 }
00078 
00079 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00080 int dumpCrsMatrixStruct(const CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>& M)
00081 {
00082   std::cout << "proc " << M.rowMap->Comm().MyPID()<<std::endl;
00083   std::cout << "numRows: " << M.numRows<<std::endl;
00084   for(LocalOrdinal i=0; i<M.numRows; ++i) {
00085     for(LocalOrdinal j=0; j<M.numEntriesPerRow[i]; ++j) {
00086       if (M.remote[i]) {
00087         std::cout << "  *"<<M.rowMap->GID(i)<<"   "
00088              <<M.importColMap->GID(M.indices[i][j])<<"   "<<M.values[i][j]<<std::endl;
00089       }
00090       else {
00091         std::cout << "   "<<M.rowMap->GID(i)<<"   "
00092              <<M.colMap->GID(M.indices[i][j])<<"   "<<M.values[i][j]<<std::endl;
00093       }
00094     }
00095   }
00096   return(0);
00097 }
00098 
00099 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00100 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::CrsWrapper_CrsMatrix(CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>& crsmatrix)
00101  : crsmat_(crsmatrix)
00102 {
00103 }
00104 
00105 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00106 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::~CrsWrapper_CrsMatrix()
00107 {
00108 }
00109 
00110 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00111 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
00112 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::getRowMap() const
00113 {
00114   return crsmat_.getRowMap();
00115 }
00116 
00117 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00118 bool CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::isFillComplete()
00119 {
00120   return crsmat_.isFillComplete();
00121 }
00122 
00123 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00124 void
00125 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::insertGlobalValues(GlobalOrdinal globalRow, const Teuchos::ArrayView<const GlobalOrdinal> &indices, const Teuchos::ArrayView<const Scalar> &values)
00126 {
00127   crsmat_.insertGlobalValues(globalRow, indices, values);
00128 }
00129 
00130 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00131 void 
00132 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::sumIntoGlobalValues(GlobalOrdinal globalRow, const Teuchos::ArrayView<const GlobalOrdinal> &indices, const Teuchos::ArrayView<const Scalar> &values)
00133 {
00134   crsmat_.sumIntoGlobalValues(globalRow, indices, values);
00135 }
00136 
00137 
00138 //------------------------------------
00139 
00140 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00141 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CrsWrapper_GraphBuilder(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map)
00142  : graph_(),
00143    rowmap_(map),
00144    max_row_length_(0)
00145 {
00146   Teuchos::ArrayView<const GlobalOrdinal> rows= map->getNodeElementList();
00147 
00148   for(int i=0; i<rows.size(); ++i) {
00149     graph_[rows[i]] = new std::set<GlobalOrdinal>;
00150   }
00151 }
00152 
00153 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00154 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~CrsWrapper_GraphBuilder()
00155 {
00156   typename std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>::iterator
00157     iter = graph_.begin(), iter_end = graph_.end();
00158   for(; iter!=iter_end; ++iter) {
00159     delete iter->second;
00160   }
00161 
00162   graph_.clear();
00163 }
00164 
00165 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00166 bool CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isFillComplete()
00167 {
00168   return false;
00169 }
00170 
00171 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00172 void
00173 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::insertGlobalValues(GlobalOrdinal globalRow, const Teuchos::ArrayView<const GlobalOrdinal> &indices, const Teuchos::ArrayView<const Scalar> &values)
00174 {
00175   /*std::cout << "inserting" << std::endl;
00176   std::cout << "  indices: " << indices << std::endl;
00177   std::cout << "  values: " << values << std::endl;*/
00178   typename std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>::iterator
00179     iter = graph_.find(globalRow);
00180 
00181   TEUCHOS_TEST_FOR_EXCEPTION(iter == graph_.end(), std::runtime_error,
00182   Teuchos::typeName(*this) << "::insertGlobalValues could not find row " << globalRow << " in the graph. Super bummer man. Hope you figure it out.\n");
00183 
00184   std::set<GlobalOrdinal>& cols = *(iter->second);
00185 
00186   for(int i=0; i<indices.size(); ++i) {
00187     cols.insert(indices[i]);
00188   }
00189 
00190   global_size_t row_length = cols.size();
00191   if (row_length > max_row_length_) max_row_length_ = row_length;
00192 
00193   
00194 }
00195 
00196 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00197 void
00198 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sumIntoGlobalValues(GlobalOrdinal globalRow, const Teuchos::ArrayView<const GlobalOrdinal> &indices, const Teuchos::ArrayView<const Scalar> &values)
00199 {
00200   insertGlobalValues(globalRow, indices, values);
00201 }
00202 
00203 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00204 std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>&
00205 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::get_graph()
00206 {
00207   return graph_;
00208 }
00209 
00210 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class SpMatOps>
00211 void insert_matrix_locations(CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node> & graphbuilder,
00212   CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>& C)
00213 {
00214   global_size_t max_row_length = graphbuilder.get_max_row_length();
00215   if (max_row_length < 1) return;
00216 
00217   Array<GlobalOrdinal> indices(max_row_length);
00218   Array<Scalar> zeros(max_row_length, 0.0);
00219 
00220   typedef std::map<GlobalOrdinal,std::set<GlobalOrdinal>*> Graph;
00221   typedef typename Graph::iterator    GraphIter;
00222   typedef std::set<GlobalOrdinal>     Set;
00223   typedef typename Set::iterator      SetIter;
00224   Graph& graph = graphbuilder.get_graph();
00225 
00226   const GraphIter iter_end = graph.end();
00227   for(GraphIter iter=graph.begin(); iter!=iter_end; ++iter) {
00228     const GlobalOrdinal row = iter->first;
00229     const std::set<GlobalOrdinal> &cols = *(iter->second);
00230     // "copy" entries out of set into contiguous array storage
00231     const size_t num_entries = std::copy(cols.begin(), cols.end(), indices.begin()) - indices.begin();
00232     // insert zeros into the result matrix at the appropriate locations
00233     C.insertGlobalValues(row, indices(0,num_entries), zeros(0,num_entries));
00234   }
00235 }
00236 
00237 } //namespace Tpetra
00238 
00239 //
00240 // Explicit instantiation macro
00241 //
00242 // Must be expanded from within the Tpetra namespace!
00243 //
00244 
00245 #define TPETRA_CRSMATRIXSTRUCT_INSTANT(SCALAR,LO,GO,NODE) \
00246   \
00247   template class CrsMatrixStruct< SCALAR , LO , GO , NODE >;
00248 
00249 #define TPETRA_CRSWRAPPER_INSTANT(SCALAR,LO,GO,NODE) \
00250   \
00251   template class CrsWrapper< SCALAR , LO , GO , NODE >;
00252 
00253 #define TPETRA_CRSWRAPPER_CRSMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
00254   \
00255   template class CrsWrapper_CrsMatrix< SCALAR , LO , GO , NODE >;
00256 
00257 #define TPETRA_CRSWRAPPER_GRAPHBUILDER_INSTANT(SCALAR,LO,GO,NODE) \
00258   \
00259   template class CrsWrapper_GraphBuilder< SCALAR , LO , GO , NODE >;
00260 
00261 #endif // TPETRA_MMHELPERS_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines