Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TpetraExt_MMHelpers_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef TPETRA_MMHELPERS_DEF_HPP
43 #define TPETRA_MMHELPERS_DEF_HPP
44 
46 #include "Teuchos_VerboseObject.hpp"
47 
52 namespace Tpetra {
53 
54 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CrsMatrixStruct()
56 {
57 }
58 
59 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
60 CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~CrsMatrixStruct()
61 {
62  deleteContents();
63 }
64 
65 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66 void CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
67 deleteContents ()
68 {
69  importMatrix.reset();
70  origMatrix = Teuchos::null;
71 }
72 
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 BlockCrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BlockCrsMatrixStruct(const LocalOrdinal blocksize_)
75  : blocksize(blocksize_)
76 {
77 }
78 
79 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80 BlockCrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~BlockCrsMatrixStruct()
81 {
82  deleteContents();
83 }
84 
85 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86 void BlockCrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
87 deleteContents ()
88 {
89  importMatrix.reset();
90  origMatrix = Teuchos::null;
91 }
92 
93 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
94 int dumpCrsMatrixStruct (const CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& M)
95 {
96  std::cout << "proc " << M.rowMap->Comm().MyPID()<<std::endl;
97  std::cout << "numRows: " << M.numRows<<std::endl;
98  for(LocalOrdinal i=0; i<M.numRows; ++i) {
99  for(LocalOrdinal j=0; j<M.numEntriesPerRow[i]; ++j) {
100  std::cout << " "<<M.rowMap->GID(i)<<" "
101  <<M.colMap->GID(M.indices[i][j])<<" "<<M.values[i][j]<<std::endl;
102  }
103  }
104 
105  return 0;
106 }
107 
108 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
109 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
110 CrsWrapper_CrsMatrix (CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsmatrix)
111  : crsmat_ (crsmatrix)
112 {
113 }
114 
115 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
116 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~CrsWrapper_CrsMatrix()
117 {
118 }
119 
120 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
122 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getRowMap() const
123 {
124  return crsmat_.getRowMap();
125 }
126 
127 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
128 bool CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
129 isFillComplete ()
130 {
131  return crsmat_.isFillComplete ();
132 }
133 
134 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
135 void
136 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
137 insertGlobalValues (GlobalOrdinal globalRow,
138  const Teuchos::ArrayView<const GlobalOrdinal> &indices,
139  const Teuchos::ArrayView<const Scalar> &values)
140 {
141  crsmat_.insertGlobalValues (globalRow, indices, values);
142 }
143 
144 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
145 void
146 CrsWrapper_CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
147 sumIntoGlobalValues (GlobalOrdinal globalRow,
148  const Teuchos::ArrayView<const GlobalOrdinal> &indices,
149  const Teuchos::ArrayView<const Scalar> &values)
150 {
151  crsmat_.sumIntoGlobalValues (globalRow, indices, values);
152 }
153 
154 
155 
156 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
157 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
158 CrsWrapper_GraphBuilder (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map)
159  : graph_(),
160  rowmap_(map),
161  max_row_length_(0)
162 {
163  Teuchos::ArrayView<const GlobalOrdinal> rows = map->getLocalElementList ();
164  const LocalOrdinal numRows = static_cast<LocalOrdinal> (rows.size ());
165  for (LocalOrdinal i = 0; i < numRows; ++i) {
166  graph_[rows[i]] = new std::set<GlobalOrdinal>;
167  }
168 }
169 
170 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
171 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
172 ~CrsWrapper_GraphBuilder ()
173 {
174  typename std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>::iterator
175  iter = graph_.begin(), iter_end = graph_.end();
176  for (; iter != iter_end; ++iter) {
177  delete iter->second;
178  }
179  graph_.clear ();
180 }
181 
182 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
183 bool CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isFillComplete()
184 {
185  return false;
186 }
187 
188 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
189 void
190 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
191 insertGlobalValues (GlobalOrdinal globalRow,
192  const Teuchos::ArrayView<const GlobalOrdinal> &indices,
193  const Teuchos::ArrayView<const Scalar> &/* values */)
194 {
195  typename std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>::iterator
196  iter = graph_.find (globalRow);
197 
198  TEUCHOS_TEST_FOR_EXCEPTION(
199  iter == graph_.end(), std::runtime_error,
200  "Tpetra::CrsWrapper_GraphBuilder::insertGlobalValues could not find row "
201  << globalRow << " in the graph. Super bummer man. Hope you figure it out.");
202 
203  std::set<GlobalOrdinal>& cols = * (iter->second);
204 
205  for (typename Teuchos::ArrayView<const GlobalOrdinal>::size_type i = 0;
206  i < indices.size (); ++i) {
207  cols.insert (indices[i]);
208  }
209 
210  const global_size_t row_length = static_cast<global_size_t> (cols.size ());
211  if (row_length > max_row_length_) {
212  max_row_length_ = row_length;
213  }
214 }
215 
216 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
217 void
218 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
219 sumIntoGlobalValues (GlobalOrdinal globalRow,
220  const Teuchos::ArrayView<const GlobalOrdinal> &indices,
221  const Teuchos::ArrayView<const Scalar> &values)
222 {
223  insertGlobalValues (globalRow, indices, values);
224 }
225 
226 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
227 std::map<GlobalOrdinal,std::set<GlobalOrdinal>*>&
228 CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::get_graph ()
229 {
230  return graph_;
231 }
232 
233 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
234 void
235 insert_matrix_locations (CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>& graphbuilder,
236  CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& C)
237 {
238  global_size_t max_row_length = graphbuilder.get_max_row_length();
239  if (max_row_length < 1) return;
240 
241  Teuchos::Array<GlobalOrdinal> indices(max_row_length);
242  Teuchos::Array<Scalar> zeros(max_row_length, Teuchos::ScalarTraits<Scalar>::zero());
243 
244  typedef std::map<GlobalOrdinal,std::set<GlobalOrdinal>*> Graph;
245  typedef typename Graph::iterator GraphIter;
246  Graph& graph = graphbuilder.get_graph ();
247 
248  const GraphIter iter_end = graph.end ();
249  for (GraphIter iter = graph.begin (); iter != iter_end; ++iter) {
250  const GlobalOrdinal row = iter->first;
251  const std::set<GlobalOrdinal>& cols = * (iter->second);
252  // "copy" entries out of set into contiguous array storage
253  const size_t num_entries = std::copy (cols.begin (), cols.end (), indices.begin ()) - indices.begin ();
254  // insert zeros into the result matrix at the appropriate locations
255  C.insertGlobalValues (row, indices (0, num_entries), zeros (0, num_entries));
256  }
257 }
258 
259 } // namespace Tpetra
260 
261 //
262 // Explicit instantiation macro
263 //
264 // Must be expanded from within the Tpetra namespace!
265 //
266 
267 #define TPETRA_CRSMATRIXSTRUCT_INSTANT(SCALAR,LO,GO,NODE) \
268  \
269  template class CrsMatrixStruct< SCALAR , LO , GO , NODE >;
270 
271 #define TPETRA_BLOCKCRSMATRIXSTRUCT_INSTANT(SCALAR,LO,GO,NODE) \
272  \
273  template class BlockCrsMatrixStruct< SCALAR , LO , GO , NODE >;
274 
275 #define TPETRA_CRSWRAPPER_INSTANT(SCALAR,LO,GO,NODE) \
276  \
277  template class CrsWrapper< SCALAR , LO , GO , NODE >;
278 
279 #define TPETRA_CRSWRAPPER_CRSMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
280  \
281  template class CrsWrapper_CrsMatrix< SCALAR , LO , GO , NODE >;
282 
283 #define TPETRA_CRSWRAPPER_GRAPHBUILDER_INSTANT(SCALAR,LO,GO,NODE) \
284  \
285  template class CrsWrapper_GraphBuilder< SCALAR , LO , GO , NODE >;
286 
287 #endif // TPETRA_MMHELPERS_DEF_HPP
size_t global_size_t
Global size_t object.
Declaration of Tpetra::MMMultiMultiply and nonmember constructors.