Amesos Package Browser (Single Doxygen Collection) Development
NewMatNewMap.cpp
Go to the documentation of this file.
```00001
00002 #include <assert.h>
00003
00004 #include "Amesos_ConfigDefs.h"
00005 #include "Epetra_Map.h"
00006 #include "Epetra_LocalMap.h"
00007 #include "Epetra_CrsMatrix.h"
00008 #include "Epetra_Comm.h"
00009 #include "Epetra_Vector.h"
00010 #include "Teuchos_RCP.hpp"
00011
00012 using namespace Teuchos;
00013
00014 int SmallRowPermute( int in ) { return in + 2 ; }
00015 int BigRowPermute( int in ) { return 5*in + 1 ; }
00016 int NoPermute( int in ) { return in ; }
00017 int SmallColPermute( int in ) { return in + 4 ; }
00018
00019 //
00020 //  Diagonal:  0=no change, 1=eliminate entry
00021 //             from the map for the largest row element in process 0
00022 //             2=add diagonal entries to the matrix, with a zero value
00023 //             (assume row map contains all diagonal entries).
00024 //
00025 //  ReindexRowMap:
00026 //    0=no change, 1= add 2 (still contiguous), 2=non-contiguous
00027 //
00028 //  ReindexColMap
00029 //    0=same as RowMap, 1=add 4 - Different From RowMap, but contiguous)
00030 //
00031 //  RangeMap:
00032 //    0=no change, 1=serial map, 2=bizarre distribution, 3=replicated map
00033 //
00034 //  DomainMap:
00035 //    0=no change, 1=serial map, 2=bizarre distribution, 3=replicated map
00036 //
00037 RCP<Epetra_CrsMatrix> NewMatNewMap(Epetra_CrsMatrix& In,
00038              int Diagonal,
00039              int ReindexRowMap,
00040              int ReindexColMap,
00041              int RangeMapType,
00042              int DomainMapType
00043              )
00044 {
00045
00046   //
00047   //  If we are making no change, return the original matrix (which has a linear map)
00048   //
00049 #if 0
00050   std::cout << __FILE__ << "::" << __LINE__ << " "
00051        << Diagonal << " "
00052        << ReindexRowMap << " "
00053        << ReindexColMap << " "
00054        << RangeMapType << " "
00055        << DomainMapType << " " << std::endl ;
00056 #endif
00057
00058   if ( Diagonal + ReindexRowMap + ReindexColMap + RangeMapType + DomainMapType == 0 ) {
00059     RCP<Epetra_CrsMatrix> ReturnOrig = rcp( &In, false );
00060     return ReturnOrig ;
00061   }
00062
00063   //
00064   //  Diagonal==2 is used for a different purpose -
00065   //    Making sure that the diagonal of the matrix is non-empty.
00066   //  Note:  The diagonal must exist in In.RowMap().
00067   //
00068   if ( Diagonal == 2 ) {
00069     assert( ReindexRowMap==0 && ReindexColMap == 0 ) ;
00070   }
00071
00072   int (*RowPermute)(int in) = 0;
00073   int (*ColPermute)(int in) = 0;
00074
00075   assert( Diagonal >= 0  && Diagonal <= 2 );
00076   assert( ReindexRowMap>=0 && ReindexRowMap<=2 );
00077   assert( ReindexColMap>=0 && ReindexColMap<=1 );
00078   assert( RangeMapType>=0 && RangeMapType<=3 );
00079   assert( DomainMapType>=0 && DomainMapType<=3 );
00080
00081   Epetra_Map DomainMap = In.DomainMap();
00082   Epetra_Map RangeMap = In.RangeMap();
00083   Epetra_Map ColMap = In.ColMap();
00084   Epetra_Map RowMap = In.RowMap();
00085   int NumMyRowElements = RowMap.NumMyElements();
00086   int NumMyColElements = ColMap.NumMyElements();
00087   int NumMyRangeElements = RangeMap.NumMyElements();
00088   int NumMyDomainElements = DomainMap.NumMyElements();
00089
00090   int NumGlobalRowElements = RowMap.NumGlobalElements();
00091   int NumGlobalColElements = ColMap.NumGlobalElements();
00092   int NumGlobalRangeElements = RangeMap.NumGlobalElements();
00093   int NumGlobalDomainElements = DomainMap.NumGlobalElements();
00094   assert( NumGlobalRangeElements == NumGlobalDomainElements ) ;
00095
00096   std::vector<int> MyGlobalRowElements( NumMyRowElements ) ;
00097   std::vector<int> NumEntriesPerRow( NumMyRowElements ) ;
00098   std::vector<int> MyPermutedGlobalRowElements( NumMyRowElements ) ;
00099   std::vector<int> MyGlobalColElements( NumMyColElements ) ;
00100   std::vector<int> MyPermutedGlobalColElements( NumMyColElements ) ; // Used to create the column map
00101   std::vector<int> MyPermutedGlobalColElementTable( NumMyColElements ) ; // To convert local indices to global
00102   std::vector<int> MyGlobalRangeElements( NumMyRangeElements ) ;
00103   std::vector<int> MyPermutedGlobalRangeElements( NumMyRangeElements ) ;
00104   std::vector<int> MyGlobalDomainElements( NumMyDomainElements ) ;
00105   std::vector<int> MyPermutedGlobalDomainElements( NumMyDomainElements ) ;
00106   RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
00107   ColMap.MyGlobalElements(&MyGlobalColElements[0]);
00108   RangeMap.MyGlobalElements(&MyGlobalRangeElements[0]);
00109   DomainMap.MyGlobalElements(&MyGlobalDomainElements[0]);
00110
00111   switch( ReindexRowMap ) {
00112   case 0:
00113     RowPermute = &NoPermute ;
00114     break;
00115   case 1:
00116     RowPermute = &SmallRowPermute ;
00117     break;
00118   case 2:
00119     RowPermute = BigRowPermute ;
00120     break;
00121   }
00122   switch( ReindexColMap ) {
00123   case 0:
00124     ColPermute = RowPermute ;
00125     break;
00126   case 1:
00127     ColPermute = &SmallColPermute ;
00128     break;
00129   }
00130
00131   //
00132   //  Create Serial Range and Domain Maps based on the permuted indexing
00133   //
00134   int nlocal = 0;
00135   if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements;
00136   std::vector<int> AllIDs( NumGlobalRangeElements ) ;
00137   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDs[i] = (*RowPermute)( i ) ;
00138   Epetra_Map SerialRangeMap( -1, nlocal, &AllIDs[0], 0, In.Comm());
00139   std::vector<int> AllIDBs( NumGlobalRangeElements ) ;
00140   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDBs[i] = (*ColPermute)( i ) ;
00141   Epetra_Map SerialDomainMap( -1, nlocal, &AllIDBs[0], 0, In.Comm());
00142
00143   //
00144   //  Create Bizarre Range and Domain Maps based on the permuted indexing
00145   //  These are nearly serial, having all but one element on process 0
00146   //  The goal here is to make sure that we can use Domain and Range maps
00147   //  that are neither serial, nor distributed in the normal manner.
00148   //
00149   std::vector<int> AllIDCs( NumGlobalRangeElements ) ;
00150   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDCs[i] = (*ColPermute)( i ) ;
00151   if ( In.Comm().NumProc() > 1 ) {
00152     if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
00153     if (In.Comm().MyPID()==1) {
00154       nlocal = 1;
00155       AllIDCs[0] = (*ColPermute)( NumGlobalRangeElements - 1 );
00156     }
00157   }
00158   int iam = In.Comm().MyPID();
00159   Epetra_Map BizarreDomainMap( -1, nlocal, &AllIDCs[0], 0, In.Comm());
00160
00161   std::vector<int> AllIDDs( NumGlobalRangeElements ) ;
00162   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDDs[i] = (*RowPermute)( i ) ;
00163   if ( In.Comm().NumProc() > 1 ) {
00164     if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
00165     if (In.Comm().MyPID()==1) {
00166       nlocal = 1;
00167       AllIDDs[0] = (*RowPermute)( NumGlobalRangeElements -1 ) ;
00168     }
00169   }
00170   Epetra_Map BizarreRangeMap( -1, nlocal, &AllIDDs[0], 0, In.Comm());
00171
00172
00173   //
00174   //  Compute the column map
00175   //
00176   //  If Diagonal==1, remove the column corresponding to the last row owned
00177   //  by process 0.  Removing this column from a tridiagonal matrix, leaves
00178   //  a disconnected, but non-singular matrix.
00179   //
00180   int NumMyColElementsOut = 0 ;
00181   int NumGlobalColElementsOut ;
00182   if ( Diagonal == 1 )
00183     NumGlobalColElementsOut = NumGlobalColElements-1;
00184   else
00185     NumGlobalColElementsOut = NumGlobalColElements;
00186   if ( Diagonal == 1 && iam==0 ) {
00187     for ( int i=0; i < NumMyColElements  ; i++ ) {
00188       if ( MyGlobalColElements[i] != MyGlobalRowElements[NumMyRowElements-1] ) {
00189   MyPermutedGlobalColElements[NumMyColElementsOut++] =
00190     (*ColPermute)( MyGlobalColElements[i] ) ;
00191       }
00192     }
00193     assert( NumMyColElementsOut == NumMyColElements-1 );
00194   } else {
00195     for ( int i=0; i < NumMyColElements  ; i++ )
00196       MyPermutedGlobalColElements[i] =
00197   (*ColPermute)( MyGlobalColElements[i] ) ;
00198     NumMyColElementsOut = NumMyColElements ;
00199     if ( Diagonal == 2 ) {
00200       //  For each row, make sure that the column map has this row in it,
00201       //    if it doesn't, add it to the column map.
00202       //  Note:  MyPermutedGlobalColElements == MyGlobalColElements when
00203       //  Diagonal==2 because  ( Diagonal == 2 ) implies:
00204       //     ReindexRowMap==0 && ReindexColMap == 0  - see assert above
00205       for ( int i=0; i < NumMyRowElements  ; i++ ) {
00206   bool MissingDiagonal = true;
00207   for ( int j=0; j < NumMyColElements; j++ ) {
00208     if ( MyGlobalRowElements[i] == MyGlobalColElements[j] ) {
00209       MissingDiagonal = false;
00210     }
00211   }
00212   if ( MissingDiagonal ) {
00213     MyPermutedGlobalColElements.resize(NumMyColElements+1);
00214     MyPermutedGlobalColElements[NumMyColElementsOut] = MyGlobalRowElements[i];
00215     NumMyColElementsOut++;
00216   }
00217       }
00218       In.Comm().SumAll(&NumMyColElementsOut,&NumGlobalColElementsOut,1);
00219     }
00220   }
00221
00222   //
00223   //  These tables are used both as the permutation tables and to create the maps.
00224   //
00225   for ( int i=0; i < NumMyColElements  ; i++ )
00226     MyPermutedGlobalColElementTable[i] =
00227       (*ColPermute)( MyGlobalColElements[i] ) ;
00228   for ( int i=0; i < NumMyRowElements  ; i++ )
00229     MyPermutedGlobalRowElements[i] =
00230       (*RowPermute)( MyGlobalRowElements[i] ) ;
00231   for ( int i=0; i < NumMyRangeElements  ; i++ )
00232     MyPermutedGlobalRangeElements[i] =
00233       (*RowPermute)( MyGlobalRangeElements[i] ) ;
00234   for ( int i=0; i < NumMyDomainElements  ; i++ )
00235     MyPermutedGlobalDomainElements[i] =
00236       (*ColPermute)( MyGlobalDomainElements[i] ) ;
00237
00238   RCP<Epetra_Map> PermutedRowMap =
00239     rcp( new Epetra_Map( NumGlobalRowElements, NumMyRowElements,
00240        &MyPermutedGlobalRowElements[0], 0, In.Comm() ) );
00241
00242   RCP<Epetra_Map> PermutedColMap =
00243     rcp( new Epetra_Map( NumGlobalColElementsOut, NumMyColElementsOut,
00244        &MyPermutedGlobalColElements[0], 0, In.Comm() ) );
00245
00246   RCP<Epetra_Map> PermutedRangeMap =
00247     rcp( new Epetra_Map( NumGlobalRangeElements, NumMyRangeElements,
00248        &MyPermutedGlobalRangeElements[0], 0, In.Comm() ) );
00249
00250   RCP<Epetra_Map> PermutedDomainMap =
00251     rcp( new Epetra_Map( NumGlobalDomainElements, NumMyDomainElements,
00252        &MyPermutedGlobalDomainElements[0], 0, In.Comm() ) );
00253
00254   //
00255   //  These vectors are filled and then passed to InsertGlobalValues
00256   //
00257   std::vector<int> ThisRowIndices( In.MaxNumEntries() );
00258   std::vector<double> ThisRowValues( In.MaxNumEntries() );
00259   std::vector<int> PermutedGlobalColIndices( In.MaxNumEntries() );
00260
00261   //std::cout << __FILE__ << "::" <<__LINE__ << std::endl ;
00262   RCP<Epetra_CrsMatrix> Out =
00263     rcp( new Epetra_CrsMatrix( Copy, *PermutedRowMap, *PermutedColMap, 0 ) );
00264
00265   for (int i=0; i<NumMyRowElements; i++)
00266     {
00267
00268       int NumIndicesThisRow = 0;
00269       assert( In.ExtractMyRowCopy( i,
00270            In.MaxNumEntries(),
00271            NumIndicesThisRow,
00272            &ThisRowValues[0],
00273            &ThisRowIndices[0] ) == 0 ) ;
00274       for (int j = 0 ; j < NumIndicesThisRow ; j++ )
00275   {
00276     PermutedGlobalColIndices[j] = MyPermutedGlobalColElementTable[ ThisRowIndices[j] ]  ;
00277   }
00278       bool MissingDiagonal = false;
00279       if ( Diagonal==2 ) {
00280   //
00281   assert( MyGlobalRowElements[i] == MyPermutedGlobalRowElements[i] );
00282   MissingDiagonal = true;
00283   for( int j =0 ; j < NumIndicesThisRow ; j++ ) {
00284     if ( PermutedGlobalColIndices[j] == MyPermutedGlobalRowElements[i] ) {
00285       MissingDiagonal = false ;
00286     }
00287   }
00288 #if 0
00289   std::cout  << __FILE__ << "::" << __LINE__
00290         << " i = " << i
00291         << " MyPermutedGlobalRowElements[i]  = " << MyPermutedGlobalRowElements[i]
00292         <<   " MissingDiagonal = " << MissingDiagonal << std::endl ;
00293 #endif
00294
00295       }
00296       if ( MissingDiagonal ) {
00297   ThisRowValues.resize(NumIndicesThisRow+1) ;
00298   ThisRowValues[NumIndicesThisRow] = 0.0;
00299   PermutedGlobalColIndices.resize(NumIndicesThisRow+1);
00300   PermutedGlobalColIndices[NumIndicesThisRow] = MyPermutedGlobalRowElements[i] ;
00301
00302 #if 0
00303   std::cout  << __FILE__ << "::" << __LINE__
00304         << " i = " << i
00305         << "NumIndicesThisRow = " << NumIndicesThisRow
00306         << "ThisRowValues[NumIndicesThisRow = " << ThisRowValues[NumIndicesThisRow]
00307         << " PermutedGlobalColIndices[NumIndcesThisRow] = " << PermutedGlobalColIndices[NumIndicesThisRow]
00308         << std::endl ;
00309 #endif
00310
00311   NumIndicesThisRow++  ;
00312
00313       }
00314       assert( Out->InsertGlobalValues( MyPermutedGlobalRowElements[i],
00315                NumIndicesThisRow,
00316                &ThisRowValues[0],
00317                &PermutedGlobalColIndices[0] ) >= 0 );
00318     }
00319
00320   //
00321
00322   Epetra_LocalMap ReplicatedMap( NumGlobalRangeElements, 0, In.Comm() );
00323
00324   RCP<Epetra_Map> OutRangeMap ;
00325   RCP<Epetra_Map> OutDomainMap ;
00326
00327   switch( RangeMapType ) {
00328   case 0:
00329     OutRangeMap = PermutedRangeMap ;
00330     break;
00331   case 1:
00332     OutRangeMap = rcp(&SerialRangeMap, false);
00333     break;
00334   case 2:
00335     OutRangeMap = rcp(&BizarreRangeMap, false);
00336     break;
00337   case 3:
00338     OutRangeMap = rcp(&ReplicatedMap, false);
00339     break;
00340   }
00341   //  switch( DomainMapType ) {
00342   switch( DomainMapType ) {
00343   case 0:
00344     OutDomainMap = PermutedDomainMap ;
00345     break;
00346   case 1:
00347     OutDomainMap = rcp(&SerialDomainMap, false);
00348     break;
00349   case 2:
00350     OutDomainMap = rcp(&BizarreDomainMap, false);
00351     break;
00352   case 3:
00353     OutDomainMap = rcp(&ReplicatedMap, false);
00354     break;
00355   }
00356 #if 0
00357   assert(Out->FillComplete( *PermutedDomainMap, *PermutedRangeMap )==0);
00358 #else
00359   assert(Out->FillComplete( *OutDomainMap, *OutRangeMap )==0);
00360 #endif
00361
00362 #if 0
00363   std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
00364   Out->Print( std::cout ) ;
00365 #endif
00366
00367   return Out;
00368 }
00369
00370
00371
```