IFPACK Development
Ifpack_OverlappingRowMatrix.cpp
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_SUBCOMM_CODE
00039 #include "Epetra_IntVector.h"
00040 #include "Epetra_MpiComm.h"
00041 #include "Teuchos_Hashtable.hpp"
00042 #include "Teuchos_Array.hpp"
00043 #include "EpetraExt_OperatorOut.h"
00044 #include "EpetraExt_BlockMapOut.h"
00045 #else
00046 # ifdef IFPACK_NODE_AWARE_CODE
00047 # include "Epetra_IntVector.h"
00048 # include "Epetra_MpiComm.h"
00049 # include "Teuchos_Hashtable.hpp"
00050 # include "Teuchos_Array.hpp"
00051 # include "EpetraExt_OperatorOut.h"
00052 # include "EpetraExt_BlockMapOut.h"
00053   extern int ML_NODE_ID;
00054   int IFPACK_NODE_ID;
00055 # endif
00056 #endif
00057 
00058 using namespace Teuchos;
00059 
00060 // ====================================================================== 
00061 // Constructor for the case of one core per subdomain
00062 Ifpack_OverlappingRowMatrix::
00063 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00064                             int OverlapLevel_in)  :
00065   Matrix_(Matrix_in),
00066   OverlapLevel_(OverlapLevel_in)
00067 {
00068   // should not be here if no overlap
00069   if (OverlapLevel_in == 0)
00070     IFPACK_CHK_ERRV(-1);
00071 
00072   // nothing to do as well with one process
00073   if (Comm().NumProc() == 1)
00074     IFPACK_CHK_ERRV(-1);
00075   
00076   NumMyRowsA_ = A().NumMyRows();
00077 
00078   // construct the external matrix
00079   vector<int> ExtElements; 
00080 
00081   RCP<Epetra_Map> TmpMap;
00082   RCP<Epetra_CrsMatrix> TmpMatrix; 
00083   RCP<Epetra_Import> TmpImporter;
00084 
00085   // importing rows corresponding to elements that are 
00086   // in ColMap, but not in RowMap 
00087   const Epetra_Map *RowMap; 
00088   const Epetra_Map *ColMap; 
00089 
00090   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap) {
00091     if (TmpMatrix != Teuchos::null) {
00092       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00093       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00094     }
00095     else {
00096       RowMap = &(A().RowMatrixRowMap()); 
00097       ColMap = &(A().RowMatrixColMap()); 
00098     }
00099 
00100     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00101     vector<int> list(size); 
00102 
00103     int count = 0; 
00104 
00105     // define the set of rows that are in ColMap but not in RowMap 
00106     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) { 
00107       int GID = ColMap->GID(i); 
00108       if (A().RowMatrixRowMap().LID(GID) == -1) { 
00109         vector<int>::iterator pos 
00110           = find(ExtElements.begin(),ExtElements.end(),GID); 
00111         if (pos == ExtElements.end()) { 
00112           ExtElements.push_back(GID);
00113           list[count] = GID; 
00114           ++count; 
00115         } 
00116       } 
00117     } 
00118 
00119     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) ); 
00120 
00121     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
00122 
00123     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
00124 
00125     TmpMatrix->Import(A(),*TmpImporter,Insert); 
00126     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
00127 
00128   }
00129 
00130   // build the map containing all the nodes (original
00131   // matrix + extended matrix)
00132   vector<int> list(NumMyRowsA_ + ExtElements.size());
00133   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
00134     list[i] = A().RowMatrixRowMap().GID(i);
00135   for (int i = 0 ; i < (int)ExtElements.size() ; ++i)
00136     list[i + NumMyRowsA_] = ExtElements[i];
00137 
00138   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ExtElements.size(),
00139                       &list[0], 0, Comm()) );
00140 #ifdef IFPACK_SUBCOMM_CODE
00141   colMap_ = &*Map_;
00142 #else
00143 # ifdef IFPACK_NODE_AWARE_CODE
00144   colMap_ = &*Map_;
00145 # endif
00146 #endif
00147   // now build the map corresponding to all the external nodes
00148   // (with respect to A().RowMatrixRowMap().
00149   ExtMap_ = rcp( new Epetra_Map(-1,ExtElements.size(),
00150                      &ExtElements[0],0,A().Comm()) );
00151   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*Map_,0) ); 
00152 
00153   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
00154   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
00155   ExtMatrix_->FillComplete(A().OperatorDomainMap(),*Map_);
00156 
00157   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) );
00158 
00159   // fix indices for overlapping matrix
00160   NumMyRowsB_ = B().NumMyRows();
00161   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;
00162   NumMyCols_ = NumMyRows_;
00163   
00164   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
00165   
00166   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
00167   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00168   MaxNumEntries_ = A().MaxNumEntries();
00169   
00170   if (MaxNumEntries_ < B().MaxNumEntries())
00171     MaxNumEntries_ = B().MaxNumEntries();
00172 
00173 }
00174 
00175 #ifdef IFPACK_SUBCOMM_CODE
00176 
00177 // ====================================================================== 
00178 // Constructor for the case of two or more cores per subdomain
00179 Ifpack_OverlappingRowMatrix::
00180 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00181                             int OverlapLevel_in, int subdomainID)  :
00182   Matrix_(Matrix_in),
00183   OverlapLevel_(OverlapLevel_in)
00184 {
00185 
00186   //FIXME timing
00187 #ifdef IFPACK_OVA_TIME_BUILD
00188   double t0 = MPI_Wtime();
00189   double t1, true_t0=t0;
00190 #endif
00191   //FIXME timing
00192 
00193   const int votesMax = INT_MAX;
00194 
00195   // should not be here if no overlap
00196   if (OverlapLevel_in == 0)
00197     IFPACK_CHK_ERRV(-1);
00198 
00199   // nothing to do as well with one process
00200   if (Comm().NumProc() == 1)
00201     IFPACK_CHK_ERRV(-1);
00202 
00203   // subdomainID is the node (aka socket) number, and so is system dependent
00204   // nodeComm is the communicator for all the processes on a particular node
00205   // these processes will have the same subdomainID.
00206 # ifdef HAVE_MPI
00207   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &Comm() );
00208   assert(pComm != NULL);
00209   MPI_Comm nodeMPIComm;
00210   MPI_Comm_split(pComm->Comm(),subdomainID,pComm->MyPID(),&nodeMPIComm);
00211   Epetra_MpiComm *nodeComm = new Epetra_MpiComm(nodeMPIComm);
00212 # else
00213   Epetra_SerialComm *nodeComm =  dynamic_cast<Epetra_MpiComm*>( &(Matrix_in->RowMatrixRowMap().Comm()) );
00214 # endif
00215   
00216   NumMyRowsA_ = A().NumMyRows();
00217 
00218   RCP<Epetra_Map> TmpMap;
00219   RCP<Epetra_CrsMatrix> TmpMatrix; 
00220   RCP<Epetra_Import> TmpImporter;
00221 
00222   // importing rows corresponding to elements that are 
00223   // in ColMap, but not in RowMap 
00224   const Epetra_Map *RowMap; 
00225   const Epetra_Map *ColMap; 
00226   const Epetra_Map *DomainMap;
00227 
00228   // TODO Count #connections from nodes I own to each ghost node
00229 
00230 
00231   //FIXME timing
00232 #ifdef IFPACK_OVA_TIME_BUILD
00233   t1 = MPI_Wtime(); 
00234   fprintf(stderr,"[node %d]: time for initialization %2.3e\n", subdomainID, t1-t0);
00235   t0=t1;
00236 #endif
00237   //FIXME timing
00238 
00239 #ifdef IFPACK_OVA_TIME_BUILD
00240   t1 = MPI_Wtime();
00241   fprintf(stderr,"[node %d]: overlap hash table allocs %2.3e\n", subdomainID, t1-t0);
00242   t0=t1;
00243 #endif
00244 
00245   Teuchos::Hashtable<int,int> colMapTable(3 * A().RowMatrixColMap().NumMyElements() );
00246 
00247   // ghostTable holds off-node GIDs that are connected to on-node rows and can potentially be this PID's overlap
00248   // TODO hopefully 3 times the # column entries is a reasonable table size
00249   Teuchos::Hashtable<int,int> ghostTable(3 * A().RowMatrixColMap().NumMyElements() );
00250 
00251   /* ** ************************************************************************** ** */
00252   /* ** ********************** start of main overlap loop ************************ ** */
00253   /* ** ************************************************************************** ** */
00254   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00255   {
00256     // nbTable holds node buddy GIDs that are connected to current PID's rows, i.e., GIDs that should be in the overlapped
00257     // matrix's column map
00258 
00259     if (TmpMatrix != Teuchos::null) {
00260       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00261       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00262       DomainMap = &(TmpMatrix->OperatorDomainMap());
00263     }
00264     else {
00265       RowMap = &(A().RowMatrixRowMap()); 
00266       ColMap = &(A().RowMatrixColMap()); 
00267       DomainMap = &(A().OperatorDomainMap());
00268     }
00269 
00270 #ifdef IFPACK_OVA_TIME_BUILD
00271     t1 = MPI_Wtime();
00272     fprintf(stderr,"[node %d]: overlap pointer pulls %2.3e\n", subdomainID, t1-t0);
00273     t0=t1;
00274 #endif
00275     
00276     // For each column ID, determine the owning node (as opposed to core)
00277     // ID of the corresponding row.
00278     Epetra_IntVector colIdList( *ColMap );
00279     Epetra_IntVector rowIdList(*DomainMap);
00280     rowIdList.PutValue(subdomainID);  
00281     Teuchos::RCP<Epetra_Import> nodeIdImporter = rcp(new Epetra_Import( *ColMap, *DomainMap ));
00282 
00283 #ifdef IFPACK_OVA_TIME_BUILD
00284     t1 = MPI_Wtime();
00285     fprintf(stderr,"[node %d]: overlap intvector/importer allocs %2.3e\n", subdomainID, t1-t0);
00286     t0=t1;
00287 #endif
00288     
00289     colIdList.Import(rowIdList,*nodeIdImporter,Insert);
00290     
00291     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00292     int count = 0; 
00293 
00294 #ifdef IFPACK_OVA_TIME_BUILD
00295     t1 = MPI_Wtime();
00296     fprintf(stderr,"[node %d]: overlap col/row ID imports %2.3e\n", subdomainID, t1-t0);
00297     t0=t1;
00298 #endif
00299 
00300     
00301     // define the set of off-node rows that are in ColMap but not in RowMap
00302     // This naturally does not include off-core rows that are on the same node as me, i.e., node buddy rows.
00303     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
00304       int GID = ColMap->GID(i); 
00305       int myvotes=0;
00306       if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
00307       if ( colIdList[i] != subdomainID && myvotes < votesMax) // don't include anybody found in a previous round
00308       {
00309         int votes;
00310         if (ghostTable.containsKey(GID)) {
00311           votes = ghostTable.get(GID);
00312           votes++;
00313           ghostTable.put(GID,votes);
00314         } else {
00315           ghostTable.put(GID,1);
00316         }
00317       }
00318     } //for (int i = 0 ; i < ColMap->NumMyElements() ; ++i)
00319 
00320     Teuchos::Array<int> gidsarray,votesarray;
00321     ghostTable.arrayify(gidsarray,votesarray);
00322     int *gids = gidsarray.getRawPtr();
00323     int *votes = votesarray.getRawPtr();
00324 
00325     /*
00326        This next bit of code decides which node buddy (NB) gets which ghost points.  Everyone sends their
00327        list of ghost points to pid 0 of the local subcommunicator.  Pid 0 decides who gets what:
00328 
00329           - if a ghost point is touched by only one NB, that NB gets the ghost point
00330           - if two or more NBs share a ghost point, whichever NB has the most connections to the ghost
00331             point gets it.
00332     */
00333 
00334 #   ifdef HAVE_MPI  //FIXME What if we build in serial?!  This file will likely not compile.
00335     int *cullList;
00336     int ncull;
00337     //int mypid = nodeComm->MyPID();
00338 
00339     //FIXME timing
00340 #ifdef IFPACK_OVA_TIME_BUILD
00341     t1 = MPI_Wtime();
00342     fprintf(stderr,"[node %d]: overlap before culling %2.3e\n", subdomainID, t1-t0);
00343     t0=t1;
00344 #endif
00345     //FIXME timing
00346 
00347 
00348     if (nodeComm->MyPID() == 0)
00349     {
00350       // Figure out how much pid 0 is to receive
00351       MPI_Status status;
00352       int *lengths = new int[nodeComm->NumProc()+1];
00353       lengths[0] = 0;
00354       lengths[1] = ghostTable.size();
00355       for (int i=1; i<nodeComm->NumProc(); i++) {
00356         int leng;
00357         MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00358         lengths[i+1] = lengths[i] + leng;
00359       }
00360       int total = lengths[nodeComm->NumProc()];
00361 
00362       int* ghosts = new int[total];
00363       for (int i=0; i<total; i++) ghosts[i] = -9;
00364       int *round  = new int[total];
00365       int *owningpid  = new int[total];
00366 
00367       for (int i=1; i<nodeComm->NumProc(); i++) {
00368         int count = lengths[i+1] - lengths[i];
00369         MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00370         MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
00371       }
00372 
00373       // put in pid 0's info
00374       for (int i=0; i<lengths[1]; i++) {
00375         ghosts[i] = gids[i];
00376         round[i] = votes[i];
00377         owningpid[i] = 0;
00378       }
00379 
00380       // put in the pid associated with each ghost
00381       for (int j=1; j<nodeComm->NumProc(); j++) {
00382         for (int i=lengths[j]; i<lengths[j+1]; i++) {
00383           owningpid[i] = j;
00384         }
00385       }
00386 
00387       // sort everything based on the ghost gids
00388       int* roundpid[2];
00389       roundpid[0] = round; roundpid[1] = owningpid;
00390       Epetra_Util epetraUtil;
00391       epetraUtil.Sort(true,total,ghosts,0,0,2,roundpid);
00392 
00393       //set up arrays that get sent back to node buddies and that tell them what ghosts to cull
00394       int *nlosers = new int[nodeComm->NumProc()];
00395       int** losers = new int*[nodeComm->NumProc()];
00396       for (int i=0; i<nodeComm->NumProc(); i++) {
00397         nlosers[i] = 0;
00398         losers[i] = new int[ lengths[i+1]-lengths[i] ];
00399       }
00400 
00401       // Walk through ghosts array and and for each sequence of duplicate ghost GIDs, choose just one NB to keep it.
00402       // The logic is pretty simple.  The ghost list is sorted, so all duplicate PIDs are together.
00403       // The list is traversed.  As duplicates are found, node pid 0 keeps track of the current "winning"
00404       // pid.  When a pid is determined to have "lost" (less votes/connections to the current GID), the
00405       // GID is added to that pid's list of GIDs to be culled.  At the end of the repeated sequence, we have
00406       // a winner, and other NBs know whether they need to delete it from their import list.
00407       int max=0;   //for duplicated ghosts, index of pid with most votes and who hence keeps the ghost.
00408                    // TODO to break ties randomly
00409 
00410       for (int i=1; i<total; i++) {
00411         if (ghosts[i] == ghosts[i-1]) {
00412           int idx = i; // pid associated with idx is current "loser"
00413           if (round[i] > round[max]) {
00414             idx = max;
00415             max=i;
00416           }
00417           int j = owningpid[idx];
00418           losers[j][nlosers[j]++] = ghosts[idx];
00419         } else {
00420           max=i;
00421         }
00422       } //for (int i=1; i<total; i++)
00423 
00424       delete [] round;
00425       delete [] ghosts;
00426       delete [] owningpid;
00427 
00428       // send the arrays of ghost GIDs to be culled back to the respective node buddies
00429       for (int i=1; i<nodeComm->NumProc(); i++) {
00430         MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
00431         MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
00432       }
00433 
00434       //FIXME Unnecessary memory allocation and copying, but makes culling code cleaner
00435       //Could we stick this info into gids and votes, since neither is being used anymore?
00436       //TODO Instead of using "losers" arrays, just use "cullList" as in the else clause
00437       ncull = nlosers[0];
00438       cullList = new int[ncull+1];
00439       for (int i=0; i<nlosers[0]; i++)
00440         cullList[i] = losers[0][i];
00441 
00442       for (int i=0; i<nodeComm->NumProc(); i++)
00443         delete [] losers[i];
00444 
00445       delete [] lengths;
00446       delete [] losers;
00447       delete [] nlosers;
00448 
00449     } else { //everyone but pid 0
00450 
00451       // send to node pid 0 all ghosts that this pid could potentially import
00452       int hashsize = ghostTable.size();
00453       MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
00454       MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
00455       MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
00456 
00457       // receive the ghost GIDs that should not be imported (subset of the list sent off just above)
00458       MPI_Status status;
00459       MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
00460       cullList = new int[ncull+1];
00461       MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
00462     }
00463 
00464 
00465     //TODO clean out hash table after each time through overlap loop   4/1/07 JJH done moved both hash tables to local scope
00466 
00467     // Remove from my row map all off-node ghosts that will be imported by a node buddy.
00468     for (int i=0; i<ncull; i++) {
00469       try{ghostTable.remove(cullList[i]);}
00470 
00471       catch(...) {
00472         printf("pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
00473                Comm().MyPID(),cullList[i]);
00474         fflush(stdout);
00475       }
00476     }//for
00477 
00478     delete [] cullList;
00479 
00480     // Save off the remaining ghost GIDs from the current overlap round.
00481     // These are off-node GIDs (rows) that I will import.
00482     gidsarray.clear(); votesarray.clear();
00483     ghostTable.arrayify(gidsarray,votesarray);
00484 
00485     vector<int> list(size); 
00486     count=0;
00487     for (int i=0; i<ghostTable.size(); i++) {
00488       // if votesarray[i] == votesmax, then this GID was found during a previous overlap round
00489       if (votesarray[i] < votesMax) {
00490         list[count] = gidsarray[i];
00491         ghostTable.put(gidsarray[i],votesMax); //set this GID's vote to the maximum so that this PID alway wins in any subsequent overlap tiebreaking.
00492         count++;
00493       }
00494     }
00495 
00496     //FIXME timing
00497 #ifdef IFPACK_OVA_TIME_BUILD
00498     t1 = MPI_Wtime();
00499     fprintf(stderr,"[node %d]: overlap duplicate removal %2.3e\n", subdomainID, t1-t0);
00500     t0=t1;
00501 #endif
00502     //FIXME timing
00503 
00504 #   endif //ifdef HAVE_MPI
00505 
00506     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) );
00507 
00508     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
00509 
00510     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
00511 
00512     TmpMatrix->Import(A(),*TmpImporter,Insert); 
00513     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
00514 
00515     //FIXME timing
00516 #ifdef IFPACK_OVA_TIME_BUILD
00517     t1 = MPI_Wtime();
00518     fprintf(stderr,"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", subdomainID, t1-t0);
00519     t0=t1;
00520 #endif
00521     //FIXME timing
00522 
00523 
00524     // These next two imports get the GIDs that need to go into the column map of the overlapped matrix.
00525 
00526     // For each column ID in the overlap, determine the owning node (as opposed to core)
00527     // ID of the corresponding row.  Save those column IDs whose owning node is the current one.
00528     // This should get all the imported ghost GIDs.
00529     Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
00530     ov_colIdList.PutValue(-1);
00531     Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
00532     ov_rowIdList.PutValue(subdomainID);  
00533     Teuchos::RCP<Epetra_Import> ov_nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), TmpMatrix->RowMap()));
00534     ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
00535 
00536     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
00537        if (ov_colIdList[i] == subdomainID) {
00538          int GID = (ov_colIdList.Map()).GID(i);
00539          colMapTable.put(GID,1);
00540        }
00541     }
00542 
00543     // Do a second import of the owning node ID from A's rowmap to TmpMat's column map.  This ensures that
00544     // all GIDs that belong to a node buddy and are in a ghost row's sparsity pattern will be in the final
00545     // overlapped matrix's column map.
00546     ov_colIdList.PutValue(-1);
00547     Epetra_IntVector ArowIdList( A().RowMatrixRowMap() );
00548     ArowIdList.PutValue(subdomainID);
00549     nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), A().RowMatrixRowMap() ));
00550     ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
00551 
00552     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
00553        if (ov_colIdList[i] == subdomainID) {
00554          int GID = (ov_colIdList.Map()).GID(i);
00555          colMapTable.put(GID,1);
00556        }
00557     }
00558 
00559     //FIXME timing
00560 #ifdef IFPACK_OVA_TIME_BUILD
00561     t1 = MPI_Wtime();
00562     fprintf(stderr,"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", subdomainID, t1-t0);
00563     t0=t1;
00564 #endif
00565     //FIXME timing
00566 
00567   } //for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00568 
00569   /* ** ************************************************************************ ** */
00570   /* ** ********************** end of main overlap loop ************************ ** */
00571   /* ** ************************************************************************ ** */
00572 
00573   // off-node GIDs that will be in the overlap
00574   vector<int> ghostElements; 
00575 
00576   Teuchos::Array<int> gidsarray,votesarray;
00577   ghostTable.arrayify(gidsarray,votesarray);
00578   for (int i=0; i<ghostTable.size(); i++) {
00579     ghostElements.push_back(gidsarray[i]);
00580   }
00581 
00582     for (int i = 0 ; i < A().RowMatrixColMap().NumMyElements() ; ++i) {
00583       int GID = A().RowMatrixColMap().GID(i);
00584       // Remove any entries that are in A's original column map
00585       if (colMapTable.containsKey(GID)) {
00586         try{colMapTable.remove(GID);}
00587         catch(...) {
00588           printf("pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n", Comm().MyPID(),GID);
00589           fflush(stdout);
00590         }
00591       }
00592     }
00593 
00594     // GIDs that will be in the overlapped matrix's column map
00595     vector<int> colMapElements; 
00596 
00597   gidsarray.clear(); votesarray.clear();
00598   colMapTable.arrayify(gidsarray,votesarray);
00599   for (int i=0; i<colMapTable.size(); i++)
00600     colMapElements.push_back(gidsarray[i]);
00601 
00602 /*
00603    We need two maps here.  The first is the row map, which we've got by using my original rows
00604    plus whatever I've picked up in ghostElements.
00605 
00606    The second is a column map.  This map should include my rows, plus any columns that could come from node buddies.
00607    These GIDs come from the std:array colMapElements, which in turn comes from colMapTable.
00608    This map should *not* omit ghosts that have been imported by my node buddies, i.e., for any row that I own,
00609    the stencil should include all column GIDs (including imported ghosts) that are on the node.
00610 */
00611 
00612   // build the row map containing all the nodes (original matrix + off-node matrix)
00613   vector<int> rowList(NumMyRowsA_ + ghostElements.size());
00614   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
00615     rowList[i] = A().RowMatrixRowMap().GID(i);
00616   for (int i = 0 ; i < (int)ghostElements.size() ; ++i)
00617     rowList[i + NumMyRowsA_] = ghostElements[i];
00618 
00619   // row map for the overlapped matrix (local + overlap)
00620   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ghostElements.size(), &rowList[0], 0, Comm()) );
00621 
00622   // build the column map for the overlapping matrix
00623   //vector<int> colList(colMapElements.size());
00624   // column map for the overlapped matrix (local + overlap)
00625   //colMap_ = rcp( new Epetra_Map(-1, colMapElements.size(), &colList[0], 0, Comm()) );
00626   //for (int i = 0 ; i < (int)colMapElements.size() ; i++)
00627   //  colList[i] = colMapElements[i];
00628   vector<int> colList(A().RowMatrixColMap().NumMyElements() + colMapElements.size());
00629   int nc = A().RowMatrixColMap().NumMyElements();
00630   for (int i = 0 ; i < nc; i++)
00631     colList[i] = A().RowMatrixColMap().GID(i);
00632   for (int i = 0 ; i < (int)colMapElements.size() ; i++)
00633     colList[nc+i] = colMapElements[i];
00634 
00635   // column map for the overlapped matrix (local + overlap)
00636   //colMap_ = rcp( new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm()) );
00637   colMap_ = new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm());
00638     
00639 
00640   //FIXME timing
00641 #ifdef IFPACK_OVA_TIME_BUILD
00642   t1 = MPI_Wtime();
00643   fprintf(stderr,"[node %d]: time to init B() col/row maps %2.3e\n", subdomainID, t1-t0);
00644   t0=t1;
00645 #endif
00646   //FIXME timing
00647 
00648   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00649   //++++ start of sort
00650   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00651   // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
00652   // different from that of Matrix.
00653   try {
00654     // build row map based on the local communicator.  We need this temporarily to build the column map.
00655     Epetra_Map* nodeMap = new Epetra_Map(-1,NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
00656     int numMyElts = colMap_->NumMyElements();
00657     assert(numMyElts!=0);
00658 
00659     // The column map *must* be sorted: first locals, then ghosts.  The ghosts
00660     // must be further sorted so that they are contiguous by owning processor.
00661 
00662     // For each GID in column map, determine owning PID in local communicator.
00663     int* myGlobalElts = new int[numMyElts];
00664     colMap_->MyGlobalElements(myGlobalElts);
00665     int* pidList = new int[numMyElts];
00666     nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
00667 
00668     // First sort column map GIDs based on the owning PID in local communicator.
00669     Epetra_Util Util;
00670     int *tt[1];
00671     tt[0] = myGlobalElts;
00672     Util.Sort(true, numMyElts, pidList, 0, (double**)0, 1, tt);
00673 
00674     // For each remote PID, sort the column map GIDs in ascending order.
00675     // Don't sort the local PID's entries.
00676     int localStart=0;
00677     while (localStart<numMyElts) {
00678       int currPID = (pidList)[localStart];
00679       int i=localStart;
00680       while (i<numMyElts && (pidList)[i] == currPID) i++;
00681       if (currPID != nodeComm->MyPID())
00682         Util.Sort(true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
00683       localStart = i;
00684     }
00685 
00686     // now move the local entries to the front of the list
00687     localStart=0;
00688     while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
00689     assert(localStart != numMyElts);
00690     int localEnd=localStart;
00691     while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
00692     int* mySortedGlobalElts = new int[numMyElts];
00693     int leng = localEnd - localStart;
00694     /* This is a little gotcha.  It appears that the ordering of the column map's local entries
00695        must be the same as that of the domain map's local entries.  See the comment in method
00696        MakeColMap() in Epetra_CrsGraph.cpp, line 1072. */
00697     int *rowGlobalElts =  nodeMap->MyGlobalElements();
00698 
00699     /*TODO TODO TODO NTS rows 68 and 83 show up as local GIDs in rowGlobalElts for both pids 0 & 1 on node 0.
00700                     This seems to imply that there is something wrong with rowList!!
00701                     It is ok for the overlapped matrix's row map to have repeated entries (it's overlapped, after all).
00702                     But ... on a node, there should be no repeats.
00703                     Ok, here's the underlying cause.  On node 0, gpid 1 finds 83 on overlap round 0.  gpid 0 finds 83
00704                     on overlap round 1.  This shouldn't happen .. once a node buddy "owns" a row, no other node buddy
00705                     should ever import ("own") that row.
00706 
00707                     Here's a possible fix.  Extend tie breaking across overlap rounds.  This means sending to lpid 0
00708                     a list of *all* rows you've imported (this round plus previous rounds) for tie-breaking
00709                     If a nb won a ghost row in a previous round, his votes for that ghost row should be really high, i.e.,
00710                     he should probably win that ghost row forever.
00711                     7/17/09
00712 
00713       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.
00714     */
00715 
00716     //move locals to front of list
00717     for (int i=0; i<leng; i++) {
00718       /* //original code */
00719       (mySortedGlobalElts)[i] = rowGlobalElts[i];
00720       //(&*mySortedGlobalElts)[i] = (&*myGlobalElts)[localStart+i];
00721       //(&*mySortedPidList)[i] = (&*pidList)[localStart+i];
00722     }
00723     for (int i=leng; i< localEnd; i++) {
00724       (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
00725     }
00726     for (int i=localEnd; i<numMyElts; i++) {
00727       (mySortedGlobalElts)[i] = (myGlobalElts)[i];
00728     }
00729 
00730     //FIXME timing
00731 #ifdef IFPACK_OVA_TIME_BUILD
00732     t1 = MPI_Wtime();
00733     fprintf(stderr,"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", subdomainID, t1-t0);
00734     t0=t1;
00735 #endif
00736     //FIXME timing
00737 
00738     int indexBase = colMap_->IndexBase();
00739     if (colMap_) delete colMap_;
00740     colMap_ = new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,Comm());
00741 
00742     delete nodeMap;
00743     delete [] myGlobalElts;
00744     delete [] pidList;
00745     delete [] mySortedGlobalElts;
00746 
00747   } //try
00748   catch(...) {
00749     printf("** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
00750   }
00751 
00752   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00753   //++++ end of sort
00754   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00755 
00756   //FIXME timing
00757 #ifdef IFPACK_OVA_TIME_BUILD
00758   t1 = MPI_Wtime();
00759   fprintf(stderr,"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", subdomainID, t1-t0);
00760   t0=t1;
00761 #endif
00762   //FIXME timing
00763 
00764 
00765 /*
00766 
00767    FIXME
00768    Does the column map need to be sorted for the overlapping matrix?
00769 
00770    The column map *must* be sorted:
00771 
00772         first locals
00773         then ghosts
00774 
00775    The ghosts must be further sorted so that they are contiguous by owning processor
00776 
00777   int* RemoteSizeList = 0
00778   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
00779 
00780   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
00781   Epetra_Util epetraUtil;
00782   SortLists[0] = RemoteColIndices;
00783   SortLists[1] = RemoteSizeList;
00784   epetraUtil.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
00785 */
00786 
00787   // now build the map corresponding to all the external nodes
00788   // (with respect to A().RowMatrixRowMap().
00789   ExtMap_ = rcp( new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,Comm()) );
00790   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*colMap_,0) ); 
00791 
00792   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
00793   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
00794 
00795   //FIXME timing
00796 #ifdef IFPACK_OVA_TIME_BUILD
00797   t1 = MPI_Wtime();
00798   fprintf(stderr,"[node %d]: time to import and setup ExtMap_ %2.3e\n", subdomainID, t1-t0);
00799   t0=t1;
00800 #endif
00801   //FIXME timing
00802 
00803 
00804 /*
00805   Notes to self:    In FillComplete, the range map does not have to be 1-1 as long as
00806                     (row map == range map).  Ditto for the domain map not being 1-1
00807                     if (col map == domain map).
00808                     
00809 */
00810 
00811   ExtMatrix_->FillComplete( *colMap_ , *Map_ ); //FIXME wrong
00812 
00813   //FIXME timing
00814 #ifdef IFPACK_OVA_TIME_BUILD
00815   t1 = MPI_Wtime();
00816   fprintf(stderr,"[node %d]: time to FillComplete B() %2.3e\n", subdomainID, t1-t0);
00817   t0=t1;
00818 #endif
00819   //FIXME timing
00820 
00821 
00822   // Note: B() = *ExtMatrix_ .
00823 
00824   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) ); //FIXME is this right?!
00825 
00826   // fix indices for overlapping matrix
00827   NumMyRowsB_ = B().NumMyRows();
00828   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;  //TODO is this wrong for a subdomain on >1 processor? // should be ok
00829   //NumMyCols_ = NumMyRows_;  //TODO is this wrong for a subdomain on >1 processor?  // YES!!!
00830   //NumMyCols_ = A().NumMyCols() + B().NumMyCols();
00831   NumMyCols_ = B().NumMyCols();
00832 
00833   /*FIXME*/ //somehow B's NumMyCols is the entire subdomain (local + overlap)
00834 
00835   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
00836   
00837   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
00838   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00839   MaxNumEntries_ = A().MaxNumEntries();
00840   
00841   if (MaxNumEntries_ < B().MaxNumEntries())
00842     MaxNumEntries_ = B().MaxNumEntries();
00843 
00844 # ifdef HAVE_MPI
00845   delete nodeComm;
00846 # endif
00847 
00848   //FIXME timing
00849 #ifdef IFPACK_OVA_TIME_BUILD
00850   t1 = MPI_Wtime();
00851   fprintf(stderr,"[node %d]: time for final calcs %2.3e\n", subdomainID, t1-t0);
00852   fprintf(stderr,"[node %d]: total IORM ctor time %2.3e\n", subdomainID, t1-true_t0);
00853 #endif
00854   //FIXME timing
00855 
00856 
00857 } //Ifpack_OverlappingRowMatrix() ctor for more than one core
00858 
00859 #else
00860 # ifdef IFPACK_NODE_AWARE_CODE
00861 
00862 // ====================================================================== 
00863 // Constructor for the case of two or more cores per subdomain
00864 Ifpack_OverlappingRowMatrix::
00865 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00866                             int OverlapLevel_in, int nodeID)  :
00867   Matrix_(Matrix_in),
00868   OverlapLevel_(OverlapLevel_in)
00869 {
00870 
00871   //FIXME timing
00872 #ifdef IFPACK_OVA_TIME_BUILD
00873   double t0 = MPI_Wtime();
00874   double t1, true_t0=t0;
00875 #endif
00876   //FIXME timing
00877 
00878   const int votesMax = INT_MAX;
00879 
00880   // should not be here if no overlap
00881   if (OverlapLevel_in == 0)
00882     IFPACK_CHK_ERRV(-1);
00883 
00884   // nothing to do as well with one process
00885   if (Comm().NumProc() == 1)
00886     IFPACK_CHK_ERRV(-1);
00887 
00888   // nodeID is the node (aka socket) number, and so is system dependent
00889   // nodeComm is the communicator for all the processes on a particular node
00890   // these processes will have the same nodeID.
00891 # ifdef HAVE_MPI
00892   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &Comm() );
00893   assert(pComm != NULL);
00894   MPI_Comm nodeMPIComm;
00895   MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm);
00896   Epetra_MpiComm *nodeComm = new Epetra_MpiComm(nodeMPIComm);
00897 # else
00898   Epetra_SerialComm *nodeComm =  dynamic_cast<Epetra_MpiComm*>( &(Matrix_in->RowMatrixRowMap().Comm()) );
00899 # endif
00900   
00901   NumMyRowsA_ = A().NumMyRows();
00902 
00903   RCP<Epetra_Map> TmpMap;
00904   RCP<Epetra_CrsMatrix> TmpMatrix; 
00905   RCP<Epetra_Import> TmpImporter;
00906 
00907   // importing rows corresponding to elements that are 
00908   // in ColMap, but not in RowMap 
00909   const Epetra_Map *RowMap; 
00910   const Epetra_Map *ColMap; 
00911   const Epetra_Map *DomainMap;
00912 
00913   // TODO Count #connections from nodes I own to each ghost node
00914 
00915 
00916   //FIXME timing
00917 #ifdef IFPACK_OVA_TIME_BUILD
00918   t1 = MPI_Wtime(); 
00919   fprintf(stderr,"[node %d]: time for initialization %2.3e\n", nodeID, t1-t0);
00920   t0=t1;
00921 #endif
00922   //FIXME timing
00923 
00924 #ifdef IFPACK_OVA_TIME_BUILD
00925   t1 = MPI_Wtime();
00926   fprintf(stderr,"[node %d]: overlap hash table allocs %2.3e\n", nodeID, t1-t0);
00927   t0=t1;
00928 #endif
00929 
00930   Teuchos::Hashtable<int,int> colMapTable(3 * A().RowMatrixColMap().NumMyElements() );
00931 
00932   // ghostTable holds off-node GIDs that are connected to on-node rows and can potentially be this PID's overlap
00933   // TODO hopefully 3 times the # column entries is a reasonable table size
00934   Teuchos::Hashtable<int,int> ghostTable(3 * A().RowMatrixColMap().NumMyElements() );
00935 
00936   /* ** ************************************************************************** ** */
00937   /* ** ********************** start of main overlap loop ************************ ** */
00938   /* ** ************************************************************************** ** */
00939   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00940   {
00941     // nbTable holds node buddy GIDs that are connected to current PID's rows, i.e., GIDs that should be in the overlapped
00942     // matrix's column map
00943 
00944     if (TmpMatrix != Teuchos::null) {
00945       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00946       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00947       DomainMap = &(TmpMatrix->OperatorDomainMap());
00948     }
00949     else {
00950       RowMap = &(A().RowMatrixRowMap()); 
00951       ColMap = &(A().RowMatrixColMap()); 
00952       DomainMap = &(A().OperatorDomainMap());
00953     }
00954 
00955 #ifdef IFPACK_OVA_TIME_BUILD
00956     t1 = MPI_Wtime();
00957     fprintf(stderr,"[node %d]: overlap pointer pulls %2.3e\n", nodeID, t1-t0);
00958     t0=t1;
00959 #endif
00960     
00961     // For each column ID, determine the owning node (as opposed to core)
00962     // ID of the corresponding row.
00963     Epetra_IntVector colIdList( *ColMap );
00964     Epetra_IntVector rowIdList(*DomainMap);
00965     rowIdList.PutValue(nodeID);  
00966     Teuchos::RCP<Epetra_Import> nodeIdImporter = rcp(new Epetra_Import( *ColMap, *DomainMap ));
00967 
00968 #ifdef IFPACK_OVA_TIME_BUILD
00969     t1 = MPI_Wtime();
00970     fprintf(stderr,"[node %d]: overlap intvector/importer allocs %2.3e\n", nodeID, t1-t0);
00971     t0=t1;
00972 #endif
00973     
00974     colIdList.Import(rowIdList,*nodeIdImporter,Insert);
00975     
00976     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00977     int count = 0; 
00978 
00979 #ifdef IFPACK_OVA_TIME_BUILD
00980     t1 = MPI_Wtime();
00981     fprintf(stderr,"[node %d]: overlap col/row ID imports %2.3e\n", nodeID, t1-t0);
00982     t0=t1;
00983 #endif
00984 
00985     
00986     // define the set of off-node rows that are in ColMap but not in RowMap
00987     // This naturally does not include off-core rows that are on the same node as me, i.e., node buddy rows.
00988     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
00989       int GID = ColMap->GID(i); 
00990       int myvotes=0;
00991       if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
00992       if ( colIdList[i] != nodeID && myvotes < votesMax) // don't include anybody found in a previous round
00993       {
00994         int votes;
00995         if (ghostTable.containsKey(GID)) {
00996           votes = ghostTable.get(GID);
00997           votes++;
00998           ghostTable.put(GID,votes);
00999         } else {
01000           ghostTable.put(GID,1);
01001         }
01002       }
01003     } //for (int i = 0 ; i < ColMap->NumMyElements() ; ++i)
01004 
01005     Teuchos::Array<int> gidsarray,votesarray;
01006     ghostTable.arrayify(gidsarray,votesarray);
01007     int *gids = gidsarray.getRawPtr();
01008     int *votes = votesarray.getRawPtr();
01009 
01010     /*
01011        This next bit of code decides which node buddy (NB) gets which ghost points.  Everyone sends their
01012        list of ghost points to pid 0 of the local subcommunicator.  Pid 0 decides who gets what:
01013 
01014           - if a ghost point is touched by only one NB, that NB gets the ghost point
01015           - if two or more NBs share a ghost point, whichever NB has the most connections to the ghost
01016             point gets it.
01017     */
01018 
01019 #   ifdef HAVE_MPI  //FIXME What if we build in serial?!  This file will likely not compile.
01020     int *cullList;
01021     int ncull;
01022     //int mypid = nodeComm->MyPID();
01023 
01024     //FIXME timing
01025 #ifdef IFPACK_OVA_TIME_BUILD
01026     t1 = MPI_Wtime();
01027     fprintf(stderr,"[node %d]: overlap before culling %2.3e\n", nodeID, t1-t0);
01028     t0=t1;
01029 #endif
01030     //FIXME timing
01031 
01032 
01033     if (nodeComm->MyPID() == 0)
01034     {
01035       // Figure out how much pid 0 is to receive
01036       MPI_Status status;
01037       int *lengths = new int[nodeComm->NumProc()+1];
01038       lengths[0] = 0;
01039       lengths[1] = ghostTable.size();
01040       for (int i=1; i<nodeComm->NumProc(); i++) {
01041         int leng;
01042         MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01043         lengths[i+1] = lengths[i] + leng;
01044       }
01045       int total = lengths[nodeComm->NumProc()];
01046 
01047       int* ghosts = new int[total];
01048       for (int i=0; i<total; i++) ghosts[i] = -9;
01049       int *round  = new int[total];
01050       int *owningpid  = new int[total];
01051 
01052       for (int i=1; i<nodeComm->NumProc(); i++) {
01053         int count = lengths[i+1] - lengths[i];
01054         MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01055         MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01056       }
01057 
01058       // put in pid 0's info
01059       for (int i=0; i<lengths[1]; i++) {
01060         ghosts[i] = gids[i];
01061         round[i] = votes[i];
01062         owningpid[i] = 0;
01063       }
01064 
01065       // put in the pid associated with each ghost
01066       for (int j=1; j<nodeComm->NumProc(); j++) {
01067         for (int i=lengths[j]; i<lengths[j+1]; i++) {
01068           owningpid[i] = j;
01069         }
01070       }
01071 
01072       // sort everything based on the ghost gids
01073       int* roundpid[2];
01074       roundpid[0] = round; roundpid[1] = owningpid;
01075       Epetra_Util epetraUtil;
01076       epetraUtil.Sort(true,total,ghosts,0,0,2,roundpid);
01077 
01078       //set up arrays that get sent back to node buddies and that tell them what ghosts to cull
01079       int *nlosers = new int[nodeComm->NumProc()];
01080       int** losers = new int*[nodeComm->NumProc()];
01081       for (int i=0; i<nodeComm->NumProc(); i++) {
01082         nlosers[i] = 0;
01083         losers[i] = new int[ lengths[i+1]-lengths[i] ];
01084       }
01085 
01086       // Walk through ghosts array and and for each sequence of duplicate ghost GIDs, choose just one NB to keep it.
01087       // The logic is pretty simple.  The ghost list is sorted, so all duplicate PIDs are together.
01088       // The list is traversed.  As duplicates are found, node pid 0 keeps track of the current "winning"
01089       // pid.  When a pid is determined to have "lost" (less votes/connections to the current GID), the
01090       // GID is added to that pid's list of GIDs to be culled.  At the end of the repeated sequence, we have
01091       // a winner, and other NBs know whether they need to delete it from their import list.
01092       int max=0;   //for duplicated ghosts, index of pid with most votes and who hence keeps the ghost.
01093                    // TODO to break ties randomly
01094 
01095       for (int i=1; i<total; i++) {
01096         if (ghosts[i] == ghosts[i-1]) {
01097           int idx = i; // pid associated with idx is current "loser"
01098           if (round[i] > round[max]) {
01099             idx = max;
01100             max=i;
01101           }
01102           int j = owningpid[idx];
01103           losers[j][nlosers[j]++] = ghosts[idx];
01104         } else {
01105           max=i;
01106         }
01107       } //for (int i=1; i<total; i++)
01108 
01109       delete [] round;
01110       delete [] ghosts;
01111       delete [] owningpid;
01112 
01113       // send the arrays of ghost GIDs to be culled back to the respective node buddies
01114       for (int i=1; i<nodeComm->NumProc(); i++) {
01115         MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
01116         MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
01117       }
01118 
01119       //FIXME Unnecessary memory allocation and copying, but makes culling code cleaner
01120       //Could we stick this info into gids and votes, since neither is being used anymore?
01121       //TODO Instead of using "losers" arrays, just use "cullList" as in the else clause
01122       ncull = nlosers[0];
01123       cullList = new int[ncull+1];
01124       for (int i=0; i<nlosers[0]; i++)
01125         cullList[i] = losers[0][i];
01126 
01127       for (int i=0; i<nodeComm->NumProc(); i++)
01128         delete [] losers[i];
01129 
01130       delete [] lengths;
01131       delete [] losers;
01132       delete [] nlosers;
01133 
01134     } else { //everyone but pid 0
01135 
01136       // send to node pid 0 all ghosts that this pid could potentially import
01137       int hashsize = ghostTable.size();
01138       MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
01139       MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
01140       MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
01141 
01142       // receive the ghost GIDs that should not be imported (subset of the list sent off just above)
01143       MPI_Status status;
01144       MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
01145       cullList = new int[ncull+1];
01146       MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
01147     }
01148 
01149 
01150     //TODO clean out hash table after each time through overlap loop   4/1/07 JJH done moved both hash tables to local scope
01151 
01152     // Remove from my row map all off-node ghosts that will be imported by a node buddy.
01153     for (int i=0; i<ncull; i++) {
01154       try{ghostTable.remove(cullList[i]);}
01155 
01156       catch(...) {
01157         printf("pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
01158                Comm().MyPID(),cullList[i]);
01159         fflush(stdout);
01160       }
01161     }//for
01162 
01163     delete [] cullList;
01164 
01165     // Save off the remaining ghost GIDs from the current overlap round.
01166     // These are off-node GIDs (rows) that I will import.
01167     gidsarray.clear(); votesarray.clear();
01168     ghostTable.arrayify(gidsarray,votesarray);
01169 
01170     vector<int> list(size); 
01171     count=0;
01172     for (int i=0; i<ghostTable.size(); i++) {
01173       // if votesarray[i] == votesmax, then this GID was found during a previous overlap round
01174       if (votesarray[i] < votesMax) {
01175         list[count] = gidsarray[i];
01176         ghostTable.put(gidsarray[i],votesMax); //set this GID's vote to the maximum so that this PID alway wins in any subsequent overlap tiebreaking.
01177         count++;
01178       }
01179     }
01180 
01181     //FIXME timing
01182 #ifdef IFPACK_OVA_TIME_BUILD
01183     t1 = MPI_Wtime();
01184     fprintf(stderr,"[node %d]: overlap duplicate removal %2.3e\n", nodeID, t1-t0);
01185     t0=t1;
01186 #endif
01187     //FIXME timing
01188 
01189 #   endif //ifdef HAVE_MPI
01190 
01191     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) );
01192 
01193     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
01194 
01195     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
01196 
01197     TmpMatrix->Import(A(),*TmpImporter,Insert); 
01198     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
01199 
01200     //FIXME timing
01201 #ifdef IFPACK_OVA_TIME_BUILD
01202     t1 = MPI_Wtime();
01203     fprintf(stderr,"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", nodeID, t1-t0);
01204     t0=t1;
01205 #endif
01206     //FIXME timing
01207 
01208 
01209     // These next two imports get the GIDs that need to go into the column map of the overlapped matrix.
01210 
01211     // For each column ID in the overlap, determine the owning node (as opposed to core)
01212     // ID of the corresponding row.  Save those column IDs whose owning node is the current one.
01213     // This should get all the imported ghost GIDs.
01214     Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
01215     ov_colIdList.PutValue(-1);
01216     Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
01217     ov_rowIdList.PutValue(nodeID);  
01218     Teuchos::RCP<Epetra_Import> ov_nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), TmpMatrix->RowMap()));
01219     ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
01220 
01221     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
01222        if (ov_colIdList[i] == nodeID) {
01223          int GID = (ov_colIdList.Map()).GID(i);
01224          colMapTable.put(GID,1);
01225        }
01226     }
01227 
01228     // Do a second import of the owning node ID from A's rowmap to TmpMat's column map.  This ensures that
01229     // all GIDs that belong to a node buddy and are in a ghost row's sparsity pattern will be in the final
01230     // overlapped matrix's column map.
01231     ov_colIdList.PutValue(-1);
01232     Epetra_IntVector ArowIdList( A().RowMatrixRowMap() );
01233     ArowIdList.PutValue(nodeID);
01234     nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), A().RowMatrixRowMap() ));
01235     ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
01236 
01237     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
01238        if (ov_colIdList[i] == nodeID) {
01239          int GID = (ov_colIdList.Map()).GID(i);
01240          colMapTable.put(GID,1);
01241        }
01242     }
01243 
01244     //FIXME timing
01245 #ifdef IFPACK_OVA_TIME_BUILD
01246     t1 = MPI_Wtime();
01247     fprintf(stderr,"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", nodeID, t1-t0);
01248     t0=t1;
01249 #endif
01250     //FIXME timing
01251 
01252   } //for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
01253 
01254   /* ** ************************************************************************ ** */
01255   /* ** ********************** end of main overlap loop ************************ ** */
01256   /* ** ************************************************************************ ** */
01257 
01258   // off-node GIDs that will be in the overlap
01259   vector<int> ghostElements; 
01260 
01261   Teuchos::Array<int> gidsarray,votesarray;
01262   ghostTable.arrayify(gidsarray,votesarray);
01263   for (int i=0; i<ghostTable.size(); i++) {
01264     ghostElements.push_back(gidsarray[i]);
01265   }
01266 
01267     for (int i = 0 ; i < A().RowMatrixColMap().NumMyElements() ; ++i) {
01268       int GID = A().RowMatrixColMap().GID(i);
01269       // Remove any entries that are in A's original column map
01270       if (colMapTable.containsKey(GID)) {
01271         try{colMapTable.remove(GID);}
01272         catch(...) {
01273           printf("pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n", Comm().MyPID(),GID);
01274           fflush(stdout);
01275         }
01276       }
01277     }
01278 
01279     // GIDs that will be in the overlapped matrix's column map
01280     vector<int> colMapElements; 
01281 
01282   gidsarray.clear(); votesarray.clear();
01283   colMapTable.arrayify(gidsarray,votesarray);
01284   for (int i=0; i<colMapTable.size(); i++)
01285     colMapElements.push_back(gidsarray[i]);
01286 
01287 /*
01288    We need two maps here.  The first is the row map, which we've got by using my original rows
01289    plus whatever I've picked up in ghostElements.
01290 
01291    The second is a column map.  This map should include my rows, plus any columns that could come from node buddies.
01292    These GIDs come from the std:array colMapElements, which in turn comes from colMapTable.
01293    This map should *not* omit ghosts that have been imported by my node buddies, i.e., for any row that I own,
01294    the stencil should include all column GIDs (including imported ghosts) that are on the node.
01295 */
01296 
01297   // build the row map containing all the nodes (original matrix + off-node matrix)
01298   vector<int> rowList(NumMyRowsA_ + ghostElements.size());
01299   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
01300     rowList[i] = A().RowMatrixRowMap().GID(i);
01301   for (int i = 0 ; i < (int)ghostElements.size() ; ++i)
01302     rowList[i + NumMyRowsA_] = ghostElements[i];
01303 
01304   // row map for the overlapped matrix (local + overlap)
01305   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ghostElements.size(), &rowList[0], 0, Comm()) );
01306 
01307   // build the column map for the overlapping matrix
01308   //vector<int> colList(colMapElements.size());
01309   // column map for the overlapped matrix (local + overlap)
01310   //colMap_ = rcp( new Epetra_Map(-1, colMapElements.size(), &colList[0], 0, Comm()) );
01311   //for (int i = 0 ; i < (int)colMapElements.size() ; i++)
01312   //  colList[i] = colMapElements[i];
01313   vector<int> colList(A().RowMatrixColMap().NumMyElements() + colMapElements.size());
01314   int nc = A().RowMatrixColMap().NumMyElements();
01315   for (int i = 0 ; i < nc; i++)
01316     colList[i] = A().RowMatrixColMap().GID(i);
01317   for (int i = 0 ; i < (int)colMapElements.size() ; i++)
01318     colList[nc+i] = colMapElements[i];
01319 
01320   // column map for the overlapped matrix (local + overlap)
01321   //colMap_ = rcp( new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm()) );
01322   colMap_ = new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm());
01323     
01324 
01325   //FIXME timing
01326 #ifdef IFPACK_OVA_TIME_BUILD
01327   t1 = MPI_Wtime();
01328   fprintf(stderr,"[node %d]: time to init B() col/row maps %2.3e\n", nodeID, t1-t0);
01329   t0=t1;
01330 #endif
01331   //FIXME timing
01332 
01333   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01334   //++++ start of sort
01335   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01336   // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
01337   // different from that of Matrix.
01338   try {
01339     // build row map based on the local communicator.  We need this temporarily to build the column map.
01340     Epetra_Map* nodeMap = new Epetra_Map(-1,NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
01341     int numMyElts = colMap_->NumMyElements();
01342     assert(numMyElts!=0);
01343 
01344     // The column map *must* be sorted: first locals, then ghosts.  The ghosts
01345     // must be further sorted so that they are contiguous by owning processor.
01346 
01347     // For each GID in column map, determine owning PID in local communicator.
01348     int* myGlobalElts = new int[numMyElts];
01349     colMap_->MyGlobalElements(myGlobalElts);
01350     int* pidList = new int[numMyElts];
01351     nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
01352 
01353     // First sort column map GIDs based on the owning PID in local communicator.
01354     Epetra_Util Util;
01355     int *tt[1];
01356     tt[0] = myGlobalElts;
01357     Util.Sort(true, numMyElts, pidList, 0, (double**)0, 1, tt);
01358 
01359     // For each remote PID, sort the column map GIDs in ascending order.
01360     // Don't sort the local PID's entries.
01361     int localStart=0;
01362     while (localStart<numMyElts) {
01363       int currPID = (pidList)[localStart];
01364       int i=localStart;
01365       while (i<numMyElts && (pidList)[i] == currPID) i++;
01366       if (currPID != nodeComm->MyPID())
01367         Util.Sort(true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
01368       localStart = i;
01369     }
01370 
01371     // now move the local entries to the front of the list
01372     localStart=0;
01373     while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
01374     assert(localStart != numMyElts);
01375     int localEnd=localStart;
01376     while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
01377     int* mySortedGlobalElts = new int[numMyElts];
01378     int leng = localEnd - localStart;
01379     /* This is a little gotcha.  It appears that the ordering of the column map's local entries
01380        must be the same as that of the domain map's local entries.  See the comment in method
01381        MakeColMap() in Epetra_CrsGraph.cpp, line 1072. */
01382     int *rowGlobalElts =  nodeMap->MyGlobalElements();
01383 
01384     /*TODO TODO TODO NTS rows 68 and 83 show up as local GIDs in rowGlobalElts for both pids 0 & 1 on node 0.
01385                     This seems to imply that there is something wrong with rowList!!
01386                     It is ok for the overlapped matrix's row map to have repeated entries (it's overlapped, after all).
01387                     But ... on a node, there should be no repeats.
01388                     Ok, here's the underlying cause.  On node 0, gpid 1 finds 83 on overlap round 0.  gpid 0 finds 83
01389                     on overlap round 1.  This shouldn't happen .. once a node buddy "owns" a row, no other node buddy
01390                     should ever import ("own") that row.
01391 
01392                     Here's a possible fix.  Extend tie breaking across overlap rounds.  This means sending to lpid 0
01393                     a list of *all* rows you've imported (this round plus previous rounds) for tie-breaking
01394                     If a nb won a ghost row in a previous round, his votes for that ghost row should be really high, i.e.,
01395                     he should probably win that ghost row forever.
01396                     7/17/09
01397 
01398       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.
01399     */
01400 
01401     //move locals to front of list
01402     for (int i=0; i<leng; i++) {
01403       /* //original code */
01404       (mySortedGlobalElts)[i] = rowGlobalElts[i];
01405       //(&*mySortedGlobalElts)[i] = (&*myGlobalElts)[localStart+i];
01406       //(&*mySortedPidList)[i] = (&*pidList)[localStart+i];
01407     }
01408     for (int i=leng; i< localEnd; i++) {
01409       (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
01410     }
01411     for (int i=localEnd; i<numMyElts; i++) {
01412       (mySortedGlobalElts)[i] = (myGlobalElts)[i];
01413     }
01414 
01415     //FIXME timing
01416 #ifdef IFPACK_OVA_TIME_BUILD
01417     t1 = MPI_Wtime();
01418     fprintf(stderr,"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", nodeID, t1-t0);
01419     t0=t1;
01420 #endif
01421     //FIXME timing
01422 
01423     int indexBase = colMap_->IndexBase();
01424     if (colMap_) delete colMap_;
01425     colMap_ = new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,Comm());
01426 
01427     delete nodeMap;
01428     delete [] myGlobalElts;
01429     delete [] pidList;
01430     delete [] mySortedGlobalElts;
01431 
01432   } //try
01433   catch(...) {
01434     printf("** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
01435   }
01436 
01437   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01438   //++++ end of sort
01439   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01440 
01441   //FIXME timing
01442 #ifdef IFPACK_OVA_TIME_BUILD
01443   t1 = MPI_Wtime();
01444   fprintf(stderr,"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", nodeID, t1-t0);
01445   t0=t1;
01446 #endif
01447   //FIXME timing
01448 
01449 
01450 /*
01451 
01452    FIXME
01453    Does the column map need to be sorted for the overlapping matrix?
01454 
01455    The column map *must* be sorted:
01456 
01457         first locals
01458         then ghosts
01459 
01460    The ghosts must be further sorted so that they are contiguous by owning processor
01461 
01462   int* RemoteSizeList = 0
01463   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01464 
01465   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01466   Epetra_Util epetraUtil;
01467   SortLists[0] = RemoteColIndices;
01468   SortLists[1] = RemoteSizeList;
01469   epetraUtil.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
01470 */
01471 
01472   // now build the map corresponding to all the external nodes
01473   // (with respect to A().RowMatrixRowMap().
01474   ExtMap_ = rcp( new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,Comm()) );
01475   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*colMap_,0) ); 
01476 
01477   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
01478   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
01479 
01480   //FIXME timing
01481 #ifdef IFPACK_OVA_TIME_BUILD
01482   t1 = MPI_Wtime();
01483   fprintf(stderr,"[node %d]: time to import and setup ExtMap_ %2.3e\n", nodeID, t1-t0);
01484   t0=t1;
01485 #endif
01486   //FIXME timing
01487 
01488 
01489 /*
01490   Notes to self:    In FillComplete, the range map does not have to be 1-1 as long as
01491                     (row map == range map).  Ditto for the domain map not being 1-1
01492                     if (col map == domain map).
01493                     
01494 */
01495 
01496   ExtMatrix_->FillComplete( *colMap_ , *Map_ ); //FIXME wrong
01497 
01498   //FIXME timing
01499 #ifdef IFPACK_OVA_TIME_BUILD
01500   t1 = MPI_Wtime();
01501   fprintf(stderr,"[node %d]: time to FillComplete B() %2.3e\n", nodeID, t1-t0);
01502   t0=t1;
01503 #endif
01504   //FIXME timing
01505 
01506 
01507   // Note: B() = *ExtMatrix_ .
01508 
01509   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) ); //FIXME is this right?!
01510 
01511   // fix indices for overlapping matrix
01512   NumMyRowsB_ = B().NumMyRows();
01513   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;  //TODO is this wrong for a subdomain on >1 processor? // should be ok
01514   //NumMyCols_ = NumMyRows_;  //TODO is this wrong for a subdomain on >1 processor?  // YES!!!
01515   //NumMyCols_ = A().NumMyCols() + B().NumMyCols();
01516   NumMyCols_ = B().NumMyCols();
01517 
01518   /*FIXME*/ //somehow B's NumMyCols is the entire subdomain (local + overlap)
01519 
01520   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
01521   
01522   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
01523   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
01524   MaxNumEntries_ = A().MaxNumEntries();
01525   
01526   if (MaxNumEntries_ < B().MaxNumEntries())
01527     MaxNumEntries_ = B().MaxNumEntries();
01528 
01529 # ifdef HAVE_MPI
01530   delete nodeComm;
01531 # endif
01532 
01533   //FIXME timing
01534 #ifdef IFPACK_OVA_TIME_BUILD
01535   t1 = MPI_Wtime();
01536   fprintf(stderr,"[node %d]: time for final calcs %2.3e\n", nodeID, t1-t0);
01537   fprintf(stderr,"[node %d]: total IORM ctor time %2.3e\n", nodeID, t1-true_t0);
01538 #endif
01539   //FIXME timing
01540 
01541 
01542 } //Ifpack_OverlappingRowMatrix() ctor for more than one core
01543 
01544 // Destructor
01545 Ifpack_OverlappingRowMatrix::~Ifpack_OverlappingRowMatrix() {
01546   delete colMap_;
01547 }
01548 
01549 # endif //ifdef IFPACK_NODE_AWARE_CODE
01550 #endif // ifdef IFPACK_SUBCOMM_CODE
01551 
01552 
01553 // ======================================================================
01554 int Ifpack_OverlappingRowMatrix::
01555 NumMyRowEntries(int MyRow, int & NumEntries) const
01556 {
01557   if (MyRow < NumMyRowsA_)
01558     return(A().NumMyRowEntries(MyRow,NumEntries));
01559   else
01560     return(B().NumMyRowEntries(MyRow - NumMyRowsA_, NumEntries));
01561 }
01562 
01563 #ifdef IFPACK_SUBCOMM_CODE
01564 // ======================================================================
01565 int Ifpack_OverlappingRowMatrix::
01566 ExtractMyRowCopy(int LocRow, int Length, int & NumEntries, double *Values, 
01567                  int * Indices) const
01568 {
01569   //assert(1==0);
01570   int ierr;
01571   const Epetra_Map *Themap;
01572   if (LocRow < NumMyRowsA_) {
01573     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01574     Themap=&A().RowMatrixColMap();
01575   }
01576   else {
01577     ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01578     Themap=&B().RowMatrixColMap();
01579   }
01580 
01581   IFPACK_RETURN(ierr);
01582 }
01583 
01584 // ======================================================================
01585 int Ifpack_OverlappingRowMatrix::
01586 ExtractGlobalRowCopy(int GlobRow, int Length, int & NumEntries, double *Values, 
01587                  int * Indices) const
01588 {
01589   int ierr;
01590   const Epetra_Map *Themap;
01591   int LocRow = A().RowMatrixRowMap().LID(GlobRow);
01592   if (LocRow < NumMyRowsA_ && LocRow != -1) { //TODO don't need to check less than nummyrows
01593     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01594     Themap=&A().RowMatrixColMap();
01595   }
01596   else {
01597     LocRow = B().RowMatrixRowMap().LID(GlobRow);
01598     assert(LocRow!=-1);
01599     //ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01600     ierr = B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01601     Themap=&B().RowMatrixColMap();
01602   }
01603 
01604   for (int i=0; i<NumEntries; i++) {
01605     Indices[i]=Themap->GID(Indices[i]);
01606     assert(Indices[i] != -1);
01607   }
01608 
01609   IFPACK_RETURN(ierr);
01610 }
01611 #else
01612 # ifdef IFPACK_NODE_AWARE_CODE
01613 // ======================================================================
01614 int Ifpack_OverlappingRowMatrix::
01615 ExtractMyRowCopy(int LocRow, int Length, int & NumEntries, double *Values, 
01616                  int * Indices) const
01617 {
01618   assert(1==0);
01619   int ierr;
01620   const Epetra_Map *Themap;
01621   if (LocRow < NumMyRowsA_) {
01622     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01623     Themap=&A().RowMatrixColMap();
01624   }
01625   else {
01626     ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01627     Themap=&B().RowMatrixColMap();
01628   }
01629 
01630   IFPACK_RETURN(ierr);
01631 }
01632 
01633 // ======================================================================
01634 int Ifpack_OverlappingRowMatrix::
01635 ExtractGlobalRowCopy(int GlobRow, int Length, int & NumEntries, double *Values, 
01636                  int * Indices) const
01637 {
01638   int ierr;
01639   const Epetra_Map *Themap;
01640   int LocRow = A().RowMatrixRowMap().LID(GlobRow);
01641   if (LocRow < NumMyRowsA_ && LocRow != -1) { //TODO don't need to check less than nummyrows
01642     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01643     Themap=&A().RowMatrixColMap();
01644   }
01645   else {
01646     LocRow = B().RowMatrixRowMap().LID(GlobRow);
01647     assert(LocRow!=-1);
01648     //ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01649     ierr = B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01650     Themap=&B().RowMatrixColMap();
01651   }
01652 
01653   for (int i=0; i<NumEntries; i++) {
01654     Indices[i]=Themap->GID(Indices[i]);
01655     assert(Indices[i] != -1);
01656   }
01657 
01658   IFPACK_RETURN(ierr);
01659 }
01660 # else
01661 
01662 // ======================================================================
01663 int Ifpack_OverlappingRowMatrix::
01664 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values,
01665                  int * Indices) const
01666 {
01667   int ierr;
01668   if (MyRow < NumMyRowsA_)
01669     ierr = A().ExtractMyRowCopy(MyRow,Length,NumEntries,Values,Indices);
01670   else
01671     ierr = B().ExtractMyRowCopy(MyRow - NumMyRowsA_,Length,NumEntries,
01672                                 Values,Indices);
01673 
01674   IFPACK_RETURN(ierr);
01675 }
01676 # endif // ifdef IFPACK_NODE_AWARE_CODE
01677 #endif // ifdef IFPACK_SUBCOMM_CODE
01678 
01679 // ======================================================================
01680 int Ifpack_OverlappingRowMatrix::
01681 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
01682 {
01683   IFPACK_CHK_ERR(-1);
01684 }
01685 
01686 
01687 // ======================================================================
01688 int Ifpack_OverlappingRowMatrix::
01689 Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01690 {
01691   int NumVectors = X.NumVectors();
01692   vector<int> Ind(MaxNumEntries_);
01693   vector<double> Val(MaxNumEntries_);
01694 
01695   Y.PutScalar(0.0);
01696 
01697   // matvec with A (local rows)
01698   for (int i = 0 ; i < NumMyRowsA_ ; ++i) {
01699     for (int k = 0 ; k < NumVectors ; ++k) {
01700       int Nnz;
01701       IFPACK_CHK_ERR(A().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
01702                                           &Val[0], &Ind[0]));
01703       for (int j = 0 ; j < Nnz ; ++j) {
01704         Y[k][i] += Val[j] * X[k][Ind[j]];
01705       }
01706     }
01707   }
01708 
01709   // matvec with B (overlapping rows)
01710   for (int i = 0 ; i < NumMyRowsB_ ; ++i) {
01711     for (int k = 0 ; k < NumVectors ; ++k) {
01712       int Nnz;
01713       IFPACK_CHK_ERR(B().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
01714                                           &Val[0], &Ind[0]));
01715       for (int j = 0 ; j < Nnz ; ++j) {
01716         Y[k][i + NumMyRowsA_] += Val[j] * X[k][Ind[j]];
01717       }
01718     }
01719   }
01720   return(0);
01721 }
01722 
01723 // ======================================================================
01724 int Ifpack_OverlappingRowMatrix::
01725 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01726 {
01727   IFPACK_CHK_ERR(Multiply(UseTranspose(),X,Y));
01728   return(0);
01729 }
01730 
01731 // ======================================================================
01732 int Ifpack_OverlappingRowMatrix::
01733 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01734 {
01735   IFPACK_CHK_ERR(-1);
01736 }
01737 
01738 // ======================================================================
01739 #ifndef IFPACK_SUBCOMM_CODE
01740 # ifndef IFPACK_NODE_AWARE_CODE
01741 Epetra_RowMatrix& Ifpack_OverlappingRowMatrix::B() const
01742 {
01743   return(*ExtMatrix_);
01744 }
01745 # endif
01746 #endif
01747 // ======================================================================
01748 const Epetra_BlockMap& Ifpack_OverlappingRowMatrix::Map() const
01749 {
01750   return(*Map_);
01751 }
01752 
01753 // ======================================================================
01754 int Ifpack_OverlappingRowMatrix::
01755 ImportMultiVector(const Epetra_MultiVector& X, Epetra_MultiVector& OvX,
01756                   Epetra_CombineMode CM)
01757 {
01758   OvX.Import(X,*Importer_,CM);
01759   return(0);
01760 }
01761 
01762 // ======================================================================
01763 int Ifpack_OverlappingRowMatrix::
01764 ExportMultiVector(const Epetra_MultiVector& OvX, Epetra_MultiVector& X,
01765                   Epetra_CombineMode CM)
01766 {
01767   X.Export(OvX,*Importer_,CM);
01768   return(0);
01769 }
01770 
 All Classes Files Functions Variables Enumerations Friends