Teko Version of the Day
Loading...
Searching...
No Matches
Teko_BlockingEpetra.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_BlockingEpetra.hpp"
48
49#include "Epetra_IntVector.h"
50#include "Epetra_LocalMap.h"
51
52using Teuchos::RCP;
53using Teuchos::rcp;
54
55namespace Teko {
56namespace Epetra {
57namespace Blocking {
58
72const MapPair buildSubMap(const std::vector< int > & gid, const Epetra_Comm &comm)
73{
74 Teuchos::RCP<Epetra_Map> gidMap = rcp(new Epetra_Map(-1,gid.size(),&gid[0],0,comm));
75 Teuchos::RCP<Epetra_Map> contigMap = rcp(new Epetra_Map(-1,gid.size(),0,comm));
76
77 return std::make_pair(gidMap,contigMap);
78}
79
89const ImExPair buildExportImport(const Epetra_Map & baseMap,const MapPair & maps)
90{
91 return std::make_pair(rcp(new Epetra_Import(*maps.first,baseMap)),
92 rcp(new Epetra_Export(*maps.first,baseMap)));
93}
94
102void buildSubVectors(const std::vector<MapPair> & maps,
103 std::vector<RCP<Epetra_MultiVector> > & vectors,int count)
104{
105 std::vector<MapPair>::const_iterator mapItr;
106
107 // loop over all maps
108 for(mapItr=maps.begin();mapItr!=maps.end();mapItr++) {
109 // add new elements to vectors
110 RCP<Epetra_MultiVector> mv = rcp(new Epetra_MultiVector(*(*mapItr).second,count));
111 vectors.push_back(mv);
112 }
113}
114
121void one2many(std::vector<RCP<Epetra_MultiVector> > & many, const Epetra_MultiVector & single,
122 const std::vector<RCP<Epetra_Import> > & subImport)
123{
124 // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
125 std::vector<RCP<Epetra_MultiVector> >::const_iterator vecItr;
126 std::vector<RCP<Epetra_Import> >::const_iterator impItr;
127
128 // using Importers fill the sub vectors from the mama vector
129 for(vecItr=many.begin(),impItr=subImport.begin();
130 vecItr!=many.end();++vecItr,++impItr) {
131 // for ease of access to the destination
132 RCP<Epetra_MultiVector> destVec = *vecItr;
133
134 // extract the map with global indicies from the current vector
135 const Epetra_BlockMap & globalMap = (*impItr)->TargetMap();
136
137 // build the import vector as a view on the destination
138 Epetra_MultiVector importVector(View,globalMap,destVec->Values(),destVec->Stride(),destVec->NumVectors());
139
140 // perform the import
141 importVector.Import(single,**impItr,Insert);
142 }
143}
144
154void many2one(Epetra_MultiVector & one, const std::vector<RCP<const Epetra_MultiVector> > & many,
155 const std::vector<RCP<Epetra_Export> > & subExport)
156{
157 // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
158 std::vector<RCP<const Epetra_MultiVector> >::const_iterator vecItr;
159 std::vector<RCP<Epetra_Export> >::const_iterator expItr;
160
161 // using Exporters fill the empty vector from the sub-vectors
162 for(vecItr=many.begin(),expItr=subExport.begin();
163 vecItr!=many.end();++vecItr,++expItr) {
164
165 // for ease of access to the source
166 RCP<const Epetra_MultiVector> srcVec = *vecItr;
167
168 // extract the map with global indicies from the current vector
169 const Epetra_BlockMap & globalMap = (*expItr)->SourceMap();
170
171 // build the export vector as a view of the destination
172 Epetra_MultiVector exportVector(View,globalMap,srcVec->Values(),srcVec->Stride(),srcVec->NumVectors());
173 one.Export(exportVector,**expItr,Insert);
174 }
175}
176
182RCP<Epetra_IntVector> getSubBlockColumnGIDs(const Epetra_CrsMatrix & A,const MapPair & mapPair)
183{
184 RCP<const Epetra_Map> blkGIDMap = mapPair.first;
185 RCP<const Epetra_Map> blkContigMap = mapPair.second;
186
187 // fill index vector for rows
188 Epetra_IntVector rIndex(A.RowMap(),true);
189 for(int i=0;i<A.NumMyRows();i++) {
190 // LID is need to get to contiguous map
191 int lid = blkGIDMap->LID(A.GRID(i)); // this LID makes me nervous
192 if(lid>-1)
193 rIndex[i] = blkContigMap->GID(lid);
194 else
195 rIndex[i] = -1;
196 }
197
198 // get relavant column indices
199 Epetra_Import import(A.ColMap(),A.RowMap());
200 RCP<Epetra_IntVector> cIndex = rcp(new Epetra_IntVector(A.ColMap(),true));
201 cIndex->Import(rIndex,import,Insert);
202
203 return cIndex;
204}
205
206// build a single subblock Epetra_CrsMatrix
207RCP<Epetra_CrsMatrix> buildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<MapPair> & subMaps)
208{
209 // get the number of variables families
210 int numVarFamily = subMaps.size();
211
212 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
213 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
214
215 const Epetra_Map & gRowMap = *subMaps[i].first; // new GIDs
216 const Epetra_Map & rowMap = *subMaps[i].second; // contiguous GIDs
217 const Epetra_Map & colMap = *subMaps[j].second;
218
219 const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A,subMaps[j]);
220 Epetra_IntVector & local2ContigGIDs = *plocal2ContigGIDs;
221
222 RCP<Epetra_CrsMatrix> mat = rcp(new Epetra_CrsMatrix(Copy,rowMap,0));
223
224 // get entry information
225 int numMyRows = rowMap.NumMyElements();
226 int maxNumEntries = A.MaxNumEntries();
227
228 // for extraction
229 std::vector<int> indices(maxNumEntries);
230 std::vector<double> values(maxNumEntries);
231
232 // for insertion
233 std::vector<int> colIndices(maxNumEntries);
234 std::vector<double> colValues(maxNumEntries);
235
236 // insert each row into subblock
237 // let FillComplete handle column distribution
238 for(int localRow=0;localRow<numMyRows;localRow++) {
239 int numEntries = -1;
240 int globalRow = gRowMap.GID(localRow);
241 int lid = A.RowMap().LID(globalRow);
242 int contigRow = rowMap.GID(localRow);
243 TEUCHOS_ASSERT(lid>-1);
244
245 int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
246 TEUCHOS_ASSERT(err==0);
247
248 int numOwnedCols = 0;
249 for(int localCol=0;localCol<numEntries;localCol++) {
250 TEUCHOS_ASSERT(indices[localCol]>-1);
251
252 // if global id is not owned by this column
253 int gid = local2ContigGIDs[indices[localCol]];
254 if(gid==-1) continue; // in contiguous row
255
256 colIndices[numOwnedCols] = gid;
257 colValues[numOwnedCols] = values[localCol];
258 numOwnedCols++;
259 }
260
261 // insert it into the new matrix
262 mat->InsertGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
263 }
264
265 // fill it and automagically optimize the storage
266 mat->FillComplete(colMap,rowMap);
267
268 return mat;
269}
270
271// build a single subblock Epetra_CrsMatrix
272void rebuildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<MapPair> & subMaps,Epetra_CrsMatrix & mat)
273{
274 // get the number of variables families
275 int numVarFamily = subMaps.size();
276
277 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
278 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
279
280 const Epetra_Map & gRowMap = *subMaps[i].first; // new GIDs
281 const Epetra_Map & rowMap = *subMaps[i].second; // contiguous GIDs
282
283 const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A,subMaps[j]);
284 Epetra_IntVector & local2ContigGIDs = *plocal2ContigGIDs;
285
286 mat.PutScalar(0.0);
287
288 // get entry information
289 int numMyRows = rowMap.NumMyElements();
290 int maxNumEntries = A.MaxNumEntries();
291
292 // for extraction
293 std::vector<int> indices(maxNumEntries);
294 std::vector<double> values(maxNumEntries);
295
296 // for insertion
297 std::vector<int> colIndices(maxNumEntries);
298 std::vector<double> colValues(maxNumEntries);
299
300 // insert each row into subblock
301 // let FillComplete handle column distribution
302 for(int localRow=0;localRow<numMyRows;localRow++) {
303 int numEntries = -1;
304 int globalRow = gRowMap.GID(localRow);
305 int lid = A.RowMap().LID(globalRow);
306 int contigRow = rowMap.GID(localRow);
307 TEUCHOS_ASSERT(lid>-1);
308
309 int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
310 TEUCHOS_ASSERT(err==0);
311
312 int numOwnedCols = 0;
313 for(int localCol=0;localCol<numEntries;localCol++) {
314 TEUCHOS_ASSERT(indices[localCol]>-1);
315
316 // if global id is not owned by this column
317 int gid = local2ContigGIDs[indices[localCol]];
318 if(gid==-1) continue; // in contiguous row
319
320 colIndices[numOwnedCols] = gid;
321 colValues[numOwnedCols] = values[localCol];
322 numOwnedCols++;
323 }
324
325 // insert it into the new matrix
326 mat.SumIntoGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
327 }
328}
329
330} // end Blocking
331} // end Epetra
332} // end Teko