Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_OverlappingRowMatrix.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Ifpack_OverlappingRowMatrix.h"
00032 #include "Epetra_RowMatrix.h"
00033 #include "Epetra_Map.h"
00034 #include "Epetra_CrsMatrix.h"
00035 #include "Epetra_Comm.h"
00036 #include "Epetra_MultiVector.h"
00037 
00038 #ifdef IFPACK_NODE_AWARE_CODE
00039 #include "Epetra_IntVector.h"
00040 #include "Epetra_MpiComm.h"
00041 #include "Teuchos_Hashtable.hpp"
00042 //#include "Ifpack_HashTable.hpp"
00043 #include "Teuchos_Array.hpp"
00044 #include "EpetraExt_OperatorOut.h"
00045 #include "EpetraExt_BlockMapOut.h"
00046 extern int ML_NODE_ID;
00047 int IFPACK_NODE_ID;
00048 #endif
00049 
00050 using namespace Teuchos;
00051 
00052 // ====================================================================== 
00053 // Constructor for the case of one core per subdomain
00054 Ifpack_OverlappingRowMatrix::
00055 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00056                             int OverlapLevel_in)  :
00057   Matrix_(Matrix_in),
00058   OverlapLevel_(OverlapLevel_in)
00059 {
00060   // should not be here if no overlap
00061   if (OverlapLevel_in == 0)
00062     IFPACK_CHK_ERRV(-1);
00063 
00064   // nothing to do as well with one process
00065   if (Comm().NumProc() == 1)
00066     IFPACK_CHK_ERRV(-1);
00067   
00068   NumMyRowsA_ = A().NumMyRows();
00069 
00070   // construct the external matrix
00071   vector<int> ExtElements; 
00072 
00073   RCP<Epetra_Map> TmpMap;
00074   RCP<Epetra_CrsMatrix> TmpMatrix; 
00075   RCP<Epetra_Import> TmpImporter;
00076 
00077   // importing rows corresponding to elements that are 
00078   // in ColMap, but not in RowMap 
00079   const Epetra_Map *RowMap; 
00080   const Epetra_Map *ColMap; 
00081 
00082   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap) {
00083     if (TmpMatrix != Teuchos::null) {
00084       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00085       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00086     }
00087     else {
00088       RowMap = &(A().RowMatrixRowMap()); 
00089       ColMap = &(A().RowMatrixColMap()); 
00090     }
00091 
00092     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00093     vector<int> list(size); 
00094 
00095     int count = 0; 
00096 
00097     // define the set of rows that are in ColMap but not in RowMap 
00098     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) { 
00099       int GID = ColMap->GID(i); 
00100       if (A().RowMatrixRowMap().LID(GID) == -1) { 
00101         vector<int>::iterator pos 
00102           = find(ExtElements.begin(),ExtElements.end(),GID); 
00103         if (pos == ExtElements.end()) { 
00104           ExtElements.push_back(GID);
00105           list[count] = GID; 
00106           ++count; 
00107         } 
00108       } 
00109     } 
00110 
00111     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) ); 
00112 
00113     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
00114 
00115     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
00116 
00117     TmpMatrix->Import(A(),*TmpImporter,Insert); 
00118     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
00119 
00120   }
00121 
00122   // build the map containing all the nodes (original
00123   // matrix + extended matrix)
00124   vector<int> list(NumMyRowsA_ + ExtElements.size());
00125   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
00126     list[i] = A().RowMatrixRowMap().GID(i);
00127   for (int i = 0 ; i < (int)ExtElements.size() ; ++i)
00128     list[i + NumMyRowsA_] = ExtElements[i];
00129 
00130   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ExtElements.size(),
00131                 &list[0], 0, Comm()) );
00132 # ifdef IFPACK_NODE_AWARE_CODE
00133   colMap_ = &*Map_;
00134 # endif
00135   // now build the map corresponding to all the external nodes
00136   // (with respect to A().RowMatrixRowMap().
00137   ExtMap_ = rcp( new Epetra_Map(-1,ExtElements.size(),
00138            &ExtElements[0],0,A().Comm()) );
00139   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*Map_,0) ); 
00140 
00141   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
00142   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
00143   ExtMatrix_->FillComplete(A().OperatorDomainMap(),*Map_);
00144 
00145   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) );
00146 
00147   // fix indices for overlapping matrix
00148   NumMyRowsB_ = B().NumMyRows();
00149   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;
00150   NumMyCols_ = NumMyRows_;
00151   
00152   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
00153   
00154   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
00155   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00156   MaxNumEntries_ = A().MaxNumEntries();
00157   
00158   if (MaxNumEntries_ < B().MaxNumEntries())
00159     MaxNumEntries_ = B().MaxNumEntries();
00160 
00161 }
00162 
00163 #ifdef IFPACK_NODE_AWARE_CODE
00164 
00165 // ====================================================================== 
00166 // Constructor for the case of two or more cores per subdomain
00167 Ifpack_OverlappingRowMatrix::
00168 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00169                             int OverlapLevel_in, int nodeID)  :
00170   Matrix_(Matrix_in),
00171   OverlapLevel_(OverlapLevel_in)
00172 {
00173 
00174   //FIXME timing
00175 #ifdef IFPACK_OVA_TIME_BUILD
00176   double t0 = MPI_Wtime();
00177   double t1, true_t0=t0;
00178 #endif
00179   //FIXME timing
00180 
00181   const int votesMax = INT_MAX;
00182 
00183   // should not be here if no overlap
00184   if (OverlapLevel_in == 0)
00185     IFPACK_CHK_ERRV(-1);
00186 
00187   // nothing to do as well with one process
00188   if (Comm().NumProc() == 1)
00189     IFPACK_CHK_ERRV(-1);
00190 
00191   // nodeID is the node (aka socket) number, and so is system dependent
00192   // nodeComm is the communicator for all the processes on a particular node
00193   // these processes will have the same nodeID.
00194 # ifdef HAVE_MPI
00195   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &Comm() );
00196   assert(pComm != NULL);
00197   MPI_Comm nodeMPIComm;
00198   MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm);
00199   Epetra_MpiComm *nodeComm = new Epetra_MpiComm(nodeMPIComm);
00200 # else
00201   Epetra_SerialComm *nodeComm =  dynamic_cast<Epetra_MpiComm*>( &(Matrix_in->RowMatrixRowMap().Comm()) );
00202 # endif
00203   
00204   NumMyRowsA_ = A().NumMyRows();
00205 
00206   RCP<Epetra_Map> TmpMap;
00207   RCP<Epetra_CrsMatrix> TmpMatrix; 
00208   RCP<Epetra_Import> TmpImporter;
00209 
00210   // importing rows corresponding to elements that are 
00211   // in ColMap, but not in RowMap 
00212   const Epetra_Map *RowMap; 
00213   const Epetra_Map *ColMap; 
00214   const Epetra_Map *DomainMap;
00215 
00216   // TODO Count #connections from nodes I own to each ghost node
00217 
00218 
00219   //FIXME timing
00220 #ifdef IFPACK_OVA_TIME_BUILD
00221   t1 = MPI_Wtime(); 
00222   fprintf(stderr,"[node %d]: time for initialization %2.3e\n", nodeID, t1-t0);
00223   t0=t1;
00224 #endif
00225   //FIXME timing
00226 
00227 #ifdef IFPACK_OVA_TIME_BUILD
00228   t1 = MPI_Wtime();
00229   fprintf(stderr,"[node %d]: overlap hash table allocs %2.3e\n", nodeID, t1-t0);
00230   t0=t1;
00231 #endif
00232 
00233   Teuchos::Hashtable<int,int> colMapTable(3 * A().RowMatrixColMap().NumMyElements() );
00234 
00235   // ghostTable holds off-node GIDs that are connected to on-node rows and can potentially be this PID's overlap
00236   // TODO hopefully 3 times the # column entries is a reasonable table size
00237   Teuchos::Hashtable<int,int> ghostTable(3 * A().RowMatrixColMap().NumMyElements() );
00238 
00239   /* ** ************************************************************************** ** */
00240   /* ** ********************** start of main overlap loop ************************ ** */
00241   /* ** ************************************************************************** ** */
00242   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00243   {
00244     // nbTable holds node buddy GIDs that are connected to current PID's rows, i.e., GIDs that should be in the overlapped
00245     // matrix's column map
00246 
00247     if (TmpMatrix != Teuchos::null) {
00248       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00249       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00250       DomainMap = &(TmpMatrix->OperatorDomainMap());
00251     }
00252     else {
00253       RowMap = &(A().RowMatrixRowMap()); 
00254       ColMap = &(A().RowMatrixColMap()); 
00255       DomainMap = &(A().OperatorDomainMap());
00256     }
00257 
00258 #ifdef IFPACK_OVA_TIME_BUILD
00259     t1 = MPI_Wtime();
00260     fprintf(stderr,"[node %d]: overlap pointer pulls %2.3e\n", nodeID, t1-t0);
00261     t0=t1;
00262 #endif
00263     
00264     // For each column ID, determine the owning node (as opposed to core)
00265     // ID of the corresponding row.
00266     Epetra_IntVector colIdList( *ColMap );
00267     Epetra_IntVector rowIdList(*DomainMap);
00268     rowIdList.PutValue(nodeID);  
00269     Teuchos::RCP<Epetra_Import> nodeIdImporter = rcp(new Epetra_Import( *ColMap, *DomainMap ));
00270 
00271 #ifdef IFPACK_OVA_TIME_BUILD
00272     t1 = MPI_Wtime();
00273     fprintf(stderr,"[node %d]: overlap intvector/importer allocs %2.3e\n", nodeID, t1-t0);
00274     t0=t1;
00275 #endif
00276     
00277     colIdList.Import(rowIdList,*nodeIdImporter,Insert);
00278     
00279     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00280     int count = 0; 
00281 
00282 #ifdef IFPACK_OVA_TIME_BUILD
00283     t1 = MPI_Wtime();
00284     fprintf(stderr,"[node %d]: overlap col/row ID imports %2.3e\n", nodeID, t1-t0);
00285     t0=t1;
00286 #endif
00287 
00288     
00289     // define the set of off-node rows that are in ColMap but not in RowMap
00290     // This naturally does not include off-core rows that are on the same node as me, i.e., node buddy rows.
00291     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
00292       int GID = ColMap->GID(i); 
00293       int myvotes=0;
00294       if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
00295       if ( colIdList[i] != nodeID && myvotes < votesMax) // don't include anybody found in a previous round
00296       {
00297         int votes;
00298         if (ghostTable.containsKey(GID)) {
00299           votes = ghostTable.get(GID);
00300           votes++;
00301           ghostTable.put(GID,votes);
00302         } else {
00303           ghostTable.put(GID,1);
00304         }
00305       }
00306     } //for (int i = 0 ; i < ColMap->NumMyElements() ; ++i)
00307 
00308     Teuchos::Array<int> gidsarray,votesarray;
00309     ghostTable.arrayify(gidsarray,votesarray);
00310     int *gids = gidsarray.getRawPtr();
00311     int *votes = votesarray.getRawPtr();
00312 
00313     /*
00314        This next bit of code decides which node buddy (NB) gets which ghost points.  Everyone sends their
00315        list of ghost points to pid 0 of the local subcommunicator.  Pid 0 decides who gets what:
00316 
00317           - if a ghost point is touched by only one NB, that NB gets the ghost point
00318           - if two or more NBs share a ghost point, whichever NB has the most connections to the ghost
00319             point gets it.
00320     */
00321 
00322 #   ifdef HAVE_MPI  //FIXME What if we build in serial?!  This file will likely not compile.
00323     int *cullList;
00324     int ncull;
00325     //int mypid = nodeComm->MyPID();
00326 
00327     //FIXME timing
00328 #ifdef IFPACK_OVA_TIME_BUILD
00329     t1 = MPI_Wtime();
00330     fprintf(stderr,"[node %d]: overlap before culling %2.3e\n", nodeID, t1-t0);
00331     t0=t1;
00332 #endif
00333     //FIXME timing
00334 
00335 
00336     if (nodeComm->MyPID() == 0)
00337     {
00338       // Figure out how much pid 0 is to receive
00339       MPI_Status status;
00340       int *lengths = new int[nodeComm->NumProc()+1];
00341       lengths[0] = 0;
00342       lengths[1] = ghostTable.size();
00343       for (int i=1; i<nodeComm->NumProc(); i++) {
00344         int leng;
00345         MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00346         lengths[i+1] = lengths[i] + leng;
00347       }
00348       int total = lengths[nodeComm->NumProc()];
00349 
00350       int* ghosts = new int[total];
00351       for (int i=0; i<total; i++) ghosts[i] = -9;
00352       int *round  = new int[total];
00353       int *owningpid  = new int[total];
00354 
00355       for (int i=1; i<nodeComm->NumProc(); i++) {
00356         int count = lengths[i+1] - lengths[i];
00357         MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00358         MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00359       }
00360 
00361       // put in pid 0's info
00362       for (int i=0; i<lengths[1]; i++) {
00363         ghosts[i] = gids[i];
00364         round[i] = votes[i];
00365         owningpid[i] = 0;
00366       }
00367 
00368       // put in the pid associated with each ghost
00369       for (int j=1; j<nodeComm->NumProc(); j++) {
00370         for (int i=lengths[j]; i<lengths[j+1]; i++) {
00371           owningpid[i] = j;
00372         }
00373       }
00374 
00375       // sort everything based on the ghost gids
00376       int* roundpid[2];
00377       roundpid[0] = round; roundpid[1] = owningpid;
00378       Epetra_Util epetraUtil;
00379       epetraUtil.Sort(true,total,ghosts,0,0,2,roundpid);
00380 
00381       //set up arrays that get sent back to node buddies and that tell them what ghosts to cull
00382       int *nlosers = new int[nodeComm->NumProc()];
00383       int** losers = new int*[nodeComm->NumProc()];
00384       for (int i=0; i<nodeComm->NumProc(); i++) {
00385         nlosers[i] = 0;
00386         losers[i] = new int[ lengths[i+1]-lengths[i] ];
00387       }
00388 
00389       // Walk through ghosts array and and for each sequence of duplicate ghost GIDs, choose just one NB to keep it.
00390       // The logic is pretty simple.  The ghost list is sorted, so all duplicate PIDs are together.
00391       // The list is traversed.  As duplicates are found, node pid 0 keeps track of the current "winning"
00392       // pid.  When a pid is determined to have "lost" (less votes/connections to the current GID), the
00393       // GID is added to that pid's list of GIDs to be culled.  At the end of the repeated sequence, we have
00394       // a winner, and other NBs know whether they need to delete it from their import list.
00395       int max=0;   //for duplicated ghosts, index of pid with most votes and who hence keeps the ghost.
00396                    // TODO to break ties randomly
00397 
00398       for (int i=1; i<total; i++) {
00399         if (ghosts[i] == ghosts[i-1]) {
00400           int idx = i; // pid associated with idx is current "loser"
00401           if (round[i] > round[max]) {
00402             idx = max;
00403             max=i;
00404           }
00405           int j = owningpid[idx];
00406           losers[j][nlosers[j]++] = ghosts[idx];
00407         } else {
00408           max=i;
00409         }
00410       } //for (int i=1; i<total; i++)
00411 
00412       delete [] round;
00413       delete [] ghosts;
00414       delete [] owningpid;
00415 
00416       // send the arrays of ghost GIDs to be culled back to the respective node buddies
00417       for (int i=1; i<nodeComm->NumProc(); i++) {
00418         MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
00419         MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
00420       }
00421 
00422       //FIXME Unnecessary memory allocation and copying, but makes culling code cleaner
00423       //Could we stick this info into gids and votes, since neither is being used anymore?
00424       //TODO Instead of using "losers" arrays, just use "cullList" as in the else clause
00425       ncull = nlosers[0];
00426       cullList = new int[ncull+1];
00427       for (int i=0; i<nlosers[0]; i++)
00428         cullList[i] = losers[0][i];
00429 
00430       for (int i=0; i<nodeComm->NumProc(); i++)
00431         delete [] losers[i];
00432 
00433       delete [] lengths;
00434       delete [] losers;
00435       delete [] nlosers;
00436 
00437     } else { //everyone but pid 0
00438 
00439       // send to node pid 0 all ghosts that this pid could potentially import
00440       int hashsize = ghostTable.size();
00441       MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
00442       MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
00443       MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
00444 
00445       // receive the ghost GIDs that should not be imported (subset of the list sent off just above)
00446       MPI_Status status;
00447       MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
00448       cullList = new int[ncull+1];
00449       MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
00450     }
00451 
00452 
00453     //TODO clean out hash table after each time through overlap loop   4/1/07 JJH done moved both hash tables to local scope
00454 
00455     // Remove from my row map all off-node ghosts that will be imported by a node buddy.
00456     for (int i=0; i<ncull; i++) {
00457       try{ghostTable.remove(cullList[i]);}
00458 
00459       catch(...) {
00460         printf("pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
00461                Comm().MyPID(),cullList[i]);
00462         fflush(stdout);
00463       }
00464     }//for
00465 
00466     delete [] cullList;
00467 
00468     // Save off the remaining ghost GIDs from the current overlap round.
00469     // These are off-node GIDs (rows) that I will import.
00470     gidsarray.clear(); votesarray.clear();
00471     ghostTable.arrayify(gidsarray,votesarray);
00472 
00473     vector<int> list(size); 
00474     count=0;
00475     for (int i=0; i<ghostTable.size(); i++) {
00476       // if votesarray[i] == votesmax, then this GID was found during a previous overlap round
00477       if (votesarray[i] < votesMax) {
00478         list[count] = gidsarray[i];
00479         ghostTable.put(gidsarray[i],votesMax); //set this GID's vote to the maximum so that this PID alway wins in any subsequent overlap tiebreaking.
00480         count++;
00481       }
00482     }
00483 
00484     //FIXME timing
00485 #ifdef IFPACK_OVA_TIME_BUILD
00486     t1 = MPI_Wtime();
00487     fprintf(stderr,"[node %d]: overlap duplicate removal %2.3e\n", nodeID, t1-t0);
00488     t0=t1;
00489 #endif
00490     //FIXME timing
00491 
00492 #   endif //ifdef HAVE_MPI
00493 
00494     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) );
00495 
00496     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
00497 
00498     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
00499 
00500     TmpMatrix->Import(A(),*TmpImporter,Insert); 
00501     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
00502 
00503     //FIXME timing
00504 #ifdef IFPACK_OVA_TIME_BUILD
00505     t1 = MPI_Wtime();
00506     fprintf(stderr,"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", nodeID, t1-t0);
00507     t0=t1;
00508 #endif
00509     //FIXME timing
00510 
00511 
00512     // These next two imports get the GIDs that need to go into the column map of the overlapped matrix.
00513 
00514     // For each column ID in the overlap, determine the owning node (as opposed to core)
00515     // ID of the corresponding row.  Save those column IDs whose owning node is the current one.
00516     // This should get all the imported ghost GIDs.
00517     Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
00518     ov_colIdList.PutValue(-1);
00519     Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
00520     ov_rowIdList.PutValue(nodeID);  
00521     Teuchos::RCP<Epetra_Import> ov_nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), TmpMatrix->RowMap()));
00522     ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
00523 
00524     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
00525        if (ov_colIdList[i] == nodeID) {
00526          int GID = (ov_colIdList.Map()).GID(i);
00527          colMapTable.put(GID,1);
00528        }
00529     }
00530 
00531     // Do a second import of the owning node ID from A's rowmap to TmpMat's column map.  This ensures that
00532     // all GIDs that belong to a node buddy and are in a ghost row's sparsity pattern will be in the final
00533     // overlapped matrix's column map.
00534     ov_colIdList.PutValue(-1);
00535     Epetra_IntVector ArowIdList( A().RowMatrixRowMap() );
00536     ArowIdList.PutValue(nodeID);
00537     nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), A().RowMatrixRowMap() ));
00538     ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
00539 
00540     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
00541        if (ov_colIdList[i] == nodeID) {
00542          int GID = (ov_colIdList.Map()).GID(i);
00543          colMapTable.put(GID,1);
00544        }
00545     }
00546 
00547     //FIXME timing
00548 #ifdef IFPACK_OVA_TIME_BUILD
00549     t1 = MPI_Wtime();
00550     fprintf(stderr,"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", nodeID, t1-t0);
00551     t0=t1;
00552 #endif
00553     //FIXME timing
00554 
00555   } //for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00556 
00557   /* ** ************************************************************************ ** */
00558   /* ** ********************** end of main overlap loop ************************ ** */
00559   /* ** ************************************************************************ ** */
00560 
00561   // off-node GIDs that will be in the overlap
00562   vector<int> ghostElements; 
00563 
00564   Teuchos::Array<int> gidsarray,votesarray;
00565   ghostTable.arrayify(gidsarray,votesarray);
00566   for (int i=0; i<ghostTable.size(); i++) {
00567     ghostElements.push_back(gidsarray[i]);
00568   }
00569 
00570     for (int i = 0 ; i < A().RowMatrixColMap().NumMyElements() ; ++i) {
00571       int GID = A().RowMatrixColMap().GID(i);
00572       // Remove any entries that are in A's original column map
00573       if (colMapTable.containsKey(GID)) {
00574         try{colMapTable.remove(GID);}
00575         catch(...) {
00576           printf("pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n", Comm().MyPID(),GID);
00577           fflush(stdout);
00578         }
00579       }
00580     }
00581 
00582     // GIDs that will be in the overlapped matrix's column map
00583     vector<int> colMapElements; 
00584 
00585   gidsarray.clear(); votesarray.clear();
00586   colMapTable.arrayify(gidsarray,votesarray);
00587   for (int i=0; i<colMapTable.size(); i++)
00588     colMapElements.push_back(gidsarray[i]);
00589 
00590 /*
00591    We need two maps here.  The first is the row map, which we've got by using my original rows
00592    plus whatever I've picked up in ghostElements.
00593 
00594    The second is a column map.  This map should include my rows, plus any columns that could come from node buddies.
00595    These GIDs come from the std:array colMapElements, which in turn comes from colMapTable.
00596    This map should *not* omit ghosts that have been imported by my node buddies, i.e., for any row that I own,
00597    the stencil should include all column GIDs (including imported ghosts) that are on the node.
00598 */
00599 
00600   // build the row map containing all the nodes (original matrix + off-node matrix)
00601   vector<int> rowList(NumMyRowsA_ + ghostElements.size());
00602   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
00603     rowList[i] = A().RowMatrixRowMap().GID(i);
00604   for (int i = 0 ; i < (int)ghostElements.size() ; ++i)
00605     rowList[i + NumMyRowsA_] = ghostElements[i];
00606 
00607   // row map for the overlapped matrix (local + overlap)
00608   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ghostElements.size(), &rowList[0], 0, Comm()) );
00609 
00610   // build the column map for the overlapping matrix
00611   //vector<int> colList(colMapElements.size());
00612   // column map for the overlapped matrix (local + overlap)
00613   //colMap_ = rcp( new Epetra_Map(-1, colMapElements.size(), &colList[0], 0, Comm()) );
00614   //for (int i = 0 ; i < (int)colMapElements.size() ; i++)
00615   //  colList[i] = colMapElements[i];
00616   vector<int> colList(A().RowMatrixColMap().NumMyElements() + colMapElements.size());
00617   int nc = A().RowMatrixColMap().NumMyElements();
00618   for (int i = 0 ; i < nc; i++)
00619     colList[i] = A().RowMatrixColMap().GID(i);
00620   for (int i = 0 ; i < (int)colMapElements.size() ; i++)
00621     colList[nc+i] = colMapElements[i];
00622 
00623   // column map for the overlapped matrix (local + overlap)
00624   //colMap_ = rcp( new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm()) );
00625   colMap_ = new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm());
00626     
00627 
00628   //FIXME timing
00629 #ifdef IFPACK_OVA_TIME_BUILD
00630   t1 = MPI_Wtime();
00631   fprintf(stderr,"[node %d]: time to init B() col/row maps %2.3e\n", nodeID, t1-t0);
00632   t0=t1;
00633 #endif
00634   //FIXME timing
00635 
00636   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00637   //++++ start of sort
00638   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00639   // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
00640   // different from that of Matrix.
00641   try {
00642     // build row map based on the local communicator.  We need this temporarily to build the column map.
00643     Epetra_Map* nodeMap = new Epetra_Map(-1,NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
00644     int numMyElts = colMap_->NumMyElements();
00645     assert(numMyElts!=0);
00646 
00647     // The column map *must* be sorted: first locals, then ghosts.  The ghosts
00648     // must be further sorted so that they are contiguous by owning processor.
00649 
00650     // For each GID in column map, determine owning PID in local communicator.
00651     int* myGlobalElts = new int[numMyElts];
00652     colMap_->MyGlobalElements(myGlobalElts);
00653     int* pidList = new int[numMyElts];
00654     nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
00655 
00656     // First sort column map GIDs based on the owning PID in local communicator.
00657     Epetra_Util Util;
00658     int *tt[1];
00659     tt[0] = myGlobalElts;
00660     Util.Sort(true, numMyElts, pidList, 0, (double**)0, 1, tt);
00661 
00662     // For each remote PID, sort the column map GIDs in ascending order.
00663     // Don't sort the local PID's entries.
00664     int localStart=0;
00665     while (localStart<numMyElts) {
00666       int currPID = (pidList)[localStart];
00667       int i=localStart;
00668       while (i<numMyElts && (pidList)[i] == currPID) i++;
00669       if (currPID != nodeComm->MyPID())
00670         Util.Sort(true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
00671       localStart = i;
00672     }
00673 
00674     // now move the local entries to the front of the list
00675     localStart=0;
00676     while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
00677     assert(localStart != numMyElts);
00678     int localEnd=localStart;
00679     while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
00680     int* mySortedGlobalElts = new int[numMyElts];
00681     int leng = localEnd - localStart;
00682     /* This is a little gotcha.  It appears that the ordering of the column map's local entries
00683        must be the same as that of the domain map's local entries.  See the comment in method
00684        MakeColMap() in Epetra_CrsGraph.cpp, line 1072. */
00685     int *rowGlobalElts =  nodeMap->MyGlobalElements();
00686 
00687     /*TODO TODO TODO NTS rows 68 and 83 show up as local GIDs in rowGlobalElts for both pids 0 & 1 on node 0.
00688                     This seems to imply that there is something wrong with rowList!!
00689                     It is ok for the overlapped matrix's row map to have repeated entries (it's overlapped, after all).
00690                     But ... on a node, there should be no repeats.
00691                     Ok, here's the underlying cause.  On node 0, gpid 1 finds 83 on overlap round 0.  gpid 0 finds 83
00692                     on overlap round 1.  This shouldn't happen .. once a node buddy "owns" a row, no other node buddy
00693                     should ever import ("own") that row.
00694 
00695                     Here's a possible fix.  Extend tie breaking across overlap rounds.  This means sending to lpid 0
00696                     a list of *all* rows you've imported (this round plus previous rounds) for tie-breaking
00697                     If a nb won a ghost row in a previous round, his votes for that ghost row should be really high, i.e.,
00698                     he should probably win that ghost row forever.
00699                     7/17/09
00700 
00701       7/28/09 JJH I believe I've fixed this, but I thought it might be helpful to have this comment in here, in case I missed something.
00702     */
00703 
00704     //move locals to front of list
00705     for (int i=0; i<leng; i++) {
00706       /* //original code */
00707       (mySortedGlobalElts)[i] = rowGlobalElts[i];
00708       //(&*mySortedGlobalElts)[i] = (&*myGlobalElts)[localStart+i];
00709       //(&*mySortedPidList)[i] = (&*pidList)[localStart+i];
00710     }
00711     for (int i=leng; i< localEnd; i++) {
00712       (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
00713     }
00714     for (int i=localEnd; i<numMyElts; i++) {
00715       (mySortedGlobalElts)[i] = (myGlobalElts)[i];
00716     }
00717 
00718     //FIXME timing
00719 #ifdef IFPACK_OVA_TIME_BUILD
00720     t1 = MPI_Wtime();
00721     fprintf(stderr,"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", nodeID, t1-t0);
00722     t0=t1;
00723 #endif
00724     //FIXME timing
00725 
00726     int indexBase = colMap_->IndexBase();
00727     if (colMap_) delete colMap_;
00728     colMap_ = new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,Comm());
00729 
00730     delete nodeMap;
00731     delete [] myGlobalElts;
00732     delete [] pidList;
00733     delete [] mySortedGlobalElts;
00734 
00735   } //try
00736   catch(...) {
00737     printf("** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
00738   }
00739 
00740   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00741   //++++ end of sort
00742   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00743 
00744   //FIXME timing
00745 #ifdef IFPACK_OVA_TIME_BUILD
00746   t1 = MPI_Wtime();
00747   fprintf(stderr,"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", nodeID, t1-t0);
00748   t0=t1;
00749 #endif
00750   //FIXME timing
00751 
00752 
00753 /*
00754 
00755    FIXME
00756    Does the column map need to be sorted for the overlapping matrix?
00757 
00758    The column map *must* be sorted:
00759 
00760         first locals
00761         then ghosts
00762 
00763    The ghosts must be further sorted so that they are contiguous by owning processor
00764 
00765   int* RemoteSizeList = 0
00766   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
00767 
00768   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
00769   Epetra_Util epetraUtil;
00770   SortLists[0] = RemoteColIndices;
00771   SortLists[1] = RemoteSizeList;
00772   epetraUtil.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
00773 */
00774 
00775   // now build the map corresponding to all the external nodes
00776   // (with respect to A().RowMatrixRowMap().
00777   ExtMap_ = rcp( new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,Comm()) );
00778   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*colMap_,0) ); 
00779 
00780   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
00781   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
00782 
00783   //FIXME timing
00784 #ifdef IFPACK_OVA_TIME_BUILD
00785   t1 = MPI_Wtime();
00786   fprintf(stderr,"[node %d]: time to import and setup ExtMap_ %2.3e\n", nodeID, t1-t0);
00787   t0=t1;
00788 #endif
00789   //FIXME timing
00790 
00791 
00792 /*
00793   Notes to self:    In FillComplete, the range map does not have to be 1-1 as long as
00794                     (row map == range map).  Ditto for the domain map not being 1-1
00795                     if (col map == domain map).
00796                     
00797 */
00798 
00799   ExtMatrix_->FillComplete( *colMap_ , *Map_ ); //FIXME wrong
00800 
00801   //FIXME timing
00802 #ifdef IFPACK_OVA_TIME_BUILD
00803   t1 = MPI_Wtime();
00804   fprintf(stderr,"[node %d]: time to FillComplete B() %2.3e\n", nodeID, t1-t0);
00805   t0=t1;
00806 #endif
00807   //FIXME timing
00808 
00809 
00810   // Note: B() = *ExtMatrix_ .
00811 
00812   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) ); //FIXME is this right?!
00813 
00814   // fix indices for overlapping matrix
00815   NumMyRowsB_ = B().NumMyRows();
00816   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;  //TODO is this wrong for a subdomain on >1 processor? // should be ok
00817   //NumMyCols_ = NumMyRows_;  //TODO is this wrong for a subdomain on >1 processor?  // YES!!!
00818   //NumMyCols_ = A().NumMyCols() + B().NumMyCols();
00819   NumMyCols_ = B().NumMyCols();
00820 
00821   /*FIXME*/ //somehow B's NumMyCols is the entire subdomain (local + overlap)
00822 
00823   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
00824   
00825   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
00826   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00827   MaxNumEntries_ = A().MaxNumEntries();
00828   
00829   if (MaxNumEntries_ < B().MaxNumEntries())
00830     MaxNumEntries_ = B().MaxNumEntries();
00831 
00832 # ifdef HAVE_MPI
00833   delete nodeComm;
00834 # endif
00835 
00836   //FIXME timing
00837 #ifdef IFPACK_OVA_TIME_BUILD
00838   t1 = MPI_Wtime();
00839   fprintf(stderr,"[node %d]: time for final calcs %2.3e\n", nodeID, t1-t0);
00840   fprintf(stderr,"[node %d]: total IORM ctor time %2.3e\n", nodeID, t1-true_t0);
00841 #endif
00842   //FIXME timing
00843 
00844 
00845 } //Ifpack_OverlappingRowMatrix() ctor for more than one core
00846 
00847 // Destructor
00848 Ifpack_OverlappingRowMatrix::~Ifpack_OverlappingRowMatrix() {
00849   delete colMap_;
00850 }
00851 
00852 #endif //ifdef IFPACK_NODE_AWARE_CODE
00853 
00854 
00855 // ======================================================================
00856 int Ifpack_OverlappingRowMatrix::
00857 NumMyRowEntries(int MyRow, int & NumEntries) const
00858 {
00859   if (MyRow < NumMyRowsA_)
00860     return(A().NumMyRowEntries(MyRow,NumEntries));
00861   else
00862     return(B().NumMyRowEntries(MyRow - NumMyRowsA_, NumEntries));
00863 }
00864 
00865 #ifdef IFPACK_NODE_AWARE_CODE
00866 // ======================================================================
00867 int Ifpack_OverlappingRowMatrix::
00868 ExtractMyRowCopy(int LocRow, int Length, int & NumEntries, double *Values, 
00869                  int * Indices) const
00870 {
00871   assert(1==0);
00872   int ierr;
00873   const Epetra_Map *Themap;
00874   if (LocRow < NumMyRowsA_) {
00875     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
00876     Themap=&A().RowMatrixColMap();
00877   }
00878   else {
00879     ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
00880     Themap=&B().RowMatrixColMap();
00881   }
00882 
00883   IFPACK_RETURN(ierr);
00884 }
00885 
00886 // ======================================================================
00887 int Ifpack_OverlappingRowMatrix::
00888 ExtractGlobalRowCopy(int GlobRow, int Length, int & NumEntries, double *Values, 
00889                  int * Indices) const
00890 {
00891   int ierr;
00892   const Epetra_Map *Themap;
00893   int LocRow = A().RowMatrixRowMap().LID(GlobRow);
00894   if (LocRow < NumMyRowsA_ && LocRow != -1) { //TODO don't need to check less than nummyrows
00895     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
00896     Themap=&A().RowMatrixColMap();
00897   }
00898   else {
00899     LocRow = B().RowMatrixRowMap().LID(GlobRow);
00900     assert(LocRow!=-1);
00901     //ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
00902     ierr = B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
00903     Themap=&B().RowMatrixColMap();
00904   }
00905 
00906   for (int i=0; i<NumEntries; i++) {
00907     Indices[i]=Themap->GID(Indices[i]);
00908     assert(Indices[i] != -1);
00909   }
00910 
00911   IFPACK_RETURN(ierr);
00912 }
00913 #else
00914 
00915 // ======================================================================
00916 int Ifpack_OverlappingRowMatrix::
00917 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values,
00918                  int * Indices) const
00919 {
00920   int ierr;
00921   if (MyRow < NumMyRowsA_)
00922     ierr = A().ExtractMyRowCopy(MyRow,Length,NumEntries,Values,Indices);
00923   else
00924     ierr = B().ExtractMyRowCopy(MyRow - NumMyRowsA_,Length,NumEntries,
00925                                 Values,Indices);
00926 
00927   IFPACK_RETURN(ierr);
00928 }
00929 #endif
00930 
00931 // ======================================================================
00932 int Ifpack_OverlappingRowMatrix::
00933 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00934 {
00935   IFPACK_CHK_ERR(-1);
00936 }
00937 
00938 
00939 // ======================================================================
00940 int Ifpack_OverlappingRowMatrix::
00941 Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00942 {
00943   int NumVectors = X.NumVectors();
00944   vector<int> Ind(MaxNumEntries_);
00945   vector<double> Val(MaxNumEntries_);
00946 
00947   Y.PutScalar(0.0);
00948 
00949   // matvec with A (local rows)
00950   for (int i = 0 ; i < NumMyRowsA_ ; ++i) {
00951     for (int k = 0 ; k < NumVectors ; ++k) {
00952       int Nnz;
00953       IFPACK_CHK_ERR(A().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
00954                                           &Val[0], &Ind[0]));
00955       for (int j = 0 ; j < Nnz ; ++j) {
00956         Y[k][i] += Val[j] * X[k][Ind[j]];
00957       }
00958     }
00959   }
00960 
00961   // matvec with B (overlapping rows)
00962   for (int i = 0 ; i < NumMyRowsB_ ; ++i) {
00963     for (int k = 0 ; k < NumVectors ; ++k) {
00964       int Nnz;
00965       IFPACK_CHK_ERR(B().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
00966                                           &Val[0], &Ind[0]));
00967       for (int j = 0 ; j < Nnz ; ++j) {
00968         Y[k][i + NumMyRowsA_] += Val[j] * X[k][Ind[j]];
00969       }
00970     }
00971   }
00972   return(0);
00973 }
00974 
00975 // ======================================================================
00976 int Ifpack_OverlappingRowMatrix::
00977 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00978 {
00979   IFPACK_CHK_ERR(Multiply(UseTranspose(),X,Y));
00980   return(0);
00981 }
00982 
00983 // ======================================================================
00984 int Ifpack_OverlappingRowMatrix::
00985 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00986 {
00987   IFPACK_CHK_ERR(-1);
00988 }
00989 
00990 // ======================================================================
00991 #ifndef IFPACK_NODE_AWARE_CODE
00992 Epetra_RowMatrix& Ifpack_OverlappingRowMatrix::B() const
00993 {
00994   return(*ExtMatrix_);
00995 }
00996 #endif
00997 // ======================================================================
00998 const Epetra_BlockMap& Ifpack_OverlappingRowMatrix::Map() const
00999 {
01000   return(*Map_);
01001 }
01002 
01003 // ======================================================================
01004 int Ifpack_OverlappingRowMatrix::
01005 ImportMultiVector(const Epetra_MultiVector& X, Epetra_MultiVector& OvX,
01006                   Epetra_CombineMode CM)
01007 {
01008   OvX.Import(X,*Importer_,CM);
01009   return(0);
01010 }
01011 
01012 // ======================================================================
01013 int Ifpack_OverlappingRowMatrix::
01014 ExportMultiVector(const Epetra_MultiVector& OvX, Epetra_MultiVector& X,
01015                   Epetra_CombineMode CM)
01016 {
01017   X.Export(OvX,*Importer_,CM);
01018   return(0);
01019 }
01020 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines