Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_OverlappingRowMatrix.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Ifpack_OverlappingRowMatrix.h"
00032 #include "Epetra_RowMatrix.h"
00033 #include "Epetra_Map.h"
00034 #include "Epetra_CrsMatrix.h"
00035 #include "Epetra_Comm.h"
00036 #include "Epetra_MultiVector.h"
00037 
00038 #ifdef IFPACK_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   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00850   MaxNumEntries_ = A().MaxNumEntries();
00851   
00852   if (MaxNumEntries_ < B().MaxNumEntries())
00853     MaxNumEntries_ = B().MaxNumEntries();
00854 
00855 # ifdef HAVE_MPI
00856   delete nodeComm;
00857 # endif
00858 
00859   //FIXME timing
00860 #ifdef IFPACK_OVA_TIME_BUILD
00861   t1 = MPI_Wtime();
00862   fprintf(stderr,"[node %d]: time for final calcs %2.3e\n", subdomainID, t1-t0);
00863   fprintf(stderr,"[node %d]: total IORM ctor time %2.3e\n", subdomainID, t1-true_t0);
00864 #endif
00865   //FIXME timing
00866 
00867 
00868 } //Ifpack_OverlappingRowMatrix() ctor for more than one core
00869 
00870 #else
00871 # ifdef IFPACK_NODE_AWARE_CODE
00872 
00873 // ====================================================================== 
00874 // Constructor for the case of two or more cores per subdomain
00875 Ifpack_OverlappingRowMatrix::
00876 Ifpack_OverlappingRowMatrix(const RCP<const Epetra_RowMatrix>& Matrix_in,
00877                             int OverlapLevel_in, int nodeID)  :
00878   Matrix_(Matrix_in),
00879   OverlapLevel_(OverlapLevel_in)
00880 {
00881 
00882   //FIXME timing
00883 #ifdef IFPACK_OVA_TIME_BUILD
00884   double t0 = MPI_Wtime();
00885   double t1, true_t0=t0;
00886 #endif
00887   //FIXME timing
00888 
00889   const int votesMax = INT_MAX;
00890 
00891   // should not be here if no overlap
00892   if (OverlapLevel_in == 0)
00893     IFPACK_CHK_ERRV(-1);
00894 
00895   // nothing to do as well with one process
00896   if (Comm().NumProc() == 1)
00897     IFPACK_CHK_ERRV(-1);
00898 
00899   // nodeID is the node (aka socket) number, and so is system dependent
00900   // nodeComm is the communicator for all the processes on a particular node
00901   // these processes will have the same nodeID.
00902 # ifdef HAVE_MPI
00903   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &Comm() );
00904   assert(pComm != NULL);
00905   MPI_Comm nodeMPIComm;
00906   MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm);
00907   Epetra_MpiComm *nodeComm = new Epetra_MpiComm(nodeMPIComm);
00908 # else
00909   Epetra_SerialComm *nodeComm =  dynamic_cast<Epetra_MpiComm*>( &(Matrix_in->RowMatrixRowMap().Comm()) );
00910 # endif
00911   
00912   NumMyRowsA_ = A().NumMyRows();
00913 
00914   RCP<Epetra_Map> TmpMap;
00915   RCP<Epetra_CrsMatrix> TmpMatrix; 
00916   RCP<Epetra_Import> TmpImporter;
00917 
00918   // importing rows corresponding to elements that are 
00919   // in ColMap, but not in RowMap 
00920   const Epetra_Map *RowMap; 
00921   const Epetra_Map *ColMap; 
00922   const Epetra_Map *DomainMap;
00923 
00924   // TODO Count #connections from nodes I own to each ghost node
00925 
00926 
00927   //FIXME timing
00928 #ifdef IFPACK_OVA_TIME_BUILD
00929   t1 = MPI_Wtime(); 
00930   fprintf(stderr,"[node %d]: time for initialization %2.3e\n", nodeID, t1-t0);
00931   t0=t1;
00932 #endif
00933   //FIXME timing
00934 
00935 #ifdef IFPACK_OVA_TIME_BUILD
00936   t1 = MPI_Wtime();
00937   fprintf(stderr,"[node %d]: overlap hash table allocs %2.3e\n", nodeID, t1-t0);
00938   t0=t1;
00939 #endif
00940 
00941   Teuchos::Hashtable<int,int> colMapTable(3 * A().RowMatrixColMap().NumMyElements() );
00942 
00943   // ghostTable holds off-node GIDs that are connected to on-node rows and can potentially be this PID's overlap
00944   // TODO hopefully 3 times the # column entries is a reasonable table size
00945   Teuchos::Hashtable<int,int> ghostTable(3 * A().RowMatrixColMap().NumMyElements() );
00946 
00947   /* ** ************************************************************************** ** */
00948   /* ** ********************** start of main overlap loop ************************ ** */
00949   /* ** ************************************************************************** ** */
00950   for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
00951   {
00952     // nbTable holds node buddy GIDs that are connected to current PID's rows, i.e., GIDs that should be in the overlapped
00953     // matrix's column map
00954 
00955     if (TmpMatrix != Teuchos::null) {
00956       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00957       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00958       DomainMap = &(TmpMatrix->OperatorDomainMap());
00959     }
00960     else {
00961       RowMap = &(A().RowMatrixRowMap()); 
00962       ColMap = &(A().RowMatrixColMap()); 
00963       DomainMap = &(A().OperatorDomainMap());
00964     }
00965 
00966 #ifdef IFPACK_OVA_TIME_BUILD
00967     t1 = MPI_Wtime();
00968     fprintf(stderr,"[node %d]: overlap pointer pulls %2.3e\n", nodeID, t1-t0);
00969     t0=t1;
00970 #endif
00971     
00972     // For each column ID, determine the owning node (as opposed to core)
00973     // ID of the corresponding row.
00974     Epetra_IntVector colIdList( *ColMap );
00975     Epetra_IntVector rowIdList(*DomainMap);
00976     rowIdList.PutValue(nodeID);  
00977     Teuchos::RCP<Epetra_Import> nodeIdImporter = rcp(new Epetra_Import( *ColMap, *DomainMap ));
00978 
00979 #ifdef IFPACK_OVA_TIME_BUILD
00980     t1 = MPI_Wtime();
00981     fprintf(stderr,"[node %d]: overlap intvector/importer allocs %2.3e\n", nodeID, t1-t0);
00982     t0=t1;
00983 #endif
00984     
00985     colIdList.Import(rowIdList,*nodeIdImporter,Insert);
00986     
00987     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00988     int count = 0; 
00989 
00990 #ifdef IFPACK_OVA_TIME_BUILD
00991     t1 = MPI_Wtime();
00992     fprintf(stderr,"[node %d]: overlap col/row ID imports %2.3e\n", nodeID, t1-t0);
00993     t0=t1;
00994 #endif
00995 
00996     
00997     // define the set of off-node rows that are in ColMap but not in RowMap
00998     // This naturally does not include off-core rows that are on the same node as me, i.e., node buddy rows.
00999     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
01000       int GID = ColMap->GID(i); 
01001       int myvotes=0;
01002       if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
01003       if ( colIdList[i] != nodeID && myvotes < votesMax) // don't include anybody found in a previous round
01004       {
01005         int votes;
01006         if (ghostTable.containsKey(GID)) {
01007           votes = ghostTable.get(GID);
01008           votes++;
01009           ghostTable.put(GID,votes);
01010         } else {
01011           ghostTable.put(GID,1);
01012         }
01013       }
01014     } //for (int i = 0 ; i < ColMap->NumMyElements() ; ++i)
01015 
01016     Teuchos::Array<int> gidsarray,votesarray;
01017     ghostTable.arrayify(gidsarray,votesarray);
01018     int *gids = gidsarray.getRawPtr();
01019     int *votes = votesarray.getRawPtr();
01020 
01021     /*
01022        This next bit of code decides which node buddy (NB) gets which ghost points.  Everyone sends their
01023        list of ghost points to pid 0 of the local subcommunicator.  Pid 0 decides who gets what:
01024 
01025           - if a ghost point is touched by only one NB, that NB gets the ghost point
01026           - if two or more NBs share a ghost point, whichever NB has the most connections to the ghost
01027             point gets it.
01028     */
01029 
01030 #   ifdef HAVE_MPI  //FIXME What if we build in serial?!  This file will likely not compile.
01031     int *cullList;
01032     int ncull;
01033     //int mypid = nodeComm->MyPID();
01034 
01035     //FIXME timing
01036 #ifdef IFPACK_OVA_TIME_BUILD
01037     t1 = MPI_Wtime();
01038     fprintf(stderr,"[node %d]: overlap before culling %2.3e\n", nodeID, t1-t0);
01039     t0=t1;
01040 #endif
01041     //FIXME timing
01042 
01043 
01044     if (nodeComm->MyPID() == 0)
01045     {
01046       // Figure out how much pid 0 is to receive
01047       MPI_Status status;
01048       int *lengths = new int[nodeComm->NumProc()+1];
01049       lengths[0] = 0;
01050       lengths[1] = ghostTable.size();
01051       for (int i=1; i<nodeComm->NumProc(); i++) {
01052         int leng;
01053         MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01054         lengths[i+1] = lengths[i] + leng;
01055       }
01056       int total = lengths[nodeComm->NumProc()];
01057 
01058       int* ghosts = new int[total];
01059       for (int i=0; i<total; i++) ghosts[i] = -9;
01060       int *round  = new int[total];
01061       int *owningpid  = new int[total];
01062 
01063       for (int i=1; i<nodeComm->NumProc(); i++) {
01064         int count = lengths[i+1] - lengths[i];
01065         MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01066         MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
01067       }
01068 
01069       // put in pid 0's info
01070       for (int i=0; i<lengths[1]; i++) {
01071         ghosts[i] = gids[i];
01072         round[i] = votes[i];
01073         owningpid[i] = 0;
01074       }
01075 
01076       // put in the pid associated with each ghost
01077       for (int j=1; j<nodeComm->NumProc(); j++) {
01078         for (int i=lengths[j]; i<lengths[j+1]; i++) {
01079           owningpid[i] = j;
01080         }
01081       }
01082 
01083       // sort everything based on the ghost gids
01084       int* roundpid[2];
01085       roundpid[0] = round; roundpid[1] = owningpid;
01086       Epetra_Util epetraUtil;
01087       epetraUtil.Sort(true,total,ghosts,0,0,2,roundpid);
01088 
01089       //set up arrays that get sent back to node buddies and that tell them what ghosts to cull
01090       int *nlosers = new int[nodeComm->NumProc()];
01091       int** losers = new int*[nodeComm->NumProc()];
01092       for (int i=0; i<nodeComm->NumProc(); i++) {
01093         nlosers[i] = 0;
01094         losers[i] = new int[ lengths[i+1]-lengths[i] ];
01095       }
01096 
01097       // Walk through ghosts array and and for each sequence of duplicate ghost GIDs, choose just one NB to keep it.
01098       // The logic is pretty simple.  The ghost list is sorted, so all duplicate PIDs are together.
01099       // The list is traversed.  As duplicates are found, node pid 0 keeps track of the current "winning"
01100       // pid.  When a pid is determined to have "lost" (less votes/connections to the current GID), the
01101       // GID is added to that pid's list of GIDs to be culled.  At the end of the repeated sequence, we have
01102       // a winner, and other NBs know whether they need to delete it from their import list.
01103       int max=0;   //for duplicated ghosts, index of pid with most votes and who hence keeps the ghost.
01104                    // TODO to break ties randomly
01105 
01106       for (int i=1; i<total; i++) {
01107         if (ghosts[i] == ghosts[i-1]) {
01108           int idx = i; // pid associated with idx is current "loser"
01109           if (round[i] > round[max]) {
01110             idx = max;
01111             max=i;
01112           }
01113           int j = owningpid[idx];
01114           losers[j][nlosers[j]++] = ghosts[idx];
01115         } else {
01116           max=i;
01117         }
01118       } //for (int i=1; i<total; i++)
01119 
01120       delete [] round;
01121       delete [] ghosts;
01122       delete [] owningpid;
01123 
01124       // send the arrays of ghost GIDs to be culled back to the respective node buddies
01125       for (int i=1; i<nodeComm->NumProc(); i++) {
01126         MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
01127         MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
01128       }
01129 
01130       //FIXME Unnecessary memory allocation and copying, but makes culling code cleaner
01131       //Could we stick this info into gids and votes, since neither is being used anymore?
01132       //TODO Instead of using "losers" arrays, just use "cullList" as in the else clause
01133       ncull = nlosers[0];
01134       cullList = new int[ncull+1];
01135       for (int i=0; i<nlosers[0]; i++)
01136         cullList[i] = losers[0][i];
01137 
01138       for (int i=0; i<nodeComm->NumProc(); i++)
01139         delete [] losers[i];
01140 
01141       delete [] lengths;
01142       delete [] losers;
01143       delete [] nlosers;
01144 
01145     } else { //everyone but pid 0
01146 
01147       // send to node pid 0 all ghosts that this pid could potentially import
01148       int hashsize = ghostTable.size();
01149       MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
01150       MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
01151       MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
01152 
01153       // receive the ghost GIDs that should not be imported (subset of the list sent off just above)
01154       MPI_Status status;
01155       MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
01156       cullList = new int[ncull+1];
01157       MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
01158     }
01159 
01160 
01161     //TODO clean out hash table after each time through overlap loop   4/1/07 JJH done moved both hash tables to local scope
01162 
01163     // Remove from my row map all off-node ghosts that will be imported by a node buddy.
01164     for (int i=0; i<ncull; i++) {
01165       try{ghostTable.remove(cullList[i]);}
01166 
01167       catch(...) {
01168         printf("pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
01169                Comm().MyPID(),cullList[i]);
01170         fflush(stdout);
01171       }
01172     }//for
01173 
01174     delete [] cullList;
01175 
01176     // Save off the remaining ghost GIDs from the current overlap round.
01177     // These are off-node GIDs (rows) that I will import.
01178     gidsarray.clear(); votesarray.clear();
01179     ghostTable.arrayify(gidsarray,votesarray);
01180 
01181     vector<int> list(size); 
01182     count=0;
01183     for (int i=0; i<ghostTable.size(); i++) {
01184       // if votesarray[i] == votesmax, then this GID was found during a previous overlap round
01185       if (votesarray[i] < votesMax) {
01186         list[count] = gidsarray[i];
01187         ghostTable.put(gidsarray[i],votesMax); //set this GID's vote to the maximum so that this PID alway wins in any subsequent overlap tiebreaking.
01188         count++;
01189       }
01190     }
01191 
01192     //FIXME timing
01193 #ifdef IFPACK_OVA_TIME_BUILD
01194     t1 = MPI_Wtime();
01195     fprintf(stderr,"[node %d]: overlap duplicate removal %2.3e\n", nodeID, t1-t0);
01196     t0=t1;
01197 #endif
01198     //FIXME timing
01199 
01200 #   endif //ifdef HAVE_MPI
01201 
01202     TmpMap = rcp( new Epetra_Map(-1,count, &list[0],0,Comm()) );
01203 
01204     TmpMatrix = rcp( new Epetra_CrsMatrix(Copy,*TmpMap,0) ); 
01205 
01206     TmpImporter = rcp( new Epetra_Import(*TmpMap,A().RowMatrixRowMap()) ); 
01207 
01208     TmpMatrix->Import(A(),*TmpImporter,Insert); 
01209     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
01210 
01211     //FIXME timing
01212 #ifdef IFPACK_OVA_TIME_BUILD
01213     t1 = MPI_Wtime();
01214     fprintf(stderr,"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", nodeID, t1-t0);
01215     t0=t1;
01216 #endif
01217     //FIXME timing
01218 
01219 
01220     // These next two imports get the GIDs that need to go into the column map of the overlapped matrix.
01221 
01222     // For each column ID in the overlap, determine the owning node (as opposed to core)
01223     // ID of the corresponding row.  Save those column IDs whose owning node is the current one.
01224     // This should get all the imported ghost GIDs.
01225     Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
01226     ov_colIdList.PutValue(-1);
01227     Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
01228     ov_rowIdList.PutValue(nodeID);  
01229     Teuchos::RCP<Epetra_Import> ov_nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), TmpMatrix->RowMap()));
01230     ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
01231 
01232     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
01233        if (ov_colIdList[i] == nodeID) {
01234          int GID = (ov_colIdList.Map()).GID(i);
01235          colMapTable.put(GID,1);
01236        }
01237     }
01238 
01239     // Do a second import of the owning node ID from A's rowmap to TmpMat's column map.  This ensures that
01240     // all GIDs that belong to a node buddy and are in a ghost row's sparsity pattern will be in the final
01241     // overlapped matrix's column map.
01242     ov_colIdList.PutValue(-1);
01243     Epetra_IntVector ArowIdList( A().RowMatrixRowMap() );
01244     ArowIdList.PutValue(nodeID);
01245     nodeIdImporter = rcp(new Epetra_Import( TmpMatrix->ColMap(), A().RowMatrixRowMap() ));
01246     ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
01247 
01248     for (int i=0 ; i<ov_colIdList.MyLength(); i++) {
01249        if (ov_colIdList[i] == nodeID) {
01250          int GID = (ov_colIdList.Map()).GID(i);
01251          colMapTable.put(GID,1);
01252        }
01253     }
01254 
01255     //FIXME timing
01256 #ifdef IFPACK_OVA_TIME_BUILD
01257     t1 = MPI_Wtime();
01258     fprintf(stderr,"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", nodeID, t1-t0);
01259     t0=t1;
01260 #endif
01261     //FIXME timing
01262 
01263   } //for (int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
01264 
01265   /* ** ************************************************************************ ** */
01266   /* ** ********************** end of main overlap loop ************************ ** */
01267   /* ** ************************************************************************ ** */
01268 
01269   // off-node GIDs that will be in the overlap
01270   vector<int> ghostElements; 
01271 
01272   Teuchos::Array<int> gidsarray,votesarray;
01273   ghostTable.arrayify(gidsarray,votesarray);
01274   for (int i=0; i<ghostTable.size(); i++) {
01275     ghostElements.push_back(gidsarray[i]);
01276   }
01277 
01278     for (int i = 0 ; i < A().RowMatrixColMap().NumMyElements() ; ++i) {
01279       int GID = A().RowMatrixColMap().GID(i);
01280       // Remove any entries that are in A's original column map
01281       if (colMapTable.containsKey(GID)) {
01282         try{colMapTable.remove(GID);}
01283         catch(...) {
01284           printf("pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n", Comm().MyPID(),GID);
01285           fflush(stdout);
01286         }
01287       }
01288     }
01289 
01290     // GIDs that will be in the overlapped matrix's column map
01291     vector<int> colMapElements; 
01292 
01293   gidsarray.clear(); votesarray.clear();
01294   colMapTable.arrayify(gidsarray,votesarray);
01295   for (int i=0; i<colMapTable.size(); i++)
01296     colMapElements.push_back(gidsarray[i]);
01297 
01298 /*
01299    We need two maps here.  The first is the row map, which we've got by using my original rows
01300    plus whatever I've picked up in ghostElements.
01301 
01302    The second is a column map.  This map should include my rows, plus any columns that could come from node buddies.
01303    These GIDs come from the std:array colMapElements, which in turn comes from colMapTable.
01304    This map should *not* omit ghosts that have been imported by my node buddies, i.e., for any row that I own,
01305    the stencil should include all column GIDs (including imported ghosts) that are on the node.
01306 */
01307 
01308   // build the row map containing all the nodes (original matrix + off-node matrix)
01309   vector<int> rowList(NumMyRowsA_ + ghostElements.size());
01310   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
01311     rowList[i] = A().RowMatrixRowMap().GID(i);
01312   for (int i = 0 ; i < (int)ghostElements.size() ; ++i)
01313     rowList[i + NumMyRowsA_] = ghostElements[i];
01314 
01315   // row map for the overlapped matrix (local + overlap)
01316   Map_ = rcp( new Epetra_Map(-1, NumMyRowsA_ + ghostElements.size(), &rowList[0], 0, Comm()) );
01317 
01318   // build the column map for the overlapping matrix
01319   //vector<int> colList(colMapElements.size());
01320   // column map for the overlapped matrix (local + overlap)
01321   //colMap_ = rcp( new Epetra_Map(-1, colMapElements.size(), &colList[0], 0, Comm()) );
01322   //for (int i = 0 ; i < (int)colMapElements.size() ; i++)
01323   //  colList[i] = colMapElements[i];
01324   vector<int> colList(A().RowMatrixColMap().NumMyElements() + colMapElements.size());
01325   int nc = A().RowMatrixColMap().NumMyElements();
01326   for (int i = 0 ; i < nc; i++)
01327     colList[i] = A().RowMatrixColMap().GID(i);
01328   for (int i = 0 ; i < (int)colMapElements.size() ; i++)
01329     colList[nc+i] = colMapElements[i];
01330 
01331   // column map for the overlapped matrix (local + overlap)
01332   //colMap_ = rcp( new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm()) );
01333   colMap_ = new Epetra_Map(-1, A().RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0, Comm());
01334     
01335 
01336   //FIXME timing
01337 #ifdef IFPACK_OVA_TIME_BUILD
01338   t1 = MPI_Wtime();
01339   fprintf(stderr,"[node %d]: time to init B() col/row maps %2.3e\n", nodeID, t1-t0);
01340   t0=t1;
01341 #endif
01342   //FIXME timing
01343 
01344   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01345   //++++ start of sort
01346   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01347   // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
01348   // different from that of Matrix.
01349   try {
01350     // build row map based on the local communicator.  We need this temporarily to build the column map.
01351     Epetra_Map* nodeMap = new Epetra_Map(-1,NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
01352     int numMyElts = colMap_->NumMyElements();
01353     assert(numMyElts!=0);
01354 
01355     // The column map *must* be sorted: first locals, then ghosts.  The ghosts
01356     // must be further sorted so that they are contiguous by owning processor.
01357 
01358     // For each GID in column map, determine owning PID in local communicator.
01359     int* myGlobalElts = new int[numMyElts];
01360     colMap_->MyGlobalElements(myGlobalElts);
01361     int* pidList = new int[numMyElts];
01362     nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
01363 
01364     // First sort column map GIDs based on the owning PID in local communicator.
01365     Epetra_Util Util;
01366     int *tt[1];
01367     tt[0] = myGlobalElts;
01368     Util.Sort(true, numMyElts, pidList, 0, (double**)0, 1, tt);
01369 
01370     // For each remote PID, sort the column map GIDs in ascending order.
01371     // Don't sort the local PID's entries.
01372     int localStart=0;
01373     while (localStart<numMyElts) {
01374       int currPID = (pidList)[localStart];
01375       int i=localStart;
01376       while (i<numMyElts && (pidList)[i] == currPID) i++;
01377       if (currPID != nodeComm->MyPID())
01378         Util.Sort(true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
01379       localStart = i;
01380     }
01381 
01382     // now move the local entries to the front of the list
01383     localStart=0;
01384     while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
01385     assert(localStart != numMyElts);
01386     int localEnd=localStart;
01387     while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
01388     int* mySortedGlobalElts = new int[numMyElts];
01389     int leng = localEnd - localStart;
01390     /* This is a little gotcha.  It appears that the ordering of the column map's local entries
01391        must be the same as that of the domain map's local entries.  See the comment in method
01392        MakeColMap() in Epetra_CrsGraph.cpp, line 1072. */
01393     int *rowGlobalElts =  nodeMap->MyGlobalElements();
01394 
01395     /*TODO TODO TODO NTS rows 68 and 83 show up as local GIDs in rowGlobalElts for both pids 0 & 1 on node 0.
01396                     This seems to imply that there is something wrong with rowList!!
01397                     It is ok for the overlapped matrix's row map to have repeated entries (it's overlapped, after all).
01398                     But ... on a node, there should be no repeats.
01399                     Ok, here's the underlying cause.  On node 0, gpid 1 finds 83 on overlap round 0.  gpid 0 finds 83
01400                     on overlap round 1.  This shouldn't happen .. once a node buddy "owns" a row, no other node buddy
01401                     should ever import ("own") that row.
01402 
01403                     Here's a possible fix.  Extend tie breaking across overlap rounds.  This means sending to lpid 0
01404                     a list of *all* rows you've imported (this round plus previous rounds) for tie-breaking
01405                     If a nb won a ghost row in a previous round, his votes for that ghost row should be really high, i.e.,
01406                     he should probably win that ghost row forever.
01407                     7/17/09
01408 
01409       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.
01410     */
01411 
01412     //move locals to front of list
01413     for (int i=0; i<leng; i++) {
01414       /* //original code */
01415       (mySortedGlobalElts)[i] = rowGlobalElts[i];
01416       //(&*mySortedGlobalElts)[i] = (&*myGlobalElts)[localStart+i];
01417       //(&*mySortedPidList)[i] = (&*pidList)[localStart+i];
01418     }
01419     for (int i=leng; i< localEnd; i++) {
01420       (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
01421     }
01422     for (int i=localEnd; i<numMyElts; i++) {
01423       (mySortedGlobalElts)[i] = (myGlobalElts)[i];
01424     }
01425 
01426     //FIXME timing
01427 #ifdef IFPACK_OVA_TIME_BUILD
01428     t1 = MPI_Wtime();
01429     fprintf(stderr,"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", nodeID, t1-t0);
01430     t0=t1;
01431 #endif
01432     //FIXME timing
01433 
01434     int indexBase = colMap_->IndexBase();
01435     if (colMap_) delete colMap_;
01436     colMap_ = new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,Comm());
01437 
01438     delete nodeMap;
01439     delete [] myGlobalElts;
01440     delete [] pidList;
01441     delete [] mySortedGlobalElts;
01442 
01443   } //try
01444   catch(...) {
01445     printf("** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
01446   }
01447 
01448   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01449   //++++ end of sort
01450   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01451 
01452   //FIXME timing
01453 #ifdef IFPACK_OVA_TIME_BUILD
01454   t1 = MPI_Wtime();
01455   fprintf(stderr,"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", nodeID, t1-t0);
01456   t0=t1;
01457 #endif
01458   //FIXME timing
01459 
01460 
01461 /*
01462 
01463    FIXME
01464    Does the column map need to be sorted for the overlapping matrix?
01465 
01466    The column map *must* be sorted:
01467 
01468         first locals
01469         then ghosts
01470 
01471    The ghosts must be further sorted so that they are contiguous by owning processor
01472 
01473   int* RemoteSizeList = 0
01474   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01475 
01476   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01477   Epetra_Util epetraUtil;
01478   SortLists[0] = RemoteColIndices;
01479   SortLists[1] = RemoteSizeList;
01480   epetraUtil.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
01481 */
01482 
01483   // now build the map corresponding to all the external nodes
01484   // (with respect to A().RowMatrixRowMap().
01485   ExtMap_ = rcp( new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,Comm()) );
01486   ExtMatrix_ = rcp( new Epetra_CrsMatrix(Copy,*ExtMap_,*colMap_,0) ); 
01487 
01488   ExtImporter_ = rcp( new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()) ); 
01489   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
01490 
01491   //FIXME timing
01492 #ifdef IFPACK_OVA_TIME_BUILD
01493   t1 = MPI_Wtime();
01494   fprintf(stderr,"[node %d]: time to import and setup ExtMap_ %2.3e\n", nodeID, t1-t0);
01495   t0=t1;
01496 #endif
01497   //FIXME timing
01498 
01499 
01500 /*
01501   Notes to self:    In FillComplete, the range map does not have to be 1-1 as long as
01502                     (row map == range map).  Ditto for the domain map not being 1-1
01503                     if (col map == domain map).
01504                     
01505 */
01506 
01507   ExtMatrix_->FillComplete( *colMap_ , *Map_ ); //FIXME wrong
01508 
01509   //FIXME timing
01510 #ifdef IFPACK_OVA_TIME_BUILD
01511   t1 = MPI_Wtime();
01512   fprintf(stderr,"[node %d]: time to FillComplete B() %2.3e\n", nodeID, t1-t0);
01513   t0=t1;
01514 #endif
01515   //FIXME timing
01516 
01517 
01518   // Note: B() = *ExtMatrix_ .
01519 
01520   Importer_ = rcp( new Epetra_Import(*Map_,A().RowMatrixRowMap()) ); //FIXME is this right?!
01521 
01522   // fix indices for overlapping matrix
01523   NumMyRowsB_ = B().NumMyRows();
01524   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;  //TODO is this wrong for a subdomain on >1 processor? // should be ok
01525   //NumMyCols_ = NumMyRows_;  //TODO is this wrong for a subdomain on >1 processor?  // YES!!!
01526   //NumMyCols_ = A().NumMyCols() + B().NumMyCols();
01527   NumMyCols_ = B().NumMyCols();
01528 
01529   /*FIXME*/ //somehow B's NumMyCols is the entire subdomain (local + overlap)
01530 
01531   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
01532   
01533   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
01534   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
01535   MaxNumEntries_ = A().MaxNumEntries();
01536   
01537   if (MaxNumEntries_ < B().MaxNumEntries())
01538     MaxNumEntries_ = B().MaxNumEntries();
01539 
01540 # ifdef HAVE_MPI
01541   delete nodeComm;
01542 # endif
01543 
01544   //FIXME timing
01545 #ifdef IFPACK_OVA_TIME_BUILD
01546   t1 = MPI_Wtime();
01547   fprintf(stderr,"[node %d]: time for final calcs %2.3e\n", nodeID, t1-t0);
01548   fprintf(stderr,"[node %d]: total IORM ctor time %2.3e\n", nodeID, t1-true_t0);
01549 #endif
01550   //FIXME timing
01551 
01552 
01553 } //Ifpack_OverlappingRowMatrix() ctor for more than one core
01554 
01555 // Destructor
01556 Ifpack_OverlappingRowMatrix::~Ifpack_OverlappingRowMatrix() {
01557   delete colMap_;
01558 }
01559 
01560 # endif //ifdef IFPACK_NODE_AWARE_CODE
01561 #endif // ifdef IFPACK_SUBCOMM_CODE
01562 
01563 
01564 // ======================================================================
01565 int Ifpack_OverlappingRowMatrix::
01566 NumMyRowEntries(int MyRow, int & NumEntries) const
01567 {
01568   if (MyRow < NumMyRowsA_)
01569     return(A().NumMyRowEntries(MyRow,NumEntries));
01570   else
01571     return(B().NumMyRowEntries(MyRow - NumMyRowsA_, NumEntries));
01572 }
01573 
01574 #ifdef IFPACK_SUBCOMM_CODE
01575 // ======================================================================
01576 int Ifpack_OverlappingRowMatrix::
01577 ExtractMyRowCopy(int LocRow, int Length, int & NumEntries, double *Values, 
01578                  int * Indices) const
01579 {
01580   //assert(1==0);
01581   int ierr;
01582   const Epetra_Map *Themap;
01583   if (LocRow < NumMyRowsA_) {
01584     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01585     Themap=&A().RowMatrixColMap();
01586   }
01587   else {
01588     ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01589     Themap=&B().RowMatrixColMap();
01590   }
01591 
01592   IFPACK_RETURN(ierr);
01593 }
01594 
01595 // ======================================================================
01596 int Ifpack_OverlappingRowMatrix::
01597 ExtractGlobalRowCopy(int GlobRow, int Length, int & NumEntries, double *Values, 
01598                  int * Indices) const
01599 {
01600   int ierr;
01601   const Epetra_Map *Themap;
01602   int LocRow = A().RowMatrixRowMap().LID(GlobRow);
01603   if (LocRow < NumMyRowsA_ && LocRow != -1) { //TODO don't need to check less than nummyrows
01604     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01605     Themap=&A().RowMatrixColMap();
01606   }
01607   else {
01608     LocRow = B().RowMatrixRowMap().LID(GlobRow);
01609     assert(LocRow!=-1);
01610     //ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01611     ierr = B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01612     Themap=&B().RowMatrixColMap();
01613   }
01614 
01615   for (int i=0; i<NumEntries; i++) {
01616     Indices[i]=Themap->GID(Indices[i]);
01617     assert(Indices[i] != -1);
01618   }
01619 
01620   IFPACK_RETURN(ierr);
01621 }
01622 #else
01623 # ifdef IFPACK_NODE_AWARE_CODE
01624 // ======================================================================
01625 int Ifpack_OverlappingRowMatrix::
01626 ExtractMyRowCopy(int LocRow, int Length, int & NumEntries, double *Values, 
01627                  int * Indices) const
01628 {
01629   assert(1==0);
01630   int ierr;
01631   const Epetra_Map *Themap;
01632   if (LocRow < NumMyRowsA_) {
01633     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01634     Themap=&A().RowMatrixColMap();
01635   }
01636   else {
01637     ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01638     Themap=&B().RowMatrixColMap();
01639   }
01640 
01641   IFPACK_RETURN(ierr);
01642 }
01643 
01644 // ======================================================================
01645 int Ifpack_OverlappingRowMatrix::
01646 ExtractGlobalRowCopy(int GlobRow, int Length, int & NumEntries, double *Values, 
01647                  int * Indices) const
01648 {
01649   int ierr;
01650   const Epetra_Map *Themap;
01651   int LocRow = A().RowMatrixRowMap().LID(GlobRow);
01652   if (LocRow < NumMyRowsA_ && LocRow != -1) { //TODO don't need to check less than nummyrows
01653     ierr = A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01654     Themap=&A().RowMatrixColMap();
01655   }
01656   else {
01657     LocRow = B().RowMatrixRowMap().LID(GlobRow);
01658     assert(LocRow!=-1);
01659     //ierr = B().ExtractMyRowCopy(LocRow-NumMyRowsA_,Length,NumEntries,Values,Indices);
01660     ierr = B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
01661     Themap=&B().RowMatrixColMap();
01662   }
01663 
01664   for (int i=0; i<NumEntries; i++) {
01665     Indices[i]=Themap->GID(Indices[i]);
01666     assert(Indices[i] != -1);
01667   }
01668 
01669   IFPACK_RETURN(ierr);
01670 }
01671 # else
01672 
01673 // ======================================================================
01674 int Ifpack_OverlappingRowMatrix::
01675 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values,
01676                  int * Indices) const
01677 {
01678   int ierr;
01679   if (MyRow < NumMyRowsA_)
01680     ierr = A().ExtractMyRowCopy(MyRow,Length,NumEntries,Values,Indices);
01681   else
01682     ierr = B().ExtractMyRowCopy(MyRow - NumMyRowsA_,Length,NumEntries,
01683                                 Values,Indices);
01684 
01685   IFPACK_RETURN(ierr);
01686 }
01687 # endif // ifdef IFPACK_NODE_AWARE_CODE
01688 #endif // ifdef IFPACK_SUBCOMM_CODE
01689 
01690 // ======================================================================
01691 int Ifpack_OverlappingRowMatrix::
01692 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
01693 {
01694   IFPACK_CHK_ERR(-1);
01695 }
01696 
01697 
01698 // ======================================================================
01699 int Ifpack_OverlappingRowMatrix::
01700 Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01701 {
01702   int NumVectors = X.NumVectors();
01703   vector<int> Ind(MaxNumEntries_);
01704   vector<double> Val(MaxNumEntries_);
01705 
01706   Y.PutScalar(0.0);
01707 
01708   // matvec with A (local rows)
01709   for (int i = 0 ; i < NumMyRowsA_ ; ++i) {
01710     for (int k = 0 ; k < NumVectors ; ++k) {
01711       int Nnz;
01712       IFPACK_CHK_ERR(A().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
01713                                           &Val[0], &Ind[0]));
01714       for (int j = 0 ; j < Nnz ; ++j) {
01715         Y[k][i] += Val[j] * X[k][Ind[j]];
01716       }
01717     }
01718   }
01719 
01720   // matvec with B (overlapping rows)
01721   for (int i = 0 ; i < NumMyRowsB_ ; ++i) {
01722     for (int k = 0 ; k < NumVectors ; ++k) {
01723       int Nnz;
01724       IFPACK_CHK_ERR(B().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
01725                                           &Val[0], &Ind[0]));
01726       for (int j = 0 ; j < Nnz ; ++j) {
01727         Y[k][i + NumMyRowsA_] += Val[j] * X[k][Ind[j]];
01728       }
01729     }
01730   }
01731   return(0);
01732 }
01733 
01734 // ======================================================================
01735 int Ifpack_OverlappingRowMatrix::
01736 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01737 {
01738   IFPACK_CHK_ERR(Multiply(UseTranspose(),X,Y));
01739   return(0);
01740 }
01741 
01742 // ======================================================================
01743 int Ifpack_OverlappingRowMatrix::
01744 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01745 {
01746   IFPACK_CHK_ERR(-1);
01747 }
01748 
01749 // ======================================================================
01750 #ifndef IFPACK_SUBCOMM_CODE
01751 # ifndef IFPACK_NODE_AWARE_CODE
01752 Epetra_RowMatrix& Ifpack_OverlappingRowMatrix::B() const
01753 {
01754   return(*ExtMatrix_);
01755 }
01756 # endif
01757 #endif
01758 // ======================================================================
01759 const Epetra_BlockMap& Ifpack_OverlappingRowMatrix::Map() const
01760 {
01761   return(*Map_);
01762 }
01763 
01764 // ======================================================================
01765 int Ifpack_OverlappingRowMatrix::
01766 ImportMultiVector(const Epetra_MultiVector& X, Epetra_MultiVector& OvX,
01767                   Epetra_CombineMode CM)
01768 {
01769   OvX.Import(X,*Importer_,CM);
01770   return(0);
01771 }
01772 
01773 // ======================================================================
01774 int Ifpack_OverlappingRowMatrix::
01775 ExportMultiVector(const Epetra_MultiVector& OvX, Epetra_MultiVector& X,
01776                   Epetra_CombineMode CM)
01777 {
01778   X.Export(OvX,*Importer_,CM);
01779   return(0);
01780 }
01781 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines