Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
12 using namespace Teuchos;
13 
14 int SmallRowPermute( int in ) { return in + 2 ; }
15 int BigRowPermute( int in ) { return 5*in + 1 ; }
16 int NoPermute( int in ) { return in ; }
17 int 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 //
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 
72  int ierr = 0;
73 
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  ierr = In.ExtractMyRowCopy( i,
272  In.MaxNumEntries(),
273  NumIndicesThisRow,
274  &ThisRowValues[0],
275  &ThisRowIndices[0] );
276  assert( ierr == 0 ) ;
277  for (int j = 0 ; j < NumIndicesThisRow ; j++ )
278  {
279  PermutedGlobalColIndices[j] = MyPermutedGlobalColElementTable[ ThisRowIndices[j] ] ;
280  }
281  bool MissingDiagonal = false;
282  if ( Diagonal==2 ) {
283  //
284  assert( MyGlobalRowElements[i] == MyPermutedGlobalRowElements[i] );
285  MissingDiagonal = true;
286  for( int j =0 ; j < NumIndicesThisRow ; j++ ) {
287  if ( PermutedGlobalColIndices[j] == MyPermutedGlobalRowElements[i] ) {
288  MissingDiagonal = false ;
289  }
290  }
291 #if 0
292  std::cout << __FILE__ << "::" << __LINE__
293  << " i = " << i
294  << " MyPermutedGlobalRowElements[i] = " << MyPermutedGlobalRowElements[i]
295  << " MissingDiagonal = " << MissingDiagonal << std::endl ;
296 #endif
297 
298  }
299  if ( MissingDiagonal ) {
300  ThisRowValues.resize(NumIndicesThisRow+1) ;
301  ThisRowValues[NumIndicesThisRow] = 0.0;
302  PermutedGlobalColIndices.resize(NumIndicesThisRow+1);
303  PermutedGlobalColIndices[NumIndicesThisRow] = MyPermutedGlobalRowElements[i] ;
304 
305 #if 0
306  std::cout << __FILE__ << "::" << __LINE__
307  << " i = " << i
308  << "NumIndicesThisRow = " << NumIndicesThisRow
309  << "ThisRowValues[NumIndicesThisRow = " << ThisRowValues[NumIndicesThisRow]
310  << " PermutedGlobalColIndices[NumIndcesThisRow] = " << PermutedGlobalColIndices[NumIndicesThisRow]
311  << std::endl ;
312 #endif
313 
314  NumIndicesThisRow++ ;
315 
316  }
317  ierr = Out->InsertGlobalValues( MyPermutedGlobalRowElements[i],
318  NumIndicesThisRow,
319  &ThisRowValues[0],
320  &PermutedGlobalColIndices[0] );
321  assert( ierr >= 0 );
322  }
323 
324  //
325 
326  Epetra_LocalMap ReplicatedMap( NumGlobalRangeElements, 0, In.Comm() );
327 
328  RCP<Epetra_Map> OutRangeMap ;
329  RCP<Epetra_Map> OutDomainMap ;
330 
331  switch( RangeMapType ) {
332  case 0:
333  OutRangeMap = PermutedRangeMap ;
334  break;
335  case 1:
336  OutRangeMap = rcp(&SerialRangeMap, false);
337  break;
338  case 2:
339  OutRangeMap = rcp(&BizarreRangeMap, false);
340  break;
341  case 3:
342  OutRangeMap = rcp(&ReplicatedMap, false);
343  break;
344  }
345  // switch( DomainMapType ) {
346  switch( DomainMapType ) {
347  case 0:
348  OutDomainMap = PermutedDomainMap ;
349  break;
350  case 1:
351  OutDomainMap = rcp(&SerialDomainMap, false);
352  break;
353  case 2:
354  OutDomainMap = rcp(&BizarreDomainMap, false);
355  break;
356  case 3:
357  OutDomainMap = rcp(&ReplicatedMap, false);
358  break;
359  }
360 #if 0
361  ierr = Out->FillComplete( *PermutedDomainMap, *PermutedRangeMap );
362  assert( ierr == 0 );
363 #else
364  ierr = Out->FillComplete( *OutDomainMap, *OutRangeMap );
365  assert( ierr == 0 );
366 #endif
367 
368 #if 0
369  std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
370  Out->Print( std::cout ) ;
371 #endif
372 
373  return Out;
374 }
375 
376 
377 
int NumGlobalElements() const
const Epetra_Map & RangeMap() const
int MyGlobalElements(int *MyGlobalElementList) const
const Epetra_Map & ColMap() const
virtual int MyPID() const =0
int SmallRowPermute(int in)
const Epetra_Map & RowMap() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int MaxNumEntries() const
int NoPermute(int in)
int BigRowPermute(int in)
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
virtual int NumProc() const =0
const Epetra_Map & DomainMap() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
const Epetra_Comm & Comm() const
int SmallColPermute(int in)