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
```