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