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