Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
NewMatNewMap.cpp
Go to the documentation of this file.
1
2#include <assert.h>
3
4#include "Amesos_ConfigDefs.h"
5#include "Epetra_Map.h"
6#include "Epetra_LocalMap.h"
7#include "Epetra_CrsMatrix.h"
8#include "Epetra_Comm.h"
9#include "Epetra_Vector.h"
10#include "Teuchos_RCP.hpp"
11
12using namespace Teuchos;
13
14int SmallRowPermute( int in ) { return in + 2 ; }
15int BigRowPermute( int in ) { return 5*in + 1 ; }
16int NoPermute( int in ) { return in ; }
17int SmallColPermute( int in ) { return in + 4 ; }
18
19//
20// Diagonal: 0=no change, 1=eliminate entry
21// from the map for the largest row element in process 0
22// 2=add diagonal entries to the matrix, with a zero value
23// (assume row map contains all diagonal entries).
24//
25// ReindexRowMap:
26// 0=no change, 1= add 2 (still contiguous), 2=non-contiguous
27//
28// ReindexColMap
29// 0=same as RowMap, 1=add 4 - Different From RowMap, but contiguous)
30//
31// RangeMap:
32// 0=no change, 1=serial map, 2=bizarre distribution, 3=replicated map
33//
34// DomainMap:
35// 0=no change, 1=serial map, 2=bizarre distribution, 3=replicated map
36//
37RCP<Epetra_CrsMatrix> NewMatNewMap(Epetra_CrsMatrix& In,
38 int Diagonal,
39 int ReindexRowMap,
40 int ReindexColMap,
41 int RangeMapType,
42 int DomainMapType
43 )
44{
45
46 //
47 // If we are making no change, return the original matrix (which has a linear map)
48 //
49#if 0
50 std::cout << __FILE__ << "::" << __LINE__ << " "
51 << Diagonal << " "
52 << ReindexRowMap << " "
53 << ReindexColMap << " "
54 << RangeMapType << " "
55 << DomainMapType << " " << std::endl ;
56#endif
57
58 if ( Diagonal + ReindexRowMap + ReindexColMap + RangeMapType + DomainMapType == 0 ) {
59 RCP<Epetra_CrsMatrix> ReturnOrig = rcp( &In, false );
60 return ReturnOrig ;
61 }
62
63 //
64 // Diagonal==2 is used for a different purpose -
65 // Making sure that the diagonal of the matrix is non-empty.
66 // Note: The diagonal must exist in In.RowMap().
67 //
68 if ( Diagonal == 2 ) {
69 assert( ReindexRowMap==0 && ReindexColMap == 0 ) ;
70 }
71#ifndef NDEBUG
72 int ierr = 0;
73#endif
74 int (*RowPermute)(int in) = 0;
75 int (*ColPermute)(int in) = 0;
76
77 assert( Diagonal >= 0 && Diagonal <= 2 );
78 assert( ReindexRowMap>=0 && ReindexRowMap<=2 );
79 assert( ReindexColMap>=0 && ReindexColMap<=1 );
80 assert( RangeMapType>=0 && RangeMapType<=3 );
81 assert( DomainMapType>=0 && DomainMapType<=3 );
82
83 Epetra_Map DomainMap = In.DomainMap();
84 Epetra_Map RangeMap = In.RangeMap();
85 Epetra_Map ColMap = In.ColMap();
86 Epetra_Map RowMap = In.RowMap();
87 int NumMyRowElements = RowMap.NumMyElements();
88 int NumMyColElements = ColMap.NumMyElements();
89 int NumMyRangeElements = RangeMap.NumMyElements();
90 int NumMyDomainElements = DomainMap.NumMyElements();
91
92 int NumGlobalRowElements = RowMap.NumGlobalElements();
93 int NumGlobalColElements = ColMap.NumGlobalElements();
94 int NumGlobalRangeElements = RangeMap.NumGlobalElements();
95 int NumGlobalDomainElements = DomainMap.NumGlobalElements();
96 assert( NumGlobalRangeElements == NumGlobalDomainElements ) ;
97
98 std::vector<int> MyGlobalRowElements( NumMyRowElements ) ;
99 std::vector<int> NumEntriesPerRow( NumMyRowElements ) ;
100 std::vector<int> MyPermutedGlobalRowElements( NumMyRowElements ) ;
101 std::vector<int> MyGlobalColElements( NumMyColElements ) ;
102 std::vector<int> MyPermutedGlobalColElements( NumMyColElements ) ; // Used to create the column map
103 std::vector<int> MyPermutedGlobalColElementTable( NumMyColElements ) ; // To convert local indices to global
104 std::vector<int> MyGlobalRangeElements( NumMyRangeElements ) ;
105 std::vector<int> MyPermutedGlobalRangeElements( NumMyRangeElements ) ;
106 std::vector<int> MyGlobalDomainElements( NumMyDomainElements ) ;
107 std::vector<int> MyPermutedGlobalDomainElements( NumMyDomainElements ) ;
108 RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
109 ColMap.MyGlobalElements(&MyGlobalColElements[0]);
110 RangeMap.MyGlobalElements(&MyGlobalRangeElements[0]);
111 DomainMap.MyGlobalElements(&MyGlobalDomainElements[0]);
112
113 switch( ReindexRowMap ) {
114 case 0:
115 RowPermute = &NoPermute ;
116 break;
117 case 1:
118 RowPermute = &SmallRowPermute ;
119 break;
120 case 2:
121 RowPermute = BigRowPermute ;
122 break;
123 }
124 switch( ReindexColMap ) {
125 case 0:
126 ColPermute = RowPermute ;
127 break;
128 case 1:
129 ColPermute = &SmallColPermute ;
130 break;
131 }
132
133 //
134 // Create Serial Range and Domain Maps based on the permuted indexing
135 //
136 int nlocal = 0;
137 if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements;
138 std::vector<int> AllIDs( NumGlobalRangeElements ) ;
139 for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDs[i] = (*RowPermute)( i ) ;
140 Epetra_Map SerialRangeMap( -1, nlocal, &AllIDs[0], 0, In.Comm());
141 std::vector<int> AllIDBs( NumGlobalRangeElements ) ;
142 for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDBs[i] = (*ColPermute)( i ) ;
143 Epetra_Map SerialDomainMap( -1, nlocal, &AllIDBs[0], 0, In.Comm());
144
145 //
146 // Create Bizarre Range and Domain Maps based on the permuted indexing
147 // These are nearly serial, having all but one element on process 0
148 // The goal here is to make sure that we can use Domain and Range maps
149 // that are neither serial, nor distributed in the normal manner.
150 //
151 std::vector<int> AllIDCs( NumGlobalRangeElements ) ;
152 for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDCs[i] = (*ColPermute)( i ) ;
153 if ( In.Comm().NumProc() > 1 ) {
154 if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
155 if (In.Comm().MyPID()==1) {
156 nlocal = 1;
157 AllIDCs[0] = (*ColPermute)( NumGlobalRangeElements - 1 );
158 }
159 }
160 int iam = In.Comm().MyPID();
161 Epetra_Map BizarreDomainMap( -1, nlocal, &AllIDCs[0], 0, In.Comm());
162
163 std::vector<int> AllIDDs( NumGlobalRangeElements ) ;
164 for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDDs[i] = (*RowPermute)( i ) ;
165 if ( In.Comm().NumProc() > 1 ) {
166 if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
167 if (In.Comm().MyPID()==1) {
168 nlocal = 1;
169 AllIDDs[0] = (*RowPermute)( NumGlobalRangeElements -1 ) ;
170 }
171 }
172 Epetra_Map BizarreRangeMap( -1, nlocal, &AllIDDs[0], 0, In.Comm());
173
174
175 //
176 // Compute the column map
177 //
178 // If Diagonal==1, remove the column corresponding to the last row owned
179 // by process 0. Removing this column from a tridiagonal matrix, leaves
180 // a disconnected, but non-singular matrix.
181 //
182 int NumMyColElementsOut = 0 ;
183 int NumGlobalColElementsOut ;
184 if ( Diagonal == 1 )
185 NumGlobalColElementsOut = NumGlobalColElements-1;
186 else
187 NumGlobalColElementsOut = NumGlobalColElements;
188 if ( Diagonal == 1 && iam==0 ) {
189 for ( int i=0; i < NumMyColElements ; i++ ) {
190 if ( MyGlobalColElements[i] != MyGlobalRowElements[NumMyRowElements-1] ) {
191 MyPermutedGlobalColElements[NumMyColElementsOut++] =
192 (*ColPermute)( MyGlobalColElements[i] ) ;
193 }
194 }
195 assert( NumMyColElementsOut == NumMyColElements-1 );
196 } else {
197 for ( int i=0; i < NumMyColElements ; i++ )
198 MyPermutedGlobalColElements[i] =
199 (*ColPermute)( MyGlobalColElements[i] ) ;
200 NumMyColElementsOut = NumMyColElements ;
201 if ( Diagonal == 2 ) {
202 // For each row, make sure that the column map has this row in it,
203 // if it doesn't, add it to the column map.
204 // Note: MyPermutedGlobalColElements == MyGlobalColElements when
205 // Diagonal==2 because ( Diagonal == 2 ) implies:
206 // ReindexRowMap==0 && ReindexColMap == 0 - see assert above
207 for ( int i=0; i < NumMyRowElements ; i++ ) {
208 bool MissingDiagonal = true;
209 for ( int j=0; j < NumMyColElements; j++ ) {
210 if ( MyGlobalRowElements[i] == MyGlobalColElements[j] ) {
211 MissingDiagonal = false;
212 }
213 }
214 if ( MissingDiagonal ) {
215 MyPermutedGlobalColElements.resize(NumMyColElements+1);
216 MyPermutedGlobalColElements[NumMyColElementsOut] = MyGlobalRowElements[i];
217 NumMyColElementsOut++;
218 }
219 }
220 In.Comm().SumAll(&NumMyColElementsOut,&NumGlobalColElementsOut,1);
221 }
222 }
223
224 //
225 // These tables are used both as the permutation tables and to create the maps.
226 //
227 for ( int i=0; i < NumMyColElements ; i++ )
228 MyPermutedGlobalColElementTable[i] =
229 (*ColPermute)( MyGlobalColElements[i] ) ;
230 for ( int i=0; i < NumMyRowElements ; i++ )
231 MyPermutedGlobalRowElements[i] =
232 (*RowPermute)( MyGlobalRowElements[i] ) ;
233 for ( int i=0; i < NumMyRangeElements ; i++ )
234 MyPermutedGlobalRangeElements[i] =
235 (*RowPermute)( MyGlobalRangeElements[i] ) ;
236 for ( int i=0; i < NumMyDomainElements ; i++ )
237 MyPermutedGlobalDomainElements[i] =
238 (*ColPermute)( MyGlobalDomainElements[i] ) ;
239
240 RCP<Epetra_Map> PermutedRowMap =
241 rcp( new Epetra_Map( NumGlobalRowElements, NumMyRowElements,
242 &MyPermutedGlobalRowElements[0], 0, In.Comm() ) );
243
244 RCP<Epetra_Map> PermutedColMap =
245 rcp( new Epetra_Map( NumGlobalColElementsOut, NumMyColElementsOut,
246 &MyPermutedGlobalColElements[0], 0, In.Comm() ) );
247
248 RCP<Epetra_Map> PermutedRangeMap =
249 rcp( new Epetra_Map( NumGlobalRangeElements, NumMyRangeElements,
250 &MyPermutedGlobalRangeElements[0], 0, In.Comm() ) );
251
252 RCP<Epetra_Map> PermutedDomainMap =
253 rcp( new Epetra_Map( NumGlobalDomainElements, NumMyDomainElements,
254 &MyPermutedGlobalDomainElements[0], 0, In.Comm() ) );
255
256 //
257 // These vectors are filled and then passed to InsertGlobalValues
258 //
259 std::vector<int> ThisRowIndices( In.MaxNumEntries() );
260 std::vector<double> ThisRowValues( In.MaxNumEntries() );
261 std::vector<int> PermutedGlobalColIndices( In.MaxNumEntries() );
262
263 //std::cout << __FILE__ << "::" <<__LINE__ << std::endl ;
264 RCP<Epetra_CrsMatrix> Out =
265 rcp( new Epetra_CrsMatrix( Copy, *PermutedRowMap, *PermutedColMap, 0 ) );
266
267 for (int i=0; i<NumMyRowElements; i++)
268 {
269
270 int NumIndicesThisRow = 0;
271#ifndef NDEBUG
272 ierr =
273#endif
274 In.ExtractMyRowCopy( i,
275 In.MaxNumEntries(),
276 NumIndicesThisRow,
277 &ThisRowValues[0],
278 &ThisRowIndices[0] );
279 assert( ierr == 0 ) ;
280 for (int j = 0 ; j < NumIndicesThisRow ; j++ )
281 {
282 PermutedGlobalColIndices[j] = MyPermutedGlobalColElementTable[ ThisRowIndices[j] ] ;
283 }
284 bool MissingDiagonal = false;
285 if ( Diagonal==2 ) {
286 //
287 assert( MyGlobalRowElements[i] == MyPermutedGlobalRowElements[i] );
288 MissingDiagonal = true;
289 for( int j =0 ; j < NumIndicesThisRow ; j++ ) {
290 if ( PermutedGlobalColIndices[j] == MyPermutedGlobalRowElements[i] ) {
291 MissingDiagonal = false ;
292 }
293 }
294#if 0
295 std::cout << __FILE__ << "::" << __LINE__
296 << " i = " << i
297 << " MyPermutedGlobalRowElements[i] = " << MyPermutedGlobalRowElements[i]
298 << " MissingDiagonal = " << MissingDiagonal << std::endl ;
299#endif
300
301 }
302 if ( MissingDiagonal ) {
303 ThisRowValues.resize(NumIndicesThisRow+1) ;
304 ThisRowValues[NumIndicesThisRow] = 0.0;
305 PermutedGlobalColIndices.resize(NumIndicesThisRow+1);
306 PermutedGlobalColIndices[NumIndicesThisRow] = MyPermutedGlobalRowElements[i] ;
307
308#if 0
309 std::cout << __FILE__ << "::" << __LINE__
310 << " i = " << i
311 << "NumIndicesThisRow = " << NumIndicesThisRow
312 << "ThisRowValues[NumIndicesThisRow = " << ThisRowValues[NumIndicesThisRow]
313 << " PermutedGlobalColIndices[NumIndcesThisRow] = " << PermutedGlobalColIndices[NumIndicesThisRow]
314 << std::endl ;
315#endif
316
317 NumIndicesThisRow++ ;
318
319 }
320#ifndef NDEBUG
321 ierr =
322#endif
323 Out->InsertGlobalValues( MyPermutedGlobalRowElements[i],
324 NumIndicesThisRow,
325 &ThisRowValues[0],
326 &PermutedGlobalColIndices[0] );
327 assert( ierr >= 0 );
328 }
329
330 //
331
332 Epetra_LocalMap ReplicatedMap( NumGlobalRangeElements, 0, In.Comm() );
333
334 RCP<Epetra_Map> OutRangeMap ;
335 RCP<Epetra_Map> OutDomainMap ;
336
337 switch( RangeMapType ) {
338 case 0:
339 OutRangeMap = PermutedRangeMap ;
340 break;
341 case 1:
342 OutRangeMap = rcp(&SerialRangeMap, false);
343 break;
344 case 2:
345 OutRangeMap = rcp(&BizarreRangeMap, false);
346 break;
347 case 3:
348 OutRangeMap = rcp(&ReplicatedMap, false);
349 break;
350 }
351 // switch( DomainMapType ) {
352 switch( DomainMapType ) {
353 case 0:
354 OutDomainMap = PermutedDomainMap ;
355 break;
356 case 1:
357 OutDomainMap = rcp(&SerialDomainMap, false);
358 break;
359 case 2:
360 OutDomainMap = rcp(&BizarreDomainMap, false);
361 break;
362 case 3:
363 OutDomainMap = rcp(&ReplicatedMap, false);
364 break;
365 }
366#if 0
367#ifndef NDEBUG
368 ierr =
369#endif
370 Out->FillComplete( *PermutedDomainMap, *PermutedRangeMap );
371 assert( ierr == 0 );
372#else
373#ifndef NDEBUG
374 ierr =
375#endif
376 Out->FillComplete( *OutDomainMap, *OutRangeMap );
377 assert( ierr == 0 );
378#endif
379
380#if 0
381 std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
382 Out->Print( std::cout ) ;
383#endif
384
385 return Out;
386}
387
388
389
int BigRowPermute(int in)
int NoPermute(int in)
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
int SmallRowPermute(int in)
int SmallColPermute(int in)