Teko Version of the Day
Loading...
Searching...
No Matches
Teko_BlockingTpetra.cpp
1// @HEADER
2// *****************************************************************************
3// Teko: A package for block and physics based preconditioning
4//
5// Copyright 2010 NTESS and the Teko contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#include "Teko_BlockingTpetra.hpp"
11#include "Teko_Utilities.hpp"
12
13#include "Tpetra_Vector.hpp"
14#include "Tpetra_Map.hpp"
15
16using Teuchos::RCP;
17using Teuchos::rcp;
18
19namespace Teko {
20namespace TpetraHelpers {
21namespace Blocking {
22
36const MapPair buildSubMap(const std::vector<GO>& gid, const Teuchos::Comm<int>& comm) {
37 using GST = Tpetra::global_size_t;
38 const GST invalid = Teuchos::OrdinalTraits<GST>::invalid();
39 Teuchos::RCP<Tpetra::Map<LO, GO, NT> > gidMap = rcp(
40 new Tpetra::Map<LO, GO, NT>(invalid, Teuchos::ArrayView<const GO>(gid), 0, rcpFromRef(comm)));
41 Teuchos::RCP<Tpetra::Map<LO, GO, NT> > contigMap =
42 rcp(new Tpetra::Map<LO, GO, NT>(invalid, gid.size(), 0, rcpFromRef(comm)));
43
44 return std::make_pair(gidMap, contigMap);
45}
46
56const ImExPair buildExportImport(const Tpetra::Map<LO, GO, NT>& baseMap, const MapPair& maps) {
57 return std::make_pair(rcp(new Tpetra::Import<LO, GO, NT>(rcpFromRef(baseMap), maps.first)),
58 rcp(new Tpetra::Export<LO, GO, NT>(maps.first, rcpFromRef(baseMap))));
59}
60
68void buildSubVectors(const std::vector<MapPair>& maps,
69 std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& vectors, int count) {
70 std::vector<MapPair>::const_iterator mapItr;
71
72 // loop over all maps
73 for (mapItr = maps.begin(); mapItr != maps.end(); mapItr++) {
74 // add new elements to vectors
75 RCP<Tpetra::MultiVector<ST, LO, GO, NT> > mv =
76 rcp(new Tpetra::MultiVector<ST, LO, GO, NT>((*mapItr).second, count));
77 vectors.push_back(mv);
78 }
79}
80
88void one2many(std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
89 const Tpetra::MultiVector<ST, LO, GO, NT>& single,
90 const std::vector<RCP<Tpetra::Import<LO, GO, NT> > >& subImport) {
91 // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
92 std::vector<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
93 std::vector<RCP<Tpetra::Import<LO, GO, NT> > >::const_iterator impItr;
94
95 // using Importers fill the sub vectors from the mama vector
96 for (vecItr = many.begin(), impItr = subImport.begin(); vecItr != many.end();
97 ++vecItr, ++impItr) {
98 // for ease of access to the destination
99 RCP<Tpetra::MultiVector<ST, LO, GO, NT> > destVec = *vecItr;
100
101 // extract the map with global indicies from the current vector
102 const Tpetra::Map<LO, GO, NT>& globalMap = *(*impItr)->getTargetMap();
103
104 // build the import vector as a view on the destination
105 GO lda = destVec->getStride();
106 GO destSize = destVec->getGlobalLength() * destVec->getNumVectors();
107 std::vector<ST> destArray(destSize);
108 Teuchos::ArrayView<ST> destVals(destArray);
109 destVec->get1dCopy(destVals, lda);
110 Tpetra::MultiVector<ST, LO, GO, NT> importVector(rcpFromRef(globalMap), destVals, lda,
111 destVec->getNumVectors());
112
113 // perform the import
114 importVector.doImport(single, **impItr, Tpetra::INSERT);
115
116 Tpetra::Import<LO, GO, NT> importer(destVec->getMap(), destVec->getMap());
117 importVector.replaceMap(destVec->getMap());
118 destVec->doExport(importVector, importer, Tpetra::INSERT);
119 }
120}
121
132void many2one(Tpetra::MultiVector<ST, LO, GO, NT>& one,
133 const std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >& many,
134 const std::vector<RCP<Tpetra::Export<LO, GO, NT> > >& subExport) {
135 // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
136 std::vector<RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > >::const_iterator vecItr;
137 std::vector<RCP<Tpetra::Export<LO, GO, NT> > >::const_iterator expItr;
138
139 // using Exporters fill the empty vector from the sub-vectors
140 for (vecItr = many.begin(), expItr = subExport.begin(); vecItr != many.end();
141 ++vecItr, ++expItr) {
142 // for ease of access to the source
143 RCP<const Tpetra::MultiVector<ST, LO, GO, NT> > srcVec = *vecItr;
144
145 // extract the map with global indicies from the current vector
146 const Tpetra::Map<LO, GO, NT>& globalMap = *(*expItr)->getSourceMap();
147
148 // build the export vector as a view of the destination
149 GO lda = srcVec->getStride();
150 GO srcSize = srcVec->getGlobalLength() * srcVec->getNumVectors();
151 std::vector<ST> srcArray(srcSize);
152 Teuchos::ArrayView<ST> srcVals(srcArray);
153 srcVec->get1dCopy(srcVals, lda);
154 Tpetra::MultiVector<ST, LO, GO, NT> exportVector(rcpFromRef(globalMap), srcVals, lda,
155 srcVec->getNumVectors());
156
157 // perform the export
158 one.doExport(exportVector, **expItr, Tpetra::INSERT);
159 }
160}
161
167RCP<Tpetra::Vector<GO, LO, GO, NT> > getSubBlockColumnGIDs(
168 const Tpetra::CrsMatrix<ST, LO, GO, NT>& A, const MapPair& mapPair) {
169 RCP<const Tpetra::Map<LO, GO, NT> > blkGIDMap = mapPair.first;
170 RCP<const Tpetra::Map<LO, GO, NT> > blkContigMap = mapPair.second;
171
172 // fill index vector for rows
173 Tpetra::Vector<GO, LO, GO, NT> rIndex(A.getRowMap(), true);
174 for (size_t i = 0; i < A.getLocalNumRows(); i++) {
175 // LID is need to get to contiguous map
176 LO lid = blkGIDMap->getLocalElement(
177 A.getRowMap()->getGlobalElement(i)); // this LID makes me nervous
178 if (lid > -1)
179 rIndex.replaceLocalValue(i, blkContigMap->getGlobalElement(lid));
180 else
181 rIndex.replaceLocalValue(i, -1);
182 }
183
184 // get relavant column indices
185
186 Tpetra::Import<LO, GO, NT> import(A.getRowMap(), A.getColMap());
187 RCP<Tpetra::Vector<GO, LO, GO, NT> > cIndex =
188 rcp(new Tpetra::Vector<GO, LO, GO, NT>(A.getColMap(), true));
189 cIndex->doImport(rIndex, import, Tpetra::INSERT);
190
191 return cIndex;
192}
193
194// build a single subblock Epetra_CrsMatrix
195RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildSubBlock(
196 int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
197 const std::vector<MapPair>& subMaps,
198 Teuchos::RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs) {
199 // get the number of variables families
200 int numVarFamily = subMaps.size();
201
202 TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
203 TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
204 const RCP<const Tpetra::Map<LO, GO, NT> > gRowMap = subMaps[i].first; // new GIDs
205 const RCP<const Tpetra::Map<LO, GO, NT> > rowMap = subMaps[i].second; // contiguous GIDs
206 const RCP<const Tpetra::Map<LO, GO, NT> > colMap = subMaps[j].second;
207
208 if (!plocal2ContigGIDs) {
209 plocal2ContigGIDs = Blocking::getSubBlockColumnGIDs(*A, subMaps[j]);
210 }
211 Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
212
213 // get entry information
214 LO numMyRows = rowMap->getLocalNumElements();
215 LO maxNumEntries = A->getLocalMaxNumRowEntries();
216
217 // for extraction
218 auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
219 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
220 auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
221 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
222
223 // for insertion
224 std::vector<GO> colIndices(maxNumEntries);
225 std::vector<ST> colValues(maxNumEntries);
226
227 // for counting
228 std::vector<size_t> nEntriesPerRow(numMyRows);
229
230 const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
231
232 // insert each row into subblock
233 // Count the number of entries per row in the new matrix
234 for (LO localRow = 0; localRow < numMyRows; localRow++) {
235 size_t numEntries = invalid;
236 GO globalRow = gRowMap->getGlobalElement(localRow);
237 LO lid = A->getRowMap()->getLocalElement(globalRow);
238 TEUCHOS_ASSERT(lid > -1);
239
240 A->getLocalRowCopy(lid, indices, values, numEntries);
241
242 LO numOwnedCols = 0;
243 auto data = local2ContigGIDs.getLocalViewHost(Tpetra::Access::ReadOnly);
244 for (size_t localCol = 0; localCol < numEntries; localCol++) {
245 TEUCHOS_ASSERT(indices(localCol) > -1);
246
247 // if global id is not owned by this column
248
249 GO gid = data(indices[localCol], 0);
250 if (gid == -1) continue; // in contiguous row
251 numOwnedCols++;
252 }
253 nEntriesPerRow[localRow] = numOwnedCols;
254 }
255
256 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > mat = rcp(
257 new Tpetra::CrsMatrix<ST, LO, GO, NT>(rowMap, Teuchos::ArrayView<size_t>(nEntriesPerRow)));
258
259 // insert each row into subblock
260 // let FillComplete handle column distribution
261 for (LO localRow = 0; localRow < numMyRows; localRow++) {
262 size_t numEntries = invalid;
263 GO globalRow = gRowMap->getGlobalElement(localRow);
264 LO lid = A->getRowMap()->getLocalElement(globalRow);
265 GO contigRow = rowMap->getGlobalElement(localRow);
266 TEUCHOS_ASSERT(lid > -1);
267
268 A->getLocalRowCopy(lid, indices, values, numEntries);
269
270 LO numOwnedCols = 0;
271 auto data = local2ContigGIDs.getLocalViewHost(Tpetra::Access::ReadOnly);
272 for (size_t localCol = 0; localCol < numEntries; localCol++) {
273 TEUCHOS_ASSERT(indices(localCol) > -1);
274
275 // if global id is not owned by this column
276
277 GO gid = data(indices(localCol), 0);
278 if (gid == -1) continue; // in contiguous row
279
280 colIndices[numOwnedCols] = gid;
281 colValues[numOwnedCols] = values(localCol);
282 numOwnedCols++;
283 }
284
285 // insert it into the new matrix
286 mat->insertGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices.data(), numOwnedCols),
287 Teuchos::ArrayView<ST>(colValues.data(), numOwnedCols));
288 }
289
290 // fill it and automagically optimize the storage
291 mat->fillComplete(colMap, rowMap);
292
293 return mat;
294}
295
296// build a single subblock Epetra_CrsMatrix
297void rebuildSubBlock(int i, int j, const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> >& A,
298 const std::vector<MapPair>& subMaps, Tpetra::CrsMatrix<ST, LO, GO, NT>& mat,
299 Teuchos::RCP<Tpetra::Vector<GO, LO, GO, NT> > plocal2ContigGIDs) {
300 // get the number of variables families
301 int numVarFamily = subMaps.size();
302
303 TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
304 TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
305
306 const Tpetra::Map<LO, GO, NT>& gRowMap = *subMaps[i].first; // new GIDs
307 const Tpetra::Map<LO, GO, NT>& rowMap = *subMaps[i].second; // contiguous GIDs
308 const Tpetra::Map<LO, GO, NT>& colMap = *subMaps[j].second;
309
310 if (!plocal2ContigGIDs) {
311 plocal2ContigGIDs = Blocking::getSubBlockColumnGIDs(*A, subMaps[j]);
312 }
313 Tpetra::Vector<GO, LO, GO, NT>& local2ContigGIDs = *plocal2ContigGIDs;
314
315 mat.resumeFill();
316 mat.setAllToScalar(0.0);
317
318 // get entry information
319 LO numMyRows = rowMap.getLocalNumElements();
320 LO maxNumEntries = A->getLocalMaxNumRowEntries();
321
322 // for extraction
323 auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_local_inds_host_view_type(
324 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
325 auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
326 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), maxNumEntries);
327
328 // for insertion
329 std::vector<GO> colIndices(maxNumEntries);
330 std::vector<ST> colValues(maxNumEntries);
331
332 const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
333
334 // insert each row into subblock
335 // let FillComplete handle column distribution
336 auto A_rowmap = A->getRowMap();
337 for (LO localRow = 0; localRow < numMyRows; localRow++) {
338 size_t numEntries = invalid;
339 GO globalRow = gRowMap.getGlobalElement(localRow);
340 LO lid = A_rowmap->getLocalElement(globalRow);
341 GO contigRow = rowMap.getGlobalElement(localRow);
342 TEUCHOS_ASSERT(lid > -1);
343
344 A->getLocalRowCopy(lid, indices, values, numEntries);
345
346 LO numOwnedCols = 0;
347 auto data = local2ContigGIDs.getLocalViewHost(Tpetra::Access::ReadOnly);
348 for (size_t localCol = 0; localCol < numEntries; localCol++) {
349 TEUCHOS_ASSERT(indices(localCol) > -1);
350
351 // if global id is not owned by this column
352 GO gid = data(indices(localCol), 0);
353 if (gid == -1) continue; // in contiguous row
354
355 colIndices[numOwnedCols] = gid;
356 colValues[numOwnedCols] = values(localCol);
357 numOwnedCols++;
358 }
359
360 mat.sumIntoGlobalValues(contigRow, Teuchos::ArrayView<GO>(colIndices.data(), numOwnedCols),
361 Teuchos::ArrayView<ST>(colValues.data(), numOwnedCols));
362 }
363
364 mat.fillComplete(rcpFromRef(colMap), rcpFromRef(rowMap));
365}
366
367} // namespace Blocking
368} // namespace TpetraHelpers
369} // namespace Teko
Teko_Utilities.hpp
Generated on for Teko by doxygen 1.15.0