Zoltan2
Loading...
Searching...
No Matches
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
51
52#ifndef _ZOLTAN2_MODELHELPERS_HPP_
53#define _ZOLTAN2_MODELHELPERS_HPP_
54
55namespace 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.
59template <typename User>
60RCP<Tpetra::CrsMatrix<int,
64get2ndAdjsMatFromAdjs(const Teuchos::RCP<const MeshAdapter<User> > &ia,
65 const RCP<const Comm<int> > comm,
66 Zoltan2::MeshEntityType sourcetarget,
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 Array<size_t> rowlens(LocalNumIDs);
166 for (size_t localElement = 0; localElement < LocalNumIDs; localElement++)
167 rowlens[localElement] = offsets[localElement+1] - offsets[localElement];
168 adjsMatrix = rcp (new sparse_matrix_type (sourcetargetMapG,//oneToOneSTMap,
169 rowlens()));
170
171 Array<nonzero_t> justOneA(maxcols, 1);
172 ArrayView<const gno_t> adjacencyIdsAV(adjacencyIds, offsets[LocalNumIDs]);
173
174 for (size_t localElement=0; localElement<LocalNumIDs; ++localElement){
175 // Insert all columns for global row Ids[localElement]
176 size_t ncols = offsets[localElement+1] - offsets[localElement];
177 adjsMatrix->insertGlobalValues(Ids[localElement],
178 adjacencyIdsAV(offsets[localElement], ncols),
179 justOneA(0, ncols));
180 }// *** source loop ***
181
182 //Fill-complete adjs Graph
183 adjsMatrix->fillComplete (oneToOneTMap, //throughMapG,
184 adjsMatrix->getRowMap());
185
186 // Form 2ndAdjs
187 RCP<sparse_matrix_type> secondAdjs =
188 rcp (new sparse_matrix_type(adjsMatrix->getRowMap(),0));
189 Tpetra::MatrixMatrix::Multiply(*adjsMatrix,false,*adjsMatrix,
190 true,*secondAdjs);
191 return secondAdjs;
192 }
193 return RCP<sparse_matrix_type>();
194}
195
196template <typename User>
198 const Teuchos::RCP<const MeshAdapter<User> > &ia,
199 const RCP<const Comm<int> > comm,
200 Zoltan2::MeshEntityType sourcetarget,
202 Teuchos::ArrayRCP<typename MeshAdapter<User>::offset_t> &offsets,
203 Teuchos::ArrayRCP<typename MeshAdapter<User>::gno_t> &adjacencyIds)
204{
205 typedef typename MeshAdapter<User>::gno_t gno_t;
206 typedef typename MeshAdapter<User>::lno_t lno_t;
207 typedef typename MeshAdapter<User>::offset_t offset_t;
208 typedef typename MeshAdapter<User>::node_t node_t;
209
210 typedef int nonzero_t; // adjacency matrix doesn't need scalar_t
211 typedef Tpetra::CrsMatrix<nonzero_t,lno_t,gno_t,node_t> sparse_matrix_type;
212
213 RCP<sparse_matrix_type> secondAdjs = get2ndAdjsMatFromAdjs(ia,comm,sourcetarget,through);
214
215 /* Allocate memory necessary for the adjacency */
216 size_t LocalNumIDs = ia->getLocalNumOf(sourcetarget);
217 lno_t *start = new lno_t [LocalNumIDs+1];
218
219 if (secondAdjs!=RCP<sparse_matrix_type>()) {
220
221 size_t nadj = 0;
222
223 gno_t const *Ids=NULL;
224 ia->getIDsViewOf(sourcetarget, Ids);
225
226 std::vector<gno_t> adj;
227
228 for (size_t localElement=0; localElement<LocalNumIDs; ++localElement){
229 start[localElement] = nadj;
230 const gno_t globalRow = Ids[localElement];
231 size_t NumEntries = secondAdjs->getNumEntriesInGlobalRow (globalRow);
232 typename sparse_matrix_type::nonconst_global_inds_host_view_type Indices("Indices", NumEntries);
233 typename sparse_matrix_type::nonconst_values_host_view_type Values("Values", NumEntries);
234 secondAdjs->getGlobalRowCopy (globalRow,Indices,Values,NumEntries);
235
236 for (size_t j = 0; j < NumEntries; ++j) {
237 if(globalRow != Indices[j]) {
238 adj.push_back(Indices[j]);
239 nadj++;;
240 }
241 }
242 }
243
244 Ids = NULL;
245 start[LocalNumIDs] = nadj;
246
247 gno_t *adj_ = new gno_t [nadj];
248
249 for (size_t i=0; i < nadj; i++) {
250 adj_[i] = adj[i];
251 }
252
253 offsets = arcp<offset_t>(start, 0, LocalNumIDs+1, true);
254 adjacencyIds = arcp<gno_t>(adj_, 0, nadj, true);
255 }
256 else {
257 // No adjacencies could be computed; return no edges and valid offsets array
258 for (size_t i = 0; i <= LocalNumIDs; i++)
259 start[i] = 0;
260
261 offsets = arcp<offset_t>(start, 0, LocalNumIDs+1, true);
262 adjacencyIds = Teuchos::null;
263 }
264
265 //return nadj;
266}
267
268}
269
270#endif
Defines the MeshAdapter interface.
InputTraits< User >::node_t node_t
InputTraits< User >::offset_t offset_t
InputTraits< User >::lno_t lno_t
InputTraits< User >::gno_t gno_t
map_t::local_ordinal_type lno_t
map_t::global_ordinal_type gno_t
Created by mbenlioglu on Aug 31, 2020.
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)
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)
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.