47#include "Teko_BlockingTpetra.hpp"
50#include "Tpetra_Vector.hpp"
51#include "Tpetra_Map.hpp"
57namespace TpetraHelpers {
73const MapPair buildSubMap(
const std::vector< GO > & gid,
const Teuchos::Comm<int> &comm)
75 using GST = Tpetra::global_size_t;
76 const GST invalid = Teuchos::OrdinalTraits<GST>::invalid();
77 Teuchos::RCP<Tpetra::Map<LO,GO,NT> > gidMap = rcp(
new Tpetra::Map<LO,GO,NT>(invalid,Teuchos::ArrayView<const GO>(gid),0,rcpFromRef(comm)));
78 Teuchos::RCP<Tpetra::Map<LO,GO,NT> > contigMap = rcp(
new Tpetra::Map<LO,GO,NT>(invalid,gid.size(),0,rcpFromRef(comm)));
80 return std::make_pair(gidMap,contigMap);
92const ImExPair buildExportImport(
const Tpetra::Map<LO,GO,NT> & baseMap,
const MapPair & maps)
94 return std::make_pair(rcp(
new Tpetra::Import<LO,GO,NT>(rcpFromRef(baseMap),maps.first)),
95 rcp(
new Tpetra::Export<LO,GO,NT>(maps.first,rcpFromRef(baseMap))));
105void buildSubVectors(
const std::vector<MapPair> & maps,
106 std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > & vectors,
int count)
108 std::vector<MapPair>::const_iterator mapItr;
111 for(mapItr=maps.begin();mapItr!=maps.end();mapItr++) {
113 RCP<Tpetra::MultiVector<ST,LO,GO,NT> > mv = rcp(
new Tpetra::MultiVector<ST,LO,GO,NT>((*mapItr).second,count));
114 vectors.push_back(mv);
124void one2many(std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > & many,
const Tpetra::MultiVector<ST,LO,GO,NT> & single,
125 const std::vector<RCP<Tpetra::Import<LO,GO,NT> > > & subImport)
128 std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
129 std::vector<RCP<Tpetra::Import<LO,GO,NT> > >::const_iterator impItr;
132 for(vecItr=many.begin(),impItr=subImport.begin();
133 vecItr!=many.end();++vecItr,++impItr) {
135 RCP<Tpetra::MultiVector<ST,LO,GO,NT> > destVec = *vecItr;
138 const Tpetra::Map<LO,GO,NT> & globalMap = *(*impItr)->getTargetMap();
141 GO lda = destVec->getStride();
142 GO destSize = destVec->getGlobalLength()*destVec->getNumVectors();
143 std::vector<ST> destArray(destSize);
144 Teuchos::ArrayView<ST> destVals(destArray);
145 destVec->get1dCopy(destVals,lda);
146 Tpetra::MultiVector<ST,LO,GO,NT> importVector(rcpFromRef(globalMap),destVals,lda,destVec->getNumVectors());
149 importVector.doImport(single,**impItr,Tpetra::INSERT);
151 Tpetra::Import<LO,GO,NT> importer(destVec->getMap(),destVec->getMap());
152 importVector.replaceMap(destVec->getMap());
153 destVec->doExport(importVector,importer,Tpetra::INSERT);
166void many2one(Tpetra::MultiVector<ST,LO,GO,NT> & one,
const std::vector<RCP<
const Tpetra::MultiVector<ST,LO,GO,NT> > > & many,
167 const std::vector<RCP<Tpetra::Export<LO,GO,NT> > > & subExport)
170 std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > >::const_iterator vecItr;
171 std::vector<RCP<Tpetra::Export<LO,GO,NT> > >::const_iterator expItr;
174 for(vecItr=many.begin(),expItr=subExport.begin();
175 vecItr!=many.end();++vecItr,++expItr) {
178 RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > srcVec = *vecItr;
181 const Tpetra::Map<LO,GO,NT> & globalMap = *(*expItr)->getSourceMap();
184 GO lda = srcVec->getStride();
185 GO srcSize = srcVec->getGlobalLength()*srcVec->getNumVectors();
186 std::vector<ST> srcArray(srcSize);
187 Teuchos::ArrayView<ST> srcVals(srcArray);
188 srcVec->get1dCopy(srcVals,lda);
189 Tpetra::MultiVector<ST,LO,GO,NT> exportVector(rcpFromRef(globalMap),srcVals,lda,srcVec->getNumVectors());
192 one.doExport(exportVector,**expItr,Tpetra::INSERT);
201RCP<Tpetra::Vector<GO,LO,GO,NT> > getSubBlockColumnGIDs(
const Tpetra::CrsMatrix<ST,LO,GO,NT> & A,
const MapPair & mapPair)
203 RCP<const Tpetra::Map<LO,GO,NT> > blkGIDMap = mapPair.first;
204 RCP<const Tpetra::Map<LO,GO,NT> > blkContigMap = mapPair.second;
207 Tpetra::Vector<GO,LO,GO,NT> rIndex(A.getRowMap(),
true);
208 for(
size_t i=0;i<A.getLocalNumRows();i++) {
210 LO lid = blkGIDMap->getLocalElement(A.getRowMap()->getGlobalElement(i));
212 rIndex.replaceLocalValue(i,blkContigMap->getGlobalElement(lid));
214 rIndex.replaceLocalValue(i,-1);
219 Tpetra::Import<LO,GO,NT>
import(A.getRowMap(),A.getColMap());
220 RCP<Tpetra::Vector<GO,LO,GO,NT> > cIndex = rcp(
new Tpetra::Vector<GO,LO,GO,NT>(A.getColMap(),
true));
221 cIndex->doImport(rIndex,
import,Tpetra::INSERT);
227RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildSubBlock(
int i,
int j,
const RCP<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >& A,
const std::vector<MapPair> & subMaps)
230 int numVarFamily = subMaps.size();
232 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
233 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
234 const RCP<const Tpetra::Map<LO,GO,NT> > gRowMap = subMaps[i].first;
235 const RCP<const Tpetra::Map<LO,GO,NT> > rowMap = subMaps[i].second;
236 const RCP<const Tpetra::Map<LO,GO,NT> > colMap = subMaps[j].second;
238 const RCP<Tpetra::Vector<GO,LO,GO,NT> > plocal2ContigGIDs = getSubBlockColumnGIDs(*A,subMaps[j]);
239 Tpetra::Vector<GO,LO,GO,NT> & local2ContigGIDs = *plocal2ContigGIDs;
243 LO numMyRows = rowMap->getLocalNumElements();
244 LO maxNumEntries = A->getLocalMaxNumRowEntries();
247 auto indices =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_local_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),maxNumEntries);
248 auto values =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),maxNumEntries);
251 std::vector<GO> colIndices(maxNumEntries);
252 std::vector<ST> colValues(maxNumEntries);
255 std::vector<size_t> nEntriesPerRow(numMyRows);
257 const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
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 TEUCHOS_ASSERT(lid>-1);
267 A->getLocalRowCopy(lid, indices, values, numEntries);
270 for(
size_t localCol=0;localCol<numEntries;localCol++) {
271 TEUCHOS_ASSERT(indices(localCol)>-1);
275 GO gid = local2ContigGIDs.getData()[indices[localCol]];
276 if(gid==-1)
continue;
279 nEntriesPerRow[localRow] = numOwnedCols;
282 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > mat = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(rowMap,Teuchos::ArrayView<size_t>(nEntriesPerRow)));
286 for(LO localRow=0;localRow<numMyRows;localRow++) {
287 size_t numEntries = invalid;
288 GO globalRow = gRowMap->getGlobalElement(localRow);
289 LO lid = A->getRowMap()->getLocalElement(globalRow);
290 GO contigRow = rowMap->getGlobalElement(localRow);
291 TEUCHOS_ASSERT(lid>-1);
293 A->getLocalRowCopy(lid, indices, values, numEntries);
296 for(
size_t localCol=0;localCol<numEntries;localCol++) {
297 TEUCHOS_ASSERT(indices(localCol)>-1);
301 GO gid = local2ContigGIDs.getData()[indices(localCol)];
302 if(gid==-1)
continue;
304 colIndices[numOwnedCols] = gid;
305 colValues[numOwnedCols] = values(localCol);
310 colIndices.resize(numOwnedCols);
311 colValues.resize(numOwnedCols);
312 mat->insertGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
313 colIndices.resize(maxNumEntries);
314 colValues.resize(maxNumEntries);
318 mat->fillComplete(colMap,rowMap);
324void rebuildSubBlock(
int i,
int j,
const RCP<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >& A,
const std::vector<MapPair> & subMaps,Tpetra::CrsMatrix<ST,LO,GO,NT> & mat)
327 int numVarFamily = subMaps.size();
329 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
330 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
332 const Tpetra::Map<LO,GO,NT> & gRowMap = *subMaps[i].first;
333 const Tpetra::Map<LO,GO,NT> & rowMap = *subMaps[i].second;
334 const Tpetra::Map<LO,GO,NT> & colMap = *subMaps[j].second;
336 const RCP<Tpetra::Vector<GO,LO,GO,NT> > plocal2ContigGIDs = getSubBlockColumnGIDs(*A,subMaps[j]);
337 Tpetra::Vector<GO,LO,GO,NT> & local2ContigGIDs = *plocal2ContigGIDs;
340 mat.setAllToScalar(0.0);
343 LO numMyRows = rowMap.getLocalNumElements();
344 LO maxNumEntries = A->getLocalMaxNumRowEntries();
347 auto indices =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_local_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),maxNumEntries);
348 auto values =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),maxNumEntries);
351 std::vector<GO> colIndices(maxNumEntries);
352 std::vector<ST> colValues(maxNumEntries);
354 const size_t invalid = Teuchos::OrdinalTraits<size_t>::invalid();
358 for(LO localRow=0;localRow<numMyRows;localRow++) {
359 size_t numEntries = invalid;
360 GO globalRow = gRowMap.getGlobalElement(localRow);
361 LO lid = A->getRowMap()->getLocalElement(globalRow);
362 GO contigRow = rowMap.getGlobalElement(localRow);
363 TEUCHOS_ASSERT(lid>-1);
365 A->getLocalRowCopy(lid, indices, values, numEntries);
368 for(
size_t localCol=0;localCol<numEntries;localCol++) {
369 TEUCHOS_ASSERT(indices(localCol)>-1);
372 GO gid = local2ContigGIDs.getData()[indices(localCol)];
373 if(gid==-1)
continue;
375 colIndices[numOwnedCols] = gid;
376 colValues[numOwnedCols] = values(localCol);
381 colIndices.resize(numOwnedCols);
382 colValues.resize(numOwnedCols);
383 mat.sumIntoGlobalValues(contigRow,Teuchos::ArrayView<GO>(colIndices),Teuchos::ArrayView<ST>(colValues));
384 colIndices.resize(maxNumEntries);
385 colValues.resize(maxNumEntries);
387 mat.fillComplete(rcpFromRef(colMap),rcpFromRef(rowMap));