Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Zoltan2_ModelHelpers.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 
50 #include <Zoltan2_MeshAdapter.hpp>
51 
52 #ifndef _ZOLTAN2_MODELHELPERS_HPP_
53 #define _ZOLTAN2_MODELHELPERS_HPP_
54 
55 namespace Zoltan2 {
56 
57 //GFD this declaration is really messy is there a better way? I couldn't typedef outside since
58 // there is no user type until the function.
59 template <typename User>
60 RCP<Tpetra::CrsMatrix<int,
61  typename MeshAdapter<User>::lno_t,
62  typename MeshAdapter<User>::gno_t,
63  typename MeshAdapter<User>::node_t> >
64 get2ndAdjsMatFromAdjs(const Teuchos::RCP<const MeshAdapter<User> > &ia,
65  const RCP<const Comm<int> > comm,
66  Zoltan2::MeshEntityType sourcetarget,
67  Zoltan2::MeshEntityType through) {
68  typedef typename MeshAdapter<User>::gno_t gno_t;
69  typedef typename MeshAdapter<User>::lno_t lno_t;
70  typedef typename MeshAdapter<User>::node_t node_t;
71  typedef typename MeshAdapter<User>::offset_t offset_t;
72 
73  typedef int nonzero_t; // adjacency matrix doesn't need scalar_t
74  typedef Tpetra::CrsMatrix<nonzero_t,lno_t,gno_t,node_t> sparse_matrix_type;
75  typedef Tpetra::Map<lno_t, gno_t, node_t> map_type;
76  typedef Tpetra::global_size_t GST;
77  const GST dummy = Teuchos::OrdinalTraits<GST>::invalid ();
78 
79 /* Find the adjacency for a nodal based decomposition */
80  if (ia->availAdjs(sourcetarget, through)) {
81  using Teuchos::Array;
82  using Teuchos::as;
83  using Teuchos::RCP;
84  using Teuchos::rcp;
85 
86  // Get node-element connectivity
87 
88  const offset_t *offsets=NULL;
89  const gno_t *adjacencyIds=NULL;
90  ia->getAdjsView(sourcetarget, through, offsets, adjacencyIds);
91 
92  const gno_t *Ids=NULL;
93  ia->getIDsViewOf(sourcetarget, Ids);
94 
95  const gno_t *throughIds=NULL;
96  ia->getIDsViewOf(through, throughIds);
97 
98  size_t LocalNumIDs = ia->getLocalNumOf(sourcetarget);
99 
100  /***********************************************************************/
101  /************************* BUILD MAPS FOR ADJS *************************/
102  /***********************************************************************/
103 
104  RCP<const map_type> sourcetargetMapG;
105  RCP<const map_type> throughMapG;
106 
107  // count owned nodes
108  size_t LocalNumOfThrough = ia->getLocalNumOf(through);
109 
110  // Build a list of the global sourcetarget ids...
111  gno_t min[2];
112  size_t maxcols = 0;
113  min[0] = std::numeric_limits<gno_t>::max();
114  for (size_t i = 0; i < LocalNumIDs; ++i) {
115  if (Ids[i] < min[0]) {
116  min[0] = Ids[i];
117  }
118  size_t ncols = offsets[i+1] - offsets[i];
119  if (ncols > maxcols) maxcols = ncols;
120  }
121 
122  // min(throughIds[i])
123  min[1] = std::numeric_limits<gno_t>::max();
124  for (size_t i = 0; i < LocalNumOfThrough; ++i) {
125  if (throughIds[i] < min[1]) {
126  min[1] = throughIds[i];
127  }
128  }
129 
130  gno_t gmin[2];
131  Teuchos::reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MIN, 2, min, gmin);
132 
133  //Generate Map for sourcetarget.
134  ArrayView<const gno_t> sourceTargetGIDs(Ids, LocalNumIDs);
135  sourcetargetMapG = rcp(new map_type(dummy,
136  sourceTargetGIDs, gmin[0], comm));
137 
138  //Create a new map with IDs uniquely assigned to ranks (oneToOneSTMap)
139  /*RCP<const map_type> oneToOneSTMap =
140  Tpetra::createOneToOne<lno_t, gno_t, node_t>(sourcetargetMapG);*/
141 
142  //Generate Map for through.
143 // TODO
144 // TODO Could check for max through id as well, and if all through ids are
145 // TODO in gmin to gmax, then default constructors works below.
146 // TODO Otherwise, may need a constructor that is not one-to-one containing
147 // TODO all through entities on processor, followed by call to createOneToOne
148 // TODO
149 
150  ArrayView<const gno_t> throughGIDs(throughIds, LocalNumOfThrough);
151  throughMapG = rcp (new map_type(dummy,
152  throughGIDs, gmin[1], comm));
153 
154  //Create a new map with IDs uniquely assigned to ranks (oneToOneTMap)
155  RCP<const map_type> oneToOneTMap =
156  Tpetra::createOneToOne<lno_t, gno_t, node_t>(throughMapG);
157 
158  /***********************************************************************/
159  /************************* BUILD GRAPH FOR ADJS ************************/
160  /***********************************************************************/
161 
162  RCP<sparse_matrix_type> adjsMatrix;
163 
164  // Construct Tpetra::CrsGraph objects.
165  adjsMatrix = rcp (new sparse_matrix_type (sourcetargetMapG,//oneToOneSTMap,
166  0));
167 
168  Array<nonzero_t> justOneA(maxcols, 1);
169  ArrayView<const gno_t> adjacencyIdsAV(adjacencyIds, offsets[LocalNumIDs]);
170 
171  for (size_t localElement=0; localElement<LocalNumIDs; ++localElement){
172  // Insert all columns for global row Ids[localElement]
173  size_t ncols = offsets[localElement+1] - offsets[localElement];
174  adjsMatrix->insertGlobalValues(Ids[localElement],
175  adjacencyIdsAV(offsets[localElement], ncols),
176  justOneA(0, ncols));
177  }// *** source loop ***
178 
179  //Fill-complete adjs Graph
180  adjsMatrix->fillComplete (oneToOneTMap, //throughMapG,
181  adjsMatrix->getRowMap());
182 
183  // Form 2ndAdjs
184  RCP<sparse_matrix_type> secondAdjs =
185  rcp (new sparse_matrix_type(adjsMatrix->getRowMap(),0));
186  Tpetra::MatrixMatrix::Multiply(*adjsMatrix,false,*adjsMatrix,
187  true,*secondAdjs);
188  return secondAdjs;
189  }
190  return RCP<sparse_matrix_type>();
191 }
192 
193 template <typename User>
195  const Teuchos::RCP<const MeshAdapter<User> > &ia,
196  const RCP<const Comm<int> > comm,
197  Zoltan2::MeshEntityType sourcetarget,
198  Zoltan2::MeshEntityType through,
199  Teuchos::ArrayRCP<typename MeshAdapter<User>::offset_t> &offsets,
200  Teuchos::ArrayRCP<typename MeshAdapter<User>::gno_t> &adjacencyIds)
201 {
202  typedef typename MeshAdapter<User>::gno_t gno_t;
203  typedef typename MeshAdapter<User>::lno_t lno_t;
204  typedef typename MeshAdapter<User>::offset_t offset_t;
205  typedef typename MeshAdapter<User>::node_t node_t;
206 
207  typedef int nonzero_t; // adjacency matrix doesn't need scalar_t
208  typedef Tpetra::CrsMatrix<nonzero_t,lno_t,gno_t,node_t> sparse_matrix_type;
209 
210  RCP<sparse_matrix_type> secondAdjs = get2ndAdjsMatFromAdjs(ia,comm,sourcetarget,through);
211 
212  /* Allocate memory necessary for the adjacency */
213  size_t LocalNumIDs = ia->getLocalNumOf(sourcetarget);
214  lno_t *start = new lno_t [LocalNumIDs+1];
215 
216  if (secondAdjs!=RCP<sparse_matrix_type>()) {
217  Array<gno_t> Indices;
218  Array<nonzero_t> Values;
219 
220  size_t nadj = 0;
221 
222  gno_t const *Ids=NULL;
223  ia->getIDsViewOf(sourcetarget, Ids);
224 
225  std::vector<gno_t> adj;
226 
227  for (size_t localElement=0; localElement<LocalNumIDs; ++localElement){
228  start[localElement] = nadj;
229  const gno_t globalRow = Ids[localElement];
230  size_t NumEntries = secondAdjs->getNumEntriesInGlobalRow (globalRow);
231  Indices.resize (NumEntries);
232  Values.resize (NumEntries);
233  secondAdjs->getGlobalRowCopy (globalRow,Indices(),Values(),NumEntries);
234 
235  for (size_t j = 0; j < NumEntries; ++j) {
236  if(globalRow != Indices[j]) {
237  adj.push_back(Indices[j]);
238  nadj++;;
239  }
240  }
241  }
242 
243  Ids = NULL;
244  start[LocalNumIDs] = nadj;
245 
246  gno_t *adj_ = new gno_t [nadj];
247 
248  for (size_t i=0; i < nadj; i++) {
249  adj_[i] = adj[i];
250  }
251 
252  offsets = arcp<offset_t>(start, 0, LocalNumIDs+1, true);
253  adjacencyIds = arcp<gno_t>(adj_, 0, nadj, true);
254  }
255  else {
256  // No adjacencies could be computed; return no edges and valid offsets array
257  for (size_t i = 0; i <= LocalNumIDs; i++)
258  start[i] = 0;
259 
260  offsets = arcp<offset_t>(start, 0, LocalNumIDs+1, true);
261  adjacencyIds = Teuchos::null;
262  }
263 
264  //return nadj;
265 }
266 
267 }
268 
269 #endif
RCP< Tpetra::CrsMatrix< int, typename MeshAdapter< User >::lno_t, typename MeshAdapter< User >::gno_t, typename MeshAdapter< User >::node_t > > get2ndAdjsMatFromAdjs(const Teuchos::RCP< const MeshAdapter< User > > &ia, const RCP< const Comm< int > > comm, Zoltan2::MeshEntityType sourcetarget, Zoltan2::MeshEntityType through)
void get2ndAdjsViewFromAdjs(const Teuchos::RCP< const MeshAdapter< User > > &ia, const RCP< const Comm< int > > comm, Zoltan2::MeshEntityType sourcetarget, Zoltan2::MeshEntityType through, Teuchos::ArrayRCP< typename MeshAdapter< User >::offset_t > &offsets, Teuchos::ArrayRCP< typename MeshAdapter< User >::gno_t > &adjacencyIds)
InputTraits< User >::gno_t gno_t
Defines the MeshAdapter interface.
MeshAdapter defines the interface for mesh input.
InputTraits< User >::lno_t lno_t
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
Tpetra::global_size_t global_size_t
InputTraits< User >::offset_t offset_t