47 #include "Teko_BlockingEpetra.hpp"
49 #include "Epetra_IntVector.h"
50 #include "Epetra_LocalMap.h"
72 const MapPair buildSubMap(
const std::vector<int> &gid,
const Epetra_Comm &comm) {
73 Teuchos::RCP<Epetra_Map> gidMap = rcp(
new Epetra_Map(-1, gid.size(), &gid[0], 0, comm));
74 Teuchos::RCP<Epetra_Map> contigMap = rcp(
new Epetra_Map(-1, gid.size(), 0, comm));
76 return std::make_pair(gidMap, contigMap);
88 const ImExPair buildExportImport(
const Epetra_Map &baseMap,
const MapPair &maps) {
89 return std::make_pair(rcp(
new Epetra_Import(*maps.first, baseMap)),
90 rcp(
new Epetra_Export(*maps.first, baseMap)));
100 void buildSubVectors(
const std::vector<MapPair> &maps,
101 std::vector<RCP<Epetra_MultiVector> > &vectors,
int count) {
102 std::vector<MapPair>::const_iterator mapItr;
105 for (mapItr = maps.begin(); mapItr != maps.end(); mapItr++) {
107 RCP<Epetra_MultiVector> mv = rcp(
new Epetra_MultiVector(*(*mapItr).second, count));
108 vectors.push_back(mv);
119 void one2many(std::vector<RCP<Epetra_MultiVector> > &many,
const Epetra_MultiVector &single,
120 const std::vector<RCP<Epetra_Import> > &subImport) {
122 std::vector<RCP<Epetra_MultiVector> >::const_iterator vecItr;
123 std::vector<RCP<Epetra_Import> >::const_iterator impItr;
126 for (vecItr = many.begin(), impItr = subImport.begin(); vecItr != many.end();
127 ++vecItr, ++impItr) {
129 RCP<Epetra_MultiVector> destVec = *vecItr;
132 const Epetra_BlockMap &globalMap = (*impItr)->TargetMap();
135 Epetra_MultiVector importVector(View, globalMap, destVec->Values(), destVec->Stride(),
136 destVec->NumVectors());
139 importVector.Import(single, **impItr, Insert);
153 void many2one(Epetra_MultiVector &one,
const std::vector<RCP<const Epetra_MultiVector> > &many,
154 const std::vector<RCP<Epetra_Export> > &subExport) {
156 std::vector<RCP<const Epetra_MultiVector> >::const_iterator vecItr;
157 std::vector<RCP<Epetra_Export> >::const_iterator expItr;
160 for (vecItr = many.begin(), expItr = subExport.begin(); vecItr != many.end();
161 ++vecItr, ++expItr) {
163 RCP<const Epetra_MultiVector> srcVec = *vecItr;
166 const Epetra_BlockMap &globalMap = (*expItr)->SourceMap();
169 Epetra_MultiVector exportVector(View, globalMap, srcVec->Values(), srcVec->Stride(),
170 srcVec->NumVectors());
171 one.Export(exportVector, **expItr, Insert);
180 RCP<Epetra_IntVector> getSubBlockColumnGIDs(
const Epetra_CrsMatrix &A,
const MapPair &mapPair) {
181 RCP<const Epetra_Map> blkGIDMap = mapPair.first;
182 RCP<const Epetra_Map> blkContigMap = mapPair.second;
185 Epetra_IntVector rIndex(A.RowMap(),
true);
186 for (
int i = 0; i < A.NumMyRows(); i++) {
188 int lid = blkGIDMap->LID(A.GRID(i));
190 rIndex[i] = blkContigMap->GID(lid);
196 Epetra_Import
import(A.ColMap(), A.RowMap());
197 RCP<Epetra_IntVector> cIndex = rcp(
new Epetra_IntVector(A.ColMap(),
true));
198 cIndex->Import(rIndex,
import, Insert);
204 RCP<Epetra_CrsMatrix> buildSubBlock(
int i,
int j,
const Epetra_CrsMatrix &A,
205 const std::vector<MapPair> &subMaps) {
207 int numVarFamily = subMaps.size();
209 TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
210 TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
212 const Epetra_Map &gRowMap = *subMaps[i].first;
213 const Epetra_Map &rowMap = *subMaps[i].second;
214 const Epetra_Map &colMap = *subMaps[j].second;
216 const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A, subMaps[j]);
217 Epetra_IntVector &local2ContigGIDs = *plocal2ContigGIDs;
219 RCP<Epetra_CrsMatrix> mat = rcp(
new Epetra_CrsMatrix(Copy, rowMap, 0));
222 int numMyRows = rowMap.NumMyElements();
223 int maxNumEntries = A.MaxNumEntries();
226 std::vector<int> indices(maxNumEntries);
227 std::vector<double> values(maxNumEntries);
230 std::vector<int> colIndices(maxNumEntries);
231 std::vector<double> colValues(maxNumEntries);
235 for (
int localRow = 0; localRow < numMyRows; localRow++) {
237 int globalRow = gRowMap.GID(localRow);
238 int lid = A.RowMap().LID(globalRow);
239 int contigRow = rowMap.GID(localRow);
240 TEUCHOS_ASSERT(lid > -1);
242 int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
243 TEUCHOS_ASSERT(err == 0);
245 int numOwnedCols = 0;
246 for (
int localCol = 0; localCol < numEntries; localCol++) {
247 TEUCHOS_ASSERT(indices[localCol] > -1);
250 int gid = local2ContigGIDs[indices[localCol]];
251 if (gid == -1)
continue;
253 colIndices[numOwnedCols] = gid;
254 colValues[numOwnedCols] = values[localCol];
259 mat->InsertGlobalValues(contigRow, numOwnedCols, &colValues[0], &colIndices[0]);
263 mat->FillComplete(colMap, rowMap);
269 void rebuildSubBlock(
int i,
int j,
const Epetra_CrsMatrix &A,
const std::vector<MapPair> &subMaps,
270 Epetra_CrsMatrix &mat) {
272 int numVarFamily = subMaps.size();
274 TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
275 TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);
277 const Epetra_Map &gRowMap = *subMaps[i].first;
278 const Epetra_Map &rowMap = *subMaps[i].second;
280 const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A, subMaps[j]);
281 Epetra_IntVector &local2ContigGIDs = *plocal2ContigGIDs;
286 int numMyRows = rowMap.NumMyElements();
287 int maxNumEntries = A.MaxNumEntries();
290 std::vector<int> indices(maxNumEntries);
291 std::vector<double> values(maxNumEntries);
294 std::vector<int> colIndices(maxNumEntries);
295 std::vector<double> colValues(maxNumEntries);
299 for (
int localRow = 0; localRow < numMyRows; localRow++) {
301 int globalRow = gRowMap.GID(localRow);
302 int lid = A.RowMap().LID(globalRow);
303 int contigRow = rowMap.GID(localRow);
304 TEUCHOS_ASSERT(lid > -1);
306 int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
307 TEUCHOS_ASSERT(err == 0);
309 int numOwnedCols = 0;
310 for (
int localCol = 0; localCol < numEntries; localCol++) {
311 TEUCHOS_ASSERT(indices[localCol] > -1);
314 int gid = local2ContigGIDs[indices[localCol]];
315 if (gid == -1)
continue;
317 colIndices[numOwnedCols] = gid;
318 colValues[numOwnedCols] = values[localCol];
323 mat.SumIntoGlobalValues(contigRow, numOwnedCols, &colValues[0], &colIndices[0]);