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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines