Zoltan2
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 
72  typedef int nonzero_t; // adjacency matrix doesn't need scalar_t
73  typedef Tpetra::CrsMatrix<nonzero_t,lno_t,gno_t,node_t> sparse_matrix_type;
74  typedef Tpetra::Map<lno_t, gno_t, node_t> map_type;
75  typedef Tpetra::global_size_t GST;
76  const GST dummy = Teuchos::OrdinalTraits<GST>::invalid ();
77 
78 /* Find the adjacency for a nodal based decomposition */
79  if (ia->availAdjs(sourcetarget, through)) {
80  using Tpetra::DefaultPlatform;
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 lno_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>
194 void get2ndAdjsViewFromAdjs(const Teuchos::RCP<const MeshAdapter<User> > &ia,
195  const RCP<const Comm<int> > comm,
196  Zoltan2::MeshEntityType sourcetarget,
197  Zoltan2::MeshEntityType through,
198  const typename MeshAdapter<User>::lno_t *&offsets,
199  const typename MeshAdapter<User>::gno_t *&adjacencyIds)
200 {
201  typedef typename MeshAdapter<User>::gno_t gno_t;
202  typedef typename MeshAdapter<User>::lno_t lno_t;
203  typedef typename MeshAdapter<User>::node_t node_t;
204 
205  typedef int nonzero_t; // adjacency matrix doesn't need scalar_t
206  typedef Tpetra::CrsMatrix<nonzero_t,lno_t,gno_t,node_t> sparse_matrix_type;
207 
208  RCP<sparse_matrix_type> secondAdjs = get2ndAdjsMatFromAdjs(ia,comm,sourcetarget,through);
209 
210  if (secondAdjs!=RCP<sparse_matrix_type>()) {
211  Array<gno_t> Indices;
212  Array<nonzero_t> Values;
213 
214  size_t nadj = 0;
215 
216  gno_t const *Ids=NULL;
217  ia->getIDsViewOf(sourcetarget, Ids);
218 
219  /* Allocate memory necessary for the adjacency */
220  size_t LocalNumIDs = ia->getLocalNumOf(sourcetarget);
221  lno_t *start = new lno_t [LocalNumIDs+1];
222  std::vector<gno_t> adj;
223 
224  for (size_t localElement=0; localElement<LocalNumIDs; ++localElement){
225  start[localElement] = nadj;
226  const gno_t globalRow = Ids[localElement];
227  size_t NumEntries = secondAdjs->getNumEntriesInGlobalRow (globalRow);
228  Indices.resize (NumEntries);
229  Values.resize (NumEntries);
230  secondAdjs->getGlobalRowCopy (globalRow,Indices(),Values(),NumEntries);
231 
232  for (size_t j = 0; j < NumEntries; ++j) {
233  if(globalRow != Indices[j]) {
234  adj.push_back(Indices[j]);
235  nadj++;;
236  }
237  }
238  }
239 
240  Ids = NULL;
241  start[LocalNumIDs] = nadj;
242 
243  gno_t *adj_ = new gno_t [nadj];
244 
245  for (size_t i=0; i < nadj; i++) {
246  adj_[i] = adj[i];
247  }
248 
249  offsets = start;
250  adjacencyIds = adj_;
251  }
252 
253  //return nadj;
254 }
255 
256 }
257 
258 #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)
InputTraits< User >::gno_t gno_t
size_t global_size_t
Defines the MeshAdapter interface.
MeshAdapter defines the interface for mesh input.
void get2ndAdjsViewFromAdjs(const Teuchos::RCP< const MeshAdapter< User > > &ia, const RCP< const Comm< int > > comm, Zoltan2::MeshEntityType sourcetarget, Zoltan2::MeshEntityType through, const typename MeshAdapter< User >::lno_t *&offsets, const typename MeshAdapter< User >::gno_t *&adjacencyIds)
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.