12 using namespace Teuchos;
50 std::cout << __FILE__ <<
"::" << __LINE__ <<
" "
52 << ReindexRowMap <<
" "
53 << ReindexColMap <<
" "
54 << RangeMapType <<
" "
55 << DomainMapType <<
" " << std::endl ;
58 if ( Diagonal + ReindexRowMap + ReindexColMap + RangeMapType + DomainMapType == 0 ) {
68 if ( Diagonal == 2 ) {
69 assert( ReindexRowMap==0 && ReindexColMap == 0 ) ;
74 int (*RowPermute)(
int in) = 0;
75 int (*ColPermute)(
int in) = 0;
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 );
96 assert( NumGlobalRangeElements == NumGlobalDomainElements ) ;
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 ) ;
103 std::vector<int> MyPermutedGlobalColElementTable( NumMyColElements ) ;
104 std::vector<int> MyGlobalRangeElements( NumMyRangeElements ) ;
105 std::vector<int> MyPermutedGlobalRangeElements( NumMyRangeElements ) ;
106 std::vector<int> MyGlobalDomainElements( NumMyDomainElements ) ;
107 std::vector<int> MyPermutedGlobalDomainElements( NumMyDomainElements ) ;
113 switch( ReindexRowMap ) {
124 switch( ReindexColMap ) {
126 ColPermute = RowPermute ;
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 ) ;
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());
151 std::vector<int> AllIDCs( NumGlobalRangeElements ) ;
152 for (
int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDCs[i] = (*ColPermute)( i ) ;
154 if (In.
Comm().
MyPID()==0) nlocal = NumGlobalRangeElements-1;
157 AllIDCs[0] = (*ColPermute)( NumGlobalRangeElements - 1 );
161 Epetra_Map BizarreDomainMap( -1, nlocal, &AllIDCs[0], 0, In.
Comm());
163 std::vector<int> AllIDDs( NumGlobalRangeElements ) ;
164 for (
int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDDs[i] = (*RowPermute)( i ) ;
166 if (In.
Comm().
MyPID()==0) nlocal = NumGlobalRangeElements-1;
169 AllIDDs[0] = (*RowPermute)( NumGlobalRangeElements -1 ) ;
172 Epetra_Map BizarreRangeMap( -1, nlocal, &AllIDDs[0], 0, In.
Comm());
182 int NumMyColElementsOut = 0 ;
183 int NumGlobalColElementsOut ;
185 NumGlobalColElementsOut = NumGlobalColElements-1;
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] ) ;
195 assert( NumMyColElementsOut == NumMyColElements-1 );
197 for (
int i=0; i < NumMyColElements ; i++ )
198 MyPermutedGlobalColElements[i] =
199 (*ColPermute)( MyGlobalColElements[i] ) ;
200 NumMyColElementsOut = NumMyColElements ;
201 if ( Diagonal == 2 ) {
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;
214 if ( MissingDiagonal ) {
215 MyPermutedGlobalColElements.resize(NumMyColElements+1);
216 MyPermutedGlobalColElements[NumMyColElementsOut] = MyGlobalRowElements[i];
217 NumMyColElementsOut++;
220 In.
Comm().
SumAll(&NumMyColElementsOut,&NumGlobalColElementsOut,1);
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] ) ;
242 &MyPermutedGlobalRowElements[0], 0, In.
Comm() ) );
245 rcp(
new Epetra_Map( NumGlobalColElementsOut, NumMyColElementsOut,
246 &MyPermutedGlobalColElements[0], 0, In.
Comm() ) );
249 rcp(
new Epetra_Map( NumGlobalRangeElements, NumMyRangeElements,
250 &MyPermutedGlobalRangeElements[0], 0, In.
Comm() ) );
253 rcp(
new Epetra_Map( NumGlobalDomainElements, NumMyDomainElements,
254 &MyPermutedGlobalDomainElements[0], 0, In.
Comm() ) );
261 std::vector<int> PermutedGlobalColIndices( In.
MaxNumEntries() );
267 for (
int i=0; i<NumMyRowElements; i++)
270 int NumIndicesThisRow = 0;
275 &ThisRowIndices[0] );
276 assert( ierr == 0 ) ;
277 for (
int j = 0 ; j < NumIndicesThisRow ; j++ )
279 PermutedGlobalColIndices[j] = MyPermutedGlobalColElementTable[ ThisRowIndices[j] ] ;
281 bool MissingDiagonal =
false;
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 ;
292 std::cout << __FILE__ <<
"::" << __LINE__
294 <<
" MyPermutedGlobalRowElements[i] = " << MyPermutedGlobalRowElements[i]
295 <<
" MissingDiagonal = " << MissingDiagonal << std::endl ;
299 if ( MissingDiagonal ) {
300 ThisRowValues.resize(NumIndicesThisRow+1) ;
301 ThisRowValues[NumIndicesThisRow] = 0.0;
302 PermutedGlobalColIndices.resize(NumIndicesThisRow+1);
303 PermutedGlobalColIndices[NumIndicesThisRow] = MyPermutedGlobalRowElements[i] ;
306 std::cout << __FILE__ <<
"::" << __LINE__
308 <<
"NumIndicesThisRow = " << NumIndicesThisRow
309 <<
"ThisRowValues[NumIndicesThisRow = " << ThisRowValues[NumIndicesThisRow]
310 <<
" PermutedGlobalColIndices[NumIndcesThisRow] = " << PermutedGlobalColIndices[NumIndicesThisRow]
314 NumIndicesThisRow++ ;
317 ierr = Out->InsertGlobalValues( MyPermutedGlobalRowElements[i],
320 &PermutedGlobalColIndices[0] );
331 switch( RangeMapType ) {
333 OutRangeMap = PermutedRangeMap ;
336 OutRangeMap =
rcp(&SerialRangeMap,
false);
339 OutRangeMap =
rcp(&BizarreRangeMap,
false);
342 OutRangeMap =
rcp(&ReplicatedMap,
false);
346 switch( DomainMapType ) {
348 OutDomainMap = PermutedDomainMap ;
351 OutDomainMap =
rcp(&SerialDomainMap,
false);
354 OutDomainMap =
rcp(&BizarreDomainMap,
false);
357 OutDomainMap =
rcp(&ReplicatedMap,
false);
361 ierr = Out->FillComplete( *PermutedDomainMap, *PermutedRangeMap );
364 ierr = Out->FillComplete( *OutDomainMap, *OutRangeMap );
369 std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
370 Out->Print( std::cout ) ;
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 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)