Teko Version of the Day
Loading...
Searching...
No Matches
Teko_InterlacedEpetra.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_InterlacedEpetra.hpp"
48
49#include <vector>
50
51using Teuchos::RCP;
52using Teuchos::rcp;
53
54namespace Teko {
55namespace Epetra {
56namespace Strided {
57
58// this assumes that there are numGlobals with numVars each interlaced
59// i.e. for numVars = 2 (u,v) then the vector is
60// [u_0,v_0,u_1,v_1,u_2,v_2, ..., u_(numGlobals-1),v_(numGlobals-1)]
61void buildSubMaps(int numGlobals,int numVars,const Epetra_Comm & comm,std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps)
62{
63 std::vector<int> vars;
64
65 // build vector describing the sub maps
66 for(int i=0;i<numVars;i++) vars.push_back(1);
67
68 // build all the submaps
69 buildSubMaps(numGlobals,vars,comm,subMaps);
70}
71
72// build maps to make other conversions
73void buildSubMaps(const Epetra_Map & globalMap,const std::vector<int> & vars,const Epetra_Comm & comm,
74 std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
75{
76 buildSubMaps(globalMap.NumGlobalElements(),globalMap.NumMyElements(),globalMap.MinMyGID(),
77 vars,comm,subMaps);
78}
79
80// build maps to make other conversions
81void buildSubMaps(int numGlobals,const std::vector<int> & vars,const Epetra_Comm & comm,std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
82{
83 std::vector<int>::const_iterator varItr;
84
85 // compute total number of variables
86 int numGlobalVars = 0;
87 for(varItr=vars.begin();varItr!=vars.end();++varItr)
88 numGlobalVars += *varItr;
89
90 // must be an even number of globals
91 TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
92
93 Epetra_Map sampleMap(numGlobals/numGlobalVars,0,comm);
94
95 buildSubMaps(numGlobals,numGlobalVars*sampleMap.NumMyElements(),numGlobalVars*sampleMap.MinMyGID(),vars,comm,subMaps);
96}
97
98// build maps to make other conversions
99void buildSubMaps(int numGlobals,int numMyElements,int minMyGID,const std::vector<int> & vars,const Epetra_Comm & comm,
100 std::vector<std::pair<int,Teuchos::RCP<Epetra_Map> > > & subMaps)
101{
102 std::vector<int>::const_iterator varItr;
103
104 // compute total number of variables
105 int numGlobalVars = 0;
106 for(varItr=vars.begin();varItr!=vars.end();++varItr)
107 numGlobalVars += *varItr;
108
109 // must be an even number of globals
110 TEUCHOS_ASSERT((numGlobals%numGlobalVars)==0);
111 TEUCHOS_ASSERT((numMyElements%numGlobalVars)==0);
112 TEUCHOS_ASSERT((minMyGID%numGlobalVars)==0);
113
114 int numBlocks = numMyElements/numGlobalVars;
115 int minBlockID = minMyGID/numGlobalVars;
116
117 subMaps.clear();
118
119 // index into local block in strided map
120 int blockOffset = 0;
121 for(varItr=vars.begin();varItr!=vars.end();++varItr) {
122 int numLocalVars = *varItr;
123 int numAllElmts = numLocalVars*numGlobals/numGlobalVars;
124 int numMyElmts = numLocalVars * numBlocks;
125
126 // create global arrays describing the as of yet uncreated maps
127 std::vector<int> subGlobals;
128 std::vector<int> contigGlobals; // the contiguous globals
129
130 // loop over each block of variables
131 int count = 0;
132 for(int blockNum=0;blockNum<numBlocks;blockNum++) {
133
134 // loop over each local variable in the block
135 for(int local=0;local<numLocalVars;++local) {
136 // global block number = minGID+blockNum
137 // block begin global id = numGlobalVars*(minGID+blockNum)
138 // global id block offset = blockOffset+local
139 subGlobals.push_back((minBlockID+blockNum)*numGlobalVars+blockOffset+local);
140
141 // also build the contiguous IDs
142 contigGlobals.push_back(numLocalVars*minBlockID+count);
143 count++;
144 }
145 }
146
147 // sanity check
148 assert((unsigned int) numMyElmts==subGlobals.size());
149
150 // create the map with contiguous elements and the map with global elements
151 RCP<Epetra_Map> subMap = rcp(new Epetra_Map(numAllElmts,numMyElmts,&subGlobals[0],0,comm));
152 RCP<Epetra_Map> contigMap = rcp(new Epetra_Map(numAllElmts,numMyElmts,&contigGlobals[0],0,comm));
153
154 Teuchos::set_extra_data(contigMap,"contigMap",Teuchos::inOutArg(subMap));
155 subMaps.push_back(std::make_pair(numLocalVars,subMap));
156
157 // update the block offset
158 blockOffset += numLocalVars;
159 }
160}
161
162void buildExportImport(const Epetra_Map & baseMap, const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,
163 std::vector<RCP<Epetra_Export> > & subExport,
164 std::vector<RCP<Epetra_Import> > & subImport)
165{
166 std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
167
168 // build importers and exporters
169 for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
170 // exctract basic map
171 const Epetra_Map & map = *(mapItr->second);
172
173 // add new elements to vectors
174 subImport.push_back(rcp(new Epetra_Import(map,baseMap)));
175 subExport.push_back(rcp(new Epetra_Export(map,baseMap)));
176 }
177}
178
179void buildSubVectors(const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,std::vector<RCP<Epetra_MultiVector> > & subVectors,int count)
180{
181 std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
182
183 // build vectors
184 for(mapItr=subMaps.begin();mapItr!=subMaps.end();++mapItr) {
185 // exctract basic map
186 const Epetra_Map & map = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(mapItr->second,"contigMap"));
187
188 // add new elements to vectors
189 RCP<Epetra_MultiVector> mv = rcp(new Epetra_MultiVector(map,count));
190 Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(mv));
191 subVectors.push_back(mv);
192 }
193}
194
195void associateSubVectors(const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,std::vector<RCP<const Epetra_MultiVector> > & subVectors)
196{
197 std::vector<std::pair<int,RCP<Epetra_Map> > >::const_iterator mapItr;
198 std::vector<RCP<const Epetra_MultiVector> >::iterator vecItr;
199
200 TEUCHOS_ASSERT(subMaps.size()==subVectors.size());
201
202 // associate the sub vectors with the subMaps
203 for(mapItr=subMaps.begin(),vecItr=subVectors.begin();mapItr!=subMaps.end();++mapItr,++vecItr)
204 Teuchos::set_extra_data(mapItr->second,"globalMap",Teuchos::inOutArg(*vecItr));
205}
206
207// build a single subblock Epetra_CrsMatrix
208RCP<Epetra_CrsMatrix> buildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps)
209{
210 // get the number of variables families
211 int numVarFamily = subMaps.size();
212
213 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
214 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
215
216 const Epetra_Map & gRowMap = *subMaps[i].second;
217 const Epetra_Map & rowMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[i].second,"contigMap");
218 const Epetra_Map & colMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[j].second,"contigMap");
219 int colFamilyCnt = subMaps[j].first;
220
221 // compute the number of global variables
222 // and the row and column block offset
223 int numGlobalVars = 0;
224 int rowBlockOffset = 0;
225 int colBlockOffset = 0;
226 for(int k=0;k<numVarFamily;k++) {
227 numGlobalVars += subMaps[k].first;
228
229 // compute block offsets
230 if(k<i) rowBlockOffset += subMaps[k].first;
231 if(k<j) colBlockOffset += subMaps[k].first;
232 }
233
234 // copy all global rows to here
235 Epetra_Import import(gRowMap,A.RowMap());
236 Epetra_CrsMatrix localA(Copy,gRowMap,0);
237 localA.Import(A,import,Insert);
238
239 RCP<Epetra_CrsMatrix> mat = rcp(new Epetra_CrsMatrix(Copy,rowMap,0));
240
241 // get entry information
242 int numMyRows = rowMap.NumMyElements();
243 int maxNumEntries = A.GlobalMaxNumEntries();
244
245 // for extraction
246 std::vector<int> indices(maxNumEntries);
247 std::vector<double> values(maxNumEntries);
248
249 // for insertion
250 std::vector<int> colIndices(maxNumEntries);
251 std::vector<double> colValues(maxNumEntries);
252
253 // insert each row into subblock
254 // let FillComplete handle column distribution
255 for(int localRow=0;localRow<numMyRows;localRow++) {
256 int numEntries = -1;
257 int globalRow = gRowMap.GID(localRow);
258 int contigRow = rowMap.GID(localRow);
259
260 TEUCHOS_ASSERT(globalRow>=0);
261 TEUCHOS_ASSERT(contigRow>=0);
262
263 // extract a global row copy
264 int err = localA.ExtractGlobalRowCopy(globalRow, maxNumEntries, numEntries, &values[0], &indices[0]);
265 TEUCHOS_ASSERT(err==0);
266
267 int numOwnedCols = 0;
268 for(int localCol=0;localCol<numEntries;localCol++) {
269 int globalCol = indices[localCol];
270
271 // determinate which block this column ID is in
272 int block = globalCol / numGlobalVars;
273
274 bool inFamily = true;
275
276 // test the beginning of the block
277 inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
278 inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
279
280 // is this column in the variable family
281 if(inFamily) {
282 int familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
283
284 colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
285 colValues[numOwnedCols] = values[localCol];
286
287 numOwnedCols++;
288 }
289 }
290
291 // insert it into the new matrix
292 mat->InsertGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
293 }
294
295 // fill it and automagically optimize the storage
296 mat->FillComplete(colMap,rowMap);
297
298 return mat;
299}
300
301// rebuild a single subblock Epetra_CrsMatrix
302void rebuildSubBlock(int i,int j,const Epetra_CrsMatrix & A,const std::vector<std::pair<int,RCP<Epetra_Map> > > & subMaps,Epetra_CrsMatrix & mat)
303{
304 // get the number of variables families
305 int numVarFamily = subMaps.size();
306
307 TEUCHOS_ASSERT(i>=0 && i<numVarFamily);
308 TEUCHOS_ASSERT(j>=0 && j<numVarFamily);
309 TEUCHOS_ASSERT(mat.Filled());
310
311 const Epetra_Map & gRowMap = *subMaps[i].second;
312 const Epetra_Map & rowMap = *Teuchos::get_extra_data<RCP<Epetra_Map> >(subMaps[i].second,"contigMap");
313 int colFamilyCnt = subMaps[j].first;
314
315 // compute the number of global variables
316 // and the row and column block offset
317 int numGlobalVars = 0;
318 int rowBlockOffset = 0;
319 int colBlockOffset = 0;
320 for(int k=0;k<numVarFamily;k++) {
321 numGlobalVars += subMaps[k].first;
322
323 // compute block offsets
324 if(k<i) rowBlockOffset += subMaps[k].first;
325 if(k<j) colBlockOffset += subMaps[k].first;
326 }
327
328 // copy all global rows to here
329 Epetra_Import import(gRowMap,A.RowMap());
330 Epetra_CrsMatrix localA(Copy,gRowMap,0);
331 localA.Import(A,import,Insert);
332
333 // clear out the old matrix
334 mat.PutScalar(0.0);
335
336 // get entry information
337 int numMyRows = rowMap.NumMyElements();
338 int maxNumEntries = A.GlobalMaxNumEntries();
339
340 // for extraction
341 std::vector<int> indices(maxNumEntries);
342 std::vector<double> values(maxNumEntries);
343
344 // for insertion
345 std::vector<int> colIndices(maxNumEntries);
346 std::vector<double> colValues(maxNumEntries);
347
348 // insert each row into subblock
349 // let FillComplete handle column distribution
350 for(int localRow=0;localRow<numMyRows;localRow++) {
351 int numEntries = -1;
352 int globalRow = gRowMap.GID(localRow);
353 int contigRow = rowMap.GID(localRow);
354
355 TEUCHOS_ASSERT(globalRow>=0);
356 TEUCHOS_ASSERT(contigRow>=0);
357
358 // extract a global row copy
359 int err = localA.ExtractGlobalRowCopy(globalRow, maxNumEntries, numEntries, &values[0], &indices[0]);
360 TEUCHOS_ASSERT(err==0);
361
362 int numOwnedCols = 0;
363 for(int localCol=0;localCol<numEntries;localCol++) {
364 int globalCol = indices[localCol];
365
366 // determinate which block this column ID is in
367 int block = globalCol / numGlobalVars;
368
369 bool inFamily = true;
370
371 // test the beginning of the block
372 inFamily &= (block*numGlobalVars+colBlockOffset <= globalCol);
373 inFamily &= ((block*numGlobalVars+colBlockOffset+colFamilyCnt) > globalCol);
374
375 // is this column in the variable family
376 if(inFamily) {
377 int familyOffset = globalCol-(block*numGlobalVars+colBlockOffset);
378
379 colIndices[numOwnedCols] = block*colFamilyCnt + familyOffset;
380 colValues[numOwnedCols] = values[localCol];
381
382 numOwnedCols++;
383 }
384 }
385
386 // insert it into the new matrix
387 mat.SumIntoGlobalValues(contigRow,numOwnedCols,&colValues[0],&colIndices[0]);
388 }
389}
390
391
392// collect subvectors into a single global vector
393void many2one(Epetra_MultiVector & one, const std::vector<RCP<const Epetra_MultiVector> > & many,
394 const std::vector<RCP<Epetra_Export> > & subExport)
395{
396 // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
397 std::vector<RCP<const Epetra_MultiVector> >::const_iterator vecItr;
398 std::vector<RCP<Epetra_Export> >::const_iterator expItr;
399
400 // using Exporters fill the empty vector from the sub-vectors
401 for(vecItr=many.begin(),expItr=subExport.begin();
402 vecItr!=many.end();++vecItr,++expItr) {
403
404 // for ease of access to the source
405 RCP<const Epetra_MultiVector> srcVec = *vecItr;
406
407 // extract the map with global indicies from the current vector
408 const Epetra_Map & globalMap = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(srcVec,"globalMap"));
409
410 // build the export vector as a view of the destination
411 Epetra_MultiVector exportVector(View,globalMap,srcVec->Values(),srcVec->Stride(),srcVec->NumVectors());
412 one.Export(exportVector,**expItr,Insert);
413 }
414}
415
416// distribute one global vector into a many subvectors
417void one2many(std::vector<RCP<Epetra_MultiVector> > & many,const Epetra_MultiVector & single,
418 const std::vector<RCP<Epetra_Import> > & subImport)
419{
420 // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
421 std::vector<RCP<Epetra_MultiVector> >::const_iterator vecItr;
422 std::vector<RCP<Epetra_Import> >::const_iterator impItr;
423
424 // using Importers fill the sub vectors from the mama vector
425 for(vecItr=many.begin(),impItr=subImport.begin();
426 vecItr!=many.end();++vecItr,++impItr) {
427 // for ease of access to the destination
428 RCP<Epetra_MultiVector> destVec = *vecItr;
429
430 // extract the map with global indicies from the current vector
431 const Epetra_Map & globalMap = *(Teuchos::get_extra_data<RCP<Epetra_Map> >(destVec,"globalMap"));
432
433 // build the import vector as a view on the destination
434 Epetra_MultiVector importVector(View,globalMap,destVec->Values(),destVec->Stride(),destVec->NumVectors());
435
436 // perform the import
437 importVector.Import(single,**impItr,Insert);
438 }
439}
440
441}
442} // end namespace Epetra
443} // end namespace Teko