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 ierr = 0;
00073 
00074   int (*RowPermute)(int in) = 0;
00075   int (*ColPermute)(int in) = 0;
00076 
00077   assert( Diagonal >= 0  && Diagonal <= 2 ); 
00078   assert( ReindexRowMap>=0 && ReindexRowMap<=2 );
00079   assert( ReindexColMap>=0 && ReindexColMap<=1 );
00080   assert( RangeMapType>=0 && RangeMapType<=3 );
00081   assert( DomainMapType>=0 && DomainMapType<=3 );
00082 
00083   Epetra_Map DomainMap = In.DomainMap();
00084   Epetra_Map RangeMap = In.RangeMap();
00085   Epetra_Map ColMap = In.ColMap();
00086   Epetra_Map RowMap = In.RowMap();
00087   int NumMyRowElements = RowMap.NumMyElements();
00088   int NumMyColElements = ColMap.NumMyElements();
00089   int NumMyRangeElements = RangeMap.NumMyElements();
00090   int NumMyDomainElements = DomainMap.NumMyElements();
00091 
00092   int NumGlobalRowElements = RowMap.NumGlobalElements();
00093   int NumGlobalColElements = ColMap.NumGlobalElements();
00094   int NumGlobalRangeElements = RangeMap.NumGlobalElements();
00095   int NumGlobalDomainElements = DomainMap.NumGlobalElements();
00096   assert( NumGlobalRangeElements == NumGlobalDomainElements ) ; 
00097 
00098   std::vector<int> MyGlobalRowElements( NumMyRowElements ) ; 
00099   std::vector<int> NumEntriesPerRow( NumMyRowElements ) ; 
00100   std::vector<int> MyPermutedGlobalRowElements( NumMyRowElements ) ; 
00101   std::vector<int> MyGlobalColElements( NumMyColElements ) ; 
00102   std::vector<int> MyPermutedGlobalColElements( NumMyColElements ) ; // Used to create the column map
00103   std::vector<int> MyPermutedGlobalColElementTable( NumMyColElements ) ; // To convert local indices to global
00104   std::vector<int> MyGlobalRangeElements( NumMyRangeElements ) ; 
00105   std::vector<int> MyPermutedGlobalRangeElements( NumMyRangeElements ) ; 
00106   std::vector<int> MyGlobalDomainElements( NumMyDomainElements ) ; 
00107   std::vector<int> MyPermutedGlobalDomainElements( NumMyDomainElements ) ; 
00108   RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
00109   ColMap.MyGlobalElements(&MyGlobalColElements[0]);
00110   RangeMap.MyGlobalElements(&MyGlobalRangeElements[0]);
00111   DomainMap.MyGlobalElements(&MyGlobalDomainElements[0]);
00112 
00113   switch( ReindexRowMap ) {
00114   case 0:
00115     RowPermute = &NoPermute ;
00116     break; 
00117   case 1:
00118     RowPermute = &SmallRowPermute ;
00119     break; 
00120   case 2:
00121     RowPermute = BigRowPermute ;
00122     break; 
00123   }
00124   switch( ReindexColMap ) {
00125   case 0:
00126     ColPermute = RowPermute ;
00127     break; 
00128   case 1:
00129     ColPermute = &SmallColPermute ;
00130     break; 
00131   }
00132 
00133   //
00134   //  Create Serial Range and Domain Maps based on the permuted indexing
00135   //
00136   int nlocal = 0;
00137   if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements;
00138   std::vector<int> AllIDs( NumGlobalRangeElements ) ; 
00139   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDs[i] = (*RowPermute)( i ) ; 
00140   Epetra_Map SerialRangeMap( -1, nlocal, &AllIDs[0], 0, In.Comm()); 
00141   std::vector<int> AllIDBs( NumGlobalRangeElements ) ; 
00142   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDBs[i] = (*ColPermute)( i ) ; 
00143   Epetra_Map SerialDomainMap( -1, nlocal, &AllIDBs[0], 0, In.Comm()); 
00144 
00145   //
00146   //  Create Bizarre Range and Domain Maps based on the permuted indexing
00147   //  These are nearly serial, having all but one element on process 0
00148   //  The goal here is to make sure that we can use Domain and Range maps 
00149   //  that are neither serial, nor distributed in the normal manner.
00150   //
00151   std::vector<int> AllIDCs( NumGlobalRangeElements ) ; 
00152   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDCs[i] = (*ColPermute)( i ) ; 
00153   if ( In.Comm().NumProc() > 1 ) { 
00154     if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
00155     if (In.Comm().MyPID()==1) {
00156       nlocal = 1;
00157       AllIDCs[0] = (*ColPermute)( NumGlobalRangeElements - 1 );
00158     }
00159   } 
00160   int iam = In.Comm().MyPID();
00161   Epetra_Map BizarreDomainMap( -1, nlocal, &AllIDCs[0], 0, In.Comm()); 
00162 
00163   std::vector<int> AllIDDs( NumGlobalRangeElements ) ; 
00164   for ( int i = 0; i < NumGlobalRangeElements ; i++ ) AllIDDs[i] = (*RowPermute)( i ) ; 
00165   if ( In.Comm().NumProc() > 1 ) { 
00166     if (In.Comm().MyPID()==0) nlocal = NumGlobalRangeElements-1;
00167     if (In.Comm().MyPID()==1) {
00168       nlocal = 1;
00169       AllIDDs[0] = (*RowPermute)( NumGlobalRangeElements -1 ) ;
00170     }
00171   } 
00172   Epetra_Map BizarreRangeMap( -1, nlocal, &AllIDDs[0], 0, In.Comm()); 
00173 
00174 
00175   //
00176   //  Compute the column map 
00177   //
00178   //  If Diagonal==1, remove the column corresponding to the last row owned 
00179   //  by process 0.  Removing this column from a tridiagonal matrix, leaves
00180   //  a disconnected, but non-singular matrix.  
00181   //
00182   int NumMyColElementsOut = 0 ; 
00183   int NumGlobalColElementsOut ; 
00184   if ( Diagonal == 1 ) 
00185     NumGlobalColElementsOut = NumGlobalColElements-1; 
00186   else
00187     NumGlobalColElementsOut = NumGlobalColElements; 
00188   if ( Diagonal == 1 && iam==0 ) { 
00189     for ( int i=0; i < NumMyColElements  ; i++ ) {
00190       if ( MyGlobalColElements[i] != MyGlobalRowElements[NumMyRowElements-1] ) {
00191   MyPermutedGlobalColElements[NumMyColElementsOut++] = 
00192     (*ColPermute)( MyGlobalColElements[i] ) ; 
00193       }
00194     }
00195     assert( NumMyColElementsOut == NumMyColElements-1 );
00196   } else {
00197     for ( int i=0; i < NumMyColElements  ; i++ )  
00198       MyPermutedGlobalColElements[i] = 
00199   (*ColPermute)( MyGlobalColElements[i] ) ; 
00200     NumMyColElementsOut = NumMyColElements ; 
00201     if ( Diagonal == 2 ) {
00202       //  For each row, make sure that the column map has this row in it, 
00203       //    if it doesn't, add it to the column map.  
00204       //  Note:  MyPermutedGlobalColElements == MyGlobalColElements when 
00205       //  Diagonal==2 because  ( Diagonal == 2 ) implies:
00206       //     ReindexRowMap==0 && ReindexColMap == 0  - see assert above
00207       for ( int i=0; i < NumMyRowElements  ; i++ ) {
00208   bool MissingDiagonal = true; 
00209   for ( int j=0; j < NumMyColElements; j++ ) { 
00210     if ( MyGlobalRowElements[i] == MyGlobalColElements[j] ) {
00211       MissingDiagonal = false; 
00212     }
00213   }
00214   if ( MissingDiagonal ) {
00215     MyPermutedGlobalColElements.resize(NumMyColElements+1);
00216     MyPermutedGlobalColElements[NumMyColElementsOut] = MyGlobalRowElements[i];
00217     NumMyColElementsOut++;
00218   }
00219       }
00220       In.Comm().SumAll(&NumMyColElementsOut,&NumGlobalColElementsOut,1); 
00221     }
00222   }
00223 
00224   //
00225   //  These tables are used both as the permutation tables and to create the maps.
00226   //
00227   for ( int i=0; i < NumMyColElements  ; i++ ) 
00228     MyPermutedGlobalColElementTable[i] = 
00229       (*ColPermute)( MyGlobalColElements[i] ) ; 
00230   for ( int i=0; i < NumMyRowElements  ; i++ ) 
00231     MyPermutedGlobalRowElements[i] = 
00232       (*RowPermute)( MyGlobalRowElements[i] ) ; 
00233   for ( int i=0; i < NumMyRangeElements  ; i++ ) 
00234     MyPermutedGlobalRangeElements[i] = 
00235       (*RowPermute)( MyGlobalRangeElements[i] ) ; 
00236   for ( int i=0; i < NumMyDomainElements  ; i++ ) 
00237     MyPermutedGlobalDomainElements[i] = 
00238       (*ColPermute)( MyGlobalDomainElements[i] ) ; 
00239 
00240   RCP<Epetra_Map> PermutedRowMap = 
00241     rcp( new Epetra_Map( NumGlobalRowElements, NumMyRowElements, 
00242        &MyPermutedGlobalRowElements[0], 0, In.Comm() ) ); 
00243                   
00244   RCP<Epetra_Map> PermutedColMap = 
00245     rcp( new Epetra_Map( NumGlobalColElementsOut, NumMyColElementsOut, 
00246        &MyPermutedGlobalColElements[0], 0, In.Comm() ) ); 
00247                   
00248   RCP<Epetra_Map> PermutedRangeMap = 
00249     rcp( new Epetra_Map( NumGlobalRangeElements, NumMyRangeElements, 
00250        &MyPermutedGlobalRangeElements[0], 0, In.Comm() ) ); 
00251                   
00252   RCP<Epetra_Map> PermutedDomainMap = 
00253     rcp( new Epetra_Map( NumGlobalDomainElements, NumMyDomainElements, 
00254        &MyPermutedGlobalDomainElements[0], 0, In.Comm() ) ); 
00255                   
00256   //
00257   //  These vectors are filled and then passed to InsertGlobalValues 
00258   //
00259   std::vector<int> ThisRowIndices( In.MaxNumEntries() );
00260   std::vector<double> ThisRowValues( In.MaxNumEntries() );
00261   std::vector<int> PermutedGlobalColIndices( In.MaxNumEntries() );
00262 
00263   //std::cout << __FILE__ << "::" <<__LINE__ << std::endl ; 
00264   RCP<Epetra_CrsMatrix> Out = 
00265     rcp( new Epetra_CrsMatrix( Copy, *PermutedRowMap, *PermutedColMap, 0 ) );
00266 
00267   for (int i=0; i<NumMyRowElements; i++)
00268     {
00269 
00270       int NumIndicesThisRow = 0;
00271       ierr = In.ExtractMyRowCopy( i, 
00272            In.MaxNumEntries(),
00273            NumIndicesThisRow,
00274            &ThisRowValues[0],
00275            &ThisRowIndices[0] );
00276       assert( ierr == 0 ) ;
00277       for (int j = 0 ; j < NumIndicesThisRow ; j++ )
00278   {
00279     PermutedGlobalColIndices[j] = MyPermutedGlobalColElementTable[ ThisRowIndices[j] ]  ;
00280   }
00281       bool MissingDiagonal = false; 
00282       if ( Diagonal==2 ) { 
00283   //
00284   assert( MyGlobalRowElements[i] == MyPermutedGlobalRowElements[i] );
00285   MissingDiagonal = true; 
00286   for( int j =0 ; j < NumIndicesThisRow ; j++ ) {
00287     if ( PermutedGlobalColIndices[j] == MyPermutedGlobalRowElements[i] ) {
00288       MissingDiagonal = false ; 
00289     }
00290   }
00291 #if 0
00292   std::cout  << __FILE__ << "::" << __LINE__ 
00293         << " i = " << i 
00294         << " MyPermutedGlobalRowElements[i]  = " << MyPermutedGlobalRowElements[i] 
00295         <<   " MissingDiagonal = " << MissingDiagonal << std::endl ; 
00296 #endif
00297 
00298       }
00299       if ( MissingDiagonal ) { 
00300   ThisRowValues.resize(NumIndicesThisRow+1) ; 
00301   ThisRowValues[NumIndicesThisRow] = 0.0;
00302   PermutedGlobalColIndices.resize(NumIndicesThisRow+1);
00303   PermutedGlobalColIndices[NumIndicesThisRow] = MyPermutedGlobalRowElements[i] ;
00304   
00305 #if 0
00306   std::cout  << __FILE__ << "::" << __LINE__ 
00307         << " i = " << i 
00308         << "NumIndicesThisRow = " << NumIndicesThisRow 
00309         << "ThisRowValues[NumIndicesThisRow = " << ThisRowValues[NumIndicesThisRow] 
00310         << " PermutedGlobalColIndices[NumIndcesThisRow] = " << PermutedGlobalColIndices[NumIndicesThisRow] 
00311         << std::endl ; 
00312 #endif
00313 
00314   NumIndicesThisRow++  ;
00315 
00316       } 
00317       ierr = Out->InsertGlobalValues( MyPermutedGlobalRowElements[i], 
00318                NumIndicesThisRow,
00319                &ThisRowValues[0],
00320                &PermutedGlobalColIndices[0] );
00321       assert( ierr >= 0 ); 
00322     }
00323 
00324   //
00325 
00326   Epetra_LocalMap ReplicatedMap( NumGlobalRangeElements, 0, In.Comm() );
00327 
00328   RCP<Epetra_Map> OutRangeMap ;
00329   RCP<Epetra_Map> OutDomainMap ;
00330   
00331   switch( RangeMapType ) {
00332   case 0:
00333     OutRangeMap = PermutedRangeMap ;
00334     break;
00335   case 1:
00336     OutRangeMap = rcp(&SerialRangeMap, false); 
00337     break;
00338   case 2:
00339     OutRangeMap = rcp(&BizarreRangeMap, false); 
00340     break;
00341   case 3:
00342     OutRangeMap = rcp(&ReplicatedMap, false); 
00343     break;
00344   }
00345   //  switch( DomainMapType ) {
00346   switch( DomainMapType ) {
00347   case 0:
00348     OutDomainMap = PermutedDomainMap ;
00349     break;
00350   case 1:
00351     OutDomainMap = rcp(&SerialDomainMap, false); 
00352     break;
00353   case 2:
00354     OutDomainMap = rcp(&BizarreDomainMap, false); 
00355     break;
00356   case 3:
00357     OutDomainMap = rcp(&ReplicatedMap, false); 
00358     break;
00359   }
00360 #if 0
00361   ierr = Out->FillComplete( *PermutedDomainMap, *PermutedRangeMap );
00362   assert( ierr == 0 );
00363 #else
00364   ierr = Out->FillComplete( *OutDomainMap, *OutRangeMap );
00365   assert( ierr == 0 );
00366 #endif
00367 
00368 #if 0
00369   std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
00370   Out->Print( std::cout ) ; 
00371 #endif
00372 
00373   return Out;
00374 }
00375 
00376 
00377 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines