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