EpetraExt Development
EpetraExt_MMHelpers.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_ConfigDefs.h>
00043 #include <EpetraExt_MMHelpers.h>
00044 #include <Epetra_Comm.h>
00045 #include <Epetra_CrsMatrix.h>
00046 #include <Epetra_Import.h>
00047 #include <Epetra_Distributor.h>
00048 #include <Epetra_HashTable.h>
00049 #include <Epetra_Util.h>
00050 #include <Epetra_Import_Util.h>
00051 #include <Epetra_GIDTypeSerialDenseVector.h>
00052 
00053 #include <Teuchos_TimeMonitor.hpp>
00054 #include <limits>
00055 
00056 #ifdef HAVE_MPI
00057 #include "Epetra_MpiComm.h"
00058 #include "Epetra_MpiDistributor.h"
00059 #endif
00060 #define MIN(x,y)    ((x)<(y)?(x):(y))
00061 #define MIN3(x,y,z) ((x)<(y)?(MIN(x,z)):(MIN(y,z)))
00062 
00063 namespace EpetraExt {
00064 
00065 //------------------------------------
00066 // DEBUGGING ROUTINES
00067 void debug_print_distor(const char * label, const Epetra_Distributor * Distor, const Epetra_Comm & Comm) {
00068 #ifdef HAVE_MPI
00069   const Epetra_MpiDistributor * MDistor = dynamic_cast<const Epetra_MpiDistributor*>(Distor);
00070   printf("[%d] %s\n",Comm.MyPID(),label);
00071   printf("[%d] NumSends = %d NumRecvs = %d\n",Comm.MyPID(),MDistor->NumSends(),MDistor->NumReceives());
00072   printf("[%d] ProcsTo = ",Comm.MyPID());
00073   for(int ii=0; ii<MDistor->NumSends(); ii++)
00074     printf("%d ",MDistor->ProcsTo()[ii]);
00075   printf("\n[%d] ProcsFrom = ",Comm.MyPID());
00076   for(int ii=0; ii<MDistor->NumReceives(); ii++)
00077     printf("%d ",MDistor->ProcsFrom()[ii]);
00078   printf("\n");
00079   fflush(stdout);
00080 #endif
00081 }
00082 
00083 //------------------------------------
00084 // DEBUGGING ROUTINES
00085 void debug_compare_import(const Epetra_Import * Import1,const Epetra_Import * Import2) {
00086   if(!Import1 && !Import2) return;
00087   const Epetra_Comm & Comm = (Import1)? Import1->SourceMap().Comm() : Import2->SourceMap().Comm();
00088   bool flag=true;
00089   int PID=Comm.MyPID();
00090   if( (!Import1 && Import2) || (Import2 && !Import1) ) {printf("[%d] DCI: One Import exists, the other does not\n",PID);return;}
00091   if(!Import1->SourceMap().SameAs(Import2->SourceMap())) {printf("[%d] DCI: SourceMaps don't match\n",PID);return;}
00092   if(!Import1->TargetMap().SameAs(Import2->TargetMap())) {printf("[%d] DCI: TargetMaps don't match\n",PID);return;}
00093 
00094   if(Import1->NumSameIDs() != Import2->NumSameIDs()) {printf("[%d] DCI NumSameIDs() mismatch %d vs. %d\n",PID,Import1->NumSameIDs(),Import2->NumSameIDs());flag=false;}
00095 
00096   if(Import1->NumPermuteIDs() != Import2->NumPermuteIDs()) {printf("[%d] DCI NumPermuteIDs() mismatch %d vs. %d\n",PID,Import1->NumPermuteIDs(),Import2->NumPermuteIDs()); flag=false;}
00097 
00098   if(Import1->NumExportIDs() != Import2->NumExportIDs()) {printf("[%d] DCI NumExportIDs() mismatch %d vs. %d\n",PID,Import1->NumExportIDs(),Import2->NumExportIDs()); flag=false;}
00099 
00100   if(Import1->NumRemoteIDs() != Import2->NumRemoteIDs()) {printf("[%d] DCI NumRemoteIDs() mismatch %d vs. %d\n",PID,Import1->NumRemoteIDs(),Import2->NumRemoteIDs()); flag=false;}
00101 
00102   if(Import1->NumSend() != Import2->NumSend()) {printf("[%d] DCI NumSend() mismatch %d vs. %d\n",PID,Import1->NumSend(),Import2->NumSend()); flag=false;}
00103 
00104   if(Import1->NumRecv() != Import2->NumRecv()) {printf("[%d] DCI NumRecv() mismatch %d vs. %d\n",PID,Import1->NumRecv(),Import2->NumRecv()); flag=false;}
00105 
00106 
00107   if(flag) printf("[%d] DCI Importers compare OK\n",PID);    
00108   fflush(stdout);
00109   Import1->SourceMap().Comm().Barrier();
00110   Import1->SourceMap().Comm().Barrier();
00111   Import1->SourceMap().Comm().Barrier();
00112   if(!flag) exit(1);
00113 }
00114 
00115 
00116 
00117 //------------------------------------
00118 CrsMatrixStruct::CrsMatrixStruct()
00119  : numRows(0), numEntriesPerRow(NULL), indices(NULL), values(NULL),
00120    remote(NULL), numRemote(0), importColMap(NULL), rowMap(NULL), colMap(NULL),
00121    domainMap(NULL), importMatrix(NULL), origMatrix(NULL)
00122 {
00123 }
00124 
00125 CrsMatrixStruct::~CrsMatrixStruct()
00126 {
00127   deleteContents();
00128 }
00129 
00130 void CrsMatrixStruct::deleteContents()
00131 {
00132   numRows = 0;
00133   delete [] numEntriesPerRow; numEntriesPerRow = NULL;
00134   delete [] indices; indices = NULL;
00135   delete [] values; values = NULL;
00136   delete [] remote; remote = NULL;
00137   numRemote = 0;
00138   delete importMatrix; importMatrix=0;
00139   // origMatrix is not owned by me, so don't delete
00140   origMatrix=0;
00141   targetMapToOrigRow.resize(0);
00142   targetMapToImportRow.resize(0);
00143 }
00144 
00145 int dumpCrsMatrixStruct(const CrsMatrixStruct& M)
00146 {
00147   std::cout << "proc " << M.rowMap->Comm().MyPID()<<std::endl;
00148   std::cout << "numRows: " << M.numRows<<std::endl;
00149   for(int i=0; i<M.numRows; ++i) {
00150     for(int j=0; j<M.numEntriesPerRow[i]; ++j) {
00151       if (M.remote[i]) {
00152         std::cout << "  *"<<M.rowMap->GID64(i)<<"   "
00153              <<M.importColMap->GID64(M.indices[i][j])<<"   "<<M.values[i][j]<<std::endl;
00154       }
00155       else {
00156         std::cout << "   "<<M.rowMap->GID64(i)<<"   "
00157              <<M.colMap->GID64(M.indices[i][j])<<"   "<<M.values[i][j]<<std::endl;
00158       }
00159     }
00160   }
00161   return(0);
00162 }
00163 
00164 CrsWrapper_Epetra_CrsMatrix::CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix& epetracrsmatrix)
00165  : ecrsmat_(epetracrsmatrix)
00166 {
00167 }
00168 
00169 CrsWrapper_Epetra_CrsMatrix::~CrsWrapper_Epetra_CrsMatrix()
00170 {
00171 }
00172 
00173 const Epetra_Map&
00174 CrsWrapper_Epetra_CrsMatrix::RowMap() const
00175 {
00176   return ecrsmat_.RowMap();
00177 }
00178 
00179 bool CrsWrapper_Epetra_CrsMatrix::Filled()
00180 {
00181   return ecrsmat_.Filled();
00182 }
00183 
00184 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00185 int
00186 CrsWrapper_Epetra_CrsMatrix::InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00187 {
00188   return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
00189 }
00190 
00191 int
00192 CrsWrapper_Epetra_CrsMatrix::SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
00193 {
00194   return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
00195 }
00196 #endif
00197 
00198 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00199 int
00200 CrsWrapper_Epetra_CrsMatrix::InsertGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices)
00201 {
00202   return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
00203 }
00204 
00205 int
00206 CrsWrapper_Epetra_CrsMatrix::SumIntoGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices)
00207 {
00208   return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
00209 }
00210 #endif
00211 
00212 
00213 //------------------------------------
00214 
00215 template<typename int_type>
00216 CrsWrapper_GraphBuilder<int_type>::CrsWrapper_GraphBuilder(const Epetra_Map& emap)
00217  : graph_(),
00218    rowmap_(emap),
00219    max_row_length_(0)
00220 {
00221   int num_rows = emap.NumMyElements();
00222   int_type* rows = 0;
00223   emap.MyGlobalElementsPtr(rows);
00224 
00225   for(int i=0; i<num_rows; ++i) {
00226     graph_[rows[i]] = new std::set<int_type>;
00227   }
00228 }
00229 
00230 template<typename int_type>
00231 CrsWrapper_GraphBuilder<int_type>::~CrsWrapper_GraphBuilder()
00232 {
00233   typename std::map<int_type,std::set<int_type>*>::iterator
00234     iter = graph_.begin(), iter_end = graph_.end();
00235   for(; iter!=iter_end; ++iter) {
00236     delete iter->second;
00237   }
00238 
00239   graph_.clear();
00240 }
00241 
00242 template<typename int_type>
00243 bool CrsWrapper_GraphBuilder<int_type>::Filled()
00244 {
00245   return false;
00246 }
00247 
00248 template<typename int_type>
00249 int
00250 CrsWrapper_GraphBuilder<int_type>::InsertGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices)
00251 {
00252   typename std::map<int_type,std::set<int_type>*>::iterator
00253     iter = graph_.find(GlobalRow);
00254 
00255   if (iter == graph_.end()) return(-1);
00256 
00257   std::set<int_type>& cols = *(iter->second);
00258 
00259   for(int i=0; i<NumEntries; ++i) {
00260     cols.insert(Indices[i]);
00261   }
00262 
00263   int row_length = cols.size();
00264   if (row_length > max_row_length_) max_row_length_ = row_length;
00265 
00266   return(0);
00267 }
00268 
00269 template<typename int_type>
00270 int
00271 CrsWrapper_GraphBuilder<int_type>::SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices)
00272 {
00273   return InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
00274 }
00275 
00276 template<typename int_type>
00277 std::map<int_type,std::set<int_type>*>&
00278 CrsWrapper_GraphBuilder<int_type>::get_graph()
00279 {
00280   return graph_;
00281 }
00282 
00283 template<typename int_type>
00284 void insert_matrix_locations(CrsWrapper_GraphBuilder<int_type>& graphbuilder,
00285                               Epetra_CrsMatrix& C)
00286 {
00287   int max_row_length = graphbuilder.get_max_row_length();
00288   if (max_row_length < 1) return;
00289 
00290   std::vector<int_type> indices(max_row_length);
00291   int_type* indices_ptr = &indices[0];
00292   std::vector<double> zeros(max_row_length, 0.0);
00293   double* zeros_ptr = &zeros[0];
00294 
00295   std::map<int_type,std::set<int_type>*>& graph = graphbuilder.get_graph();
00296 
00297   typename std::map<int_type,std::set<int_type>*>::iterator
00298     iter = graph.begin(), iter_end = graph.end();
00299 
00300   for(; iter!=iter_end; ++iter) {
00301     int_type row = iter->first;
00302     std::set<int_type>& cols = *(iter->second);
00303     int num_entries = cols.size();
00304 
00305     typename std::set<int_type>::iterator
00306       col_iter = cols.begin(), col_end = cols.end();
00307     for(int j=0; col_iter!=col_end; ++col_iter, ++j) {
00308       indices_ptr[j] = *col_iter;
00309     }
00310 
00311     C.InsertGlobalValues(row, num_entries, zeros_ptr, indices_ptr);
00312   }
00313 }
00314 
00315 template<typename int_type>
00316 void Tpack_outgoing_rows(const Epetra_CrsMatrix& mtx,
00317                         const std::vector<int_type>& proc_col_ranges,
00318                         std::vector<int_type>& send_rows,
00319                         std::vector<int>& rows_per_send_proc)
00320 {
00321   const Epetra_Map& rowmap = mtx.RowMap();
00322   int numrows = mtx.NumMyRows();
00323   const Epetra_CrsGraph& graph = mtx.Graph();
00324   int rowlen = 0;
00325   int* col_indices = NULL;
00326   int_type* Tcol_indices = NULL;
00327   int num_col_ranges = proc_col_ranges.size()/2;
00328   rows_per_send_proc.resize(num_col_ranges);
00329   send_rows.clear();
00330   for(int nc=0; nc<num_col_ranges; ++nc) {
00331     int_type first_col = proc_col_ranges[nc*2];
00332     int_type last_col = proc_col_ranges[nc*2+1];
00333     int num_send_rows = 0;
00334     for(int i=0; i<numrows; ++i) {
00335       int_type grow = (int_type) rowmap.GID64(i);
00336       if (mtx.Filled()) {
00337         const Epetra_Map& colmap = mtx.ColMap();
00338         graph.ExtractMyRowView(i, rowlen, col_indices);
00339         for(int j=0; j<rowlen; ++j) {
00340           int_type col = (int_type) colmap.GID64(col_indices[j]);
00341           if (first_col <= col && last_col >= col) {
00342             ++num_send_rows;
00343             send_rows.push_back(grow);
00344             break;
00345           }
00346         }
00347       }
00348       else {
00349         graph.ExtractGlobalRowView(grow, rowlen, Tcol_indices);
00350         for(int j=0; j<rowlen; ++j) {
00351           if (first_col <= Tcol_indices[j] && last_col >= Tcol_indices[j]) {
00352             ++num_send_rows;
00353             send_rows.push_back(grow);
00354             break;
00355           }
00356         }
00357       }
00358     }
00359     rows_per_send_proc[nc] = num_send_rows;
00360   }
00361 }
00362 
00363 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00364 void pack_outgoing_rows(const Epetra_CrsMatrix& mtx,
00365                         const std::vector<int>& proc_col_ranges,
00366                         std::vector<int>& send_rows,
00367                         std::vector<int>& rows_per_send_proc)
00368 {
00369   if(mtx.RowMatrixRowMap().GlobalIndicesInt()) {
00370     Tpack_outgoing_rows<int>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
00371   }
00372   else {
00373     throw "EpetraExt::pack_outgoing_rows: Global indices not int";
00374   }
00375 }
00376 #endif
00377 
00378 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00379 void pack_outgoing_rows(const Epetra_CrsMatrix& mtx,
00380                         const std::vector<long long>& proc_col_ranges,
00381                         std::vector<long long>& send_rows,
00382                         std::vector<int>& rows_per_send_proc)
00383 {
00384   if(mtx.RowMatrixRowMap().GlobalIndicesLongLong()) {
00385     Tpack_outgoing_rows<long long>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
00386   }
00387   else {
00388     throw "EpetraExt::pack_outgoing_rows: Global indices not long long";
00389   }
00390 }
00391 #endif
00392 
00393 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00394 template<>
00395 std::pair<int,int> get_col_range<int>(const Epetra_Map& emap)
00396  {
00397   if(emap.GlobalIndicesInt())
00398     return std::make_pair(emap.MinMyGID(),emap.MaxMyGID());
00399   else
00400     throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_Map";
00401 }
00402 #endif
00403 
00404 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00405 template<>
00406 std::pair<long long,long long> get_col_range<long long>(const Epetra_Map& emap)
00407 {
00408   if(emap.GlobalIndicesLongLong())
00409     return std::make_pair(emap.MinMyGID64(),emap.MaxMyGID64());
00410   else
00411     throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_Map";
00412 }
00413 #endif
00414 
00415 template<typename int_type>
00416 std::pair<int_type,int_type> Tget_col_range(const Epetra_CrsMatrix& mtx)
00417 {
00418   std::pair<int_type,int_type> col_range;
00419   if (mtx.Filled()) {
00420     col_range = get_col_range<int_type>(mtx.ColMap());
00421   }
00422   else {
00423     const Epetra_Map& row_map = mtx.RowMap();
00424     col_range.first = row_map.MaxMyGID64();
00425     col_range.second = row_map.MinMyGID64();
00426     int rowlen = 0;
00427     int_type* col_indices = NULL;
00428     const Epetra_CrsGraph& graph = mtx.Graph();
00429     for(int i=0; i<row_map.NumMyElements(); ++i) {
00430       graph.ExtractGlobalRowView((int_type) row_map.GID64(i), rowlen, col_indices);
00431       for(int j=0; j<rowlen; ++j) {
00432         if (col_indices[j] < col_range.first) col_range.first = col_indices[j];
00433         if (col_indices[j] > col_range.second) col_range.second = col_indices[j];
00434       }
00435     }
00436   }
00437 
00438   return col_range;
00439 }
00440 
00441 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00442 template<>
00443 std::pair<int,int> get_col_range<int>(const Epetra_CrsMatrix& mtx)
00444 {
00445   if(mtx.RowMatrixRowMap().GlobalIndicesInt())
00446     return Tget_col_range<int>(mtx);
00447   else
00448     throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_CrsMatrix";
00449 }
00450 #endif
00451 
00452 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00453 template<>
00454 std::pair<long long,long long> get_col_range<long long>(const Epetra_CrsMatrix& mtx)
00455 {
00456   if(mtx.RowMatrixRowMap().GlobalIndicesLongLong())
00457     return Tget_col_range<long long>(mtx);
00458   else
00459     throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_CrsMatrix";
00460 }
00461 #endif
00462 
00463 
00464 /**********************************************************************************************/
00465 /**********************************************************************************************/
00466 /**********************************************************************************************/
00467 #ifdef HAVE_MPI
00468 template <typename MyType>
00469 void boundary_exchange(const Epetra_MpiComm Comm, MPI_Datatype DataType, 
00470            int NumSends, const int * SendProcs, const int * SendSizes, MyType* SendBuffer, 
00471            int NumRecvs, const int * RecvProcs, const int * RecvSizes, MyType* RecvBuffer,int SizeOfPacket,int msg_tag)
00472 {  
00473 
00474   MPI_Comm comm = Comm.Comm();
00475   std::vector<MPI_Request> requests(NumRecvs);
00476   std::vector<MPI_Status>  status(NumRecvs);
00477 
00478   int i,num_waits=0,MyPID=Comm.MyPID();
00479   int start, self_recv_len=-1,self_recv_start=-1, self_send_start=-1;
00480   
00481   // Default send/recv size if the Sizes arrays are NULL.
00482   int mysendsize=1, myrecvsize=1;
00483 
00484   // Post Recvs
00485   start=0;
00486   for(i=0; i<NumRecvs; i++){
00487     if(RecvSizes) myrecvsize=RecvSizes[i]*SizeOfPacket;
00488     if(RecvProcs[i] != MyPID) {
00489       MPI_Irecv(RecvBuffer + start, myrecvsize, DataType, RecvProcs[i], msg_tag, comm, &requests[num_waits]);
00490       num_waits++;
00491     }
00492     else {
00493       self_recv_len = myrecvsize;
00494       self_recv_start=start;
00495     } 
00496     start+=myrecvsize;
00497   }
00498 
00499   // Do sends
00500   start=0;
00501   for(i=0; i<NumSends; i++){
00502     if(SendSizes) mysendsize=SendSizes[i]*SizeOfPacket;
00503     if(SendProcs[i] != MyPID)
00504       MPI_Send(SendBuffer + start, mysendsize,DataType,SendProcs[i],msg_tag,comm);
00505     else 
00506       self_send_start=start;
00507     start+=mysendsize;
00508   }
00509 
00510   // Self-copy (if needed)
00511   if(self_recv_len != -1)
00512     memcpy(RecvBuffer+self_recv_start,SendBuffer+self_send_start,self_recv_len*sizeof(MyType)*SizeOfPacket);
00513 
00514   // Wait
00515   if(NumRecvs > 0)
00516     MPI_Waitall(num_waits, &requests[0],&status[0]);
00517 }
00518 #endif
00519 
00520 
00521 #ifdef HAVE_MPI
00522 template <typename MyType>
00523 void boundary_exchange_varsize(const Epetra_MpiComm Comm, MPI_Datatype DataType, 
00524              int NumSends, const int * SendProcs, const int * SendSizes, MyType* SendBuffer, 
00525              int NumRecvs, const int * RecvProcs, int * RecvSizes, MyType*& RecvBuffer,int SizeOfPacket,int msg_tag)
00526 {  
00527 
00528   int i,rbuffersize=0;
00529 
00530   // Do a first round of boundary exchange with the the SendBuffer sizes
00531   boundary_exchange<int>(Comm,MPI_INT,NumSends,SendProcs,(int*)0,const_cast<int*>(SendSizes),NumRecvs,RecvProcs,(int*)0,RecvSizes,1,msg_tag);
00532 
00533   // Allocate the RecvBuffer 
00534   for(i=0; i<NumRecvs; i++) rbuffersize+=RecvSizes[i]*SizeOfPacket;
00535   RecvBuffer = new MyType[rbuffersize];
00536 
00537   // Do a second round of boundary exchange to trade the actual values
00538   boundary_exchange<MyType>(Comm,DataType,NumSends,SendProcs,SendSizes,SendBuffer,NumRecvs,RecvProcs,RecvSizes,RecvBuffer,SizeOfPacket,msg_tag+100);
00539 }
00540 #endif
00541 
00542 
00543 //=========================================================================
00544 //=========================================================================
00545 //=========================================================================
00546 LightweightMapData::LightweightMapData():
00547   Epetra_Data(),
00548   IndexBase_(0),
00549 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00550   LIDHash_int_(0),
00551 #endif
00552 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00553   LIDHash_LL_(0),
00554 #endif
00555   CopyMap_(0)
00556 {
00557 }
00558 //=========================================================================
00559 LightweightMapData::~LightweightMapData(){
00560 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00561   delete LIDHash_int_;
00562 #endif
00563 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00564   delete LIDHash_LL_;
00565 #endif
00566   delete CopyMap_;
00567 }
00568 
00569 //=========================================================================
00570 LightweightMap::LightweightMap():Data_(0),IsLongLong(false),IsInt(false){;}
00571 
00572 //=========================================================================
00573 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00574 void LightweightMap::Construct_int(int NumGlobalElements,int NumMyElements, const int * MyGlobalElements, long long IndexBase, bool GenerateHash)
00575 {
00576   Data_=new LightweightMapData();
00577   Data_->MyGlobalElements_int_.resize(NumMyElements);
00578 
00579   // Build the hash table
00580   if(GenerateHash) Data_->LIDHash_int_ = new Epetra_HashTable<int>(NumMyElements + 1 );
00581     for(int i=0; i < NumMyElements; ++i ) {
00582       Data_->MyGlobalElements_int_[i]=MyGlobalElements[i];
00583       if(GenerateHash) Data_->LIDHash_int_->Add(MyGlobalElements[i], i);
00584     }  
00585   IsLongLong = false;
00586   IsInt = true;
00587 }
00588 #endif
00589 
00590 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00591 void LightweightMap::Construct_LL(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, long long IndexBase, bool GenerateHash)
00592 {
00593   Data_=new LightweightMapData();
00594   Data_->MyGlobalElements_LL_.resize(NumMyElements);
00595 
00596   // Build the hash table
00597   if(GenerateHash) Data_->LIDHash_LL_ = new Epetra_HashTable<long long>(NumMyElements + 1 );
00598     for(int i=0; i < NumMyElements; ++i ) {
00599       Data_->MyGlobalElements_LL_[i]=MyGlobalElements[i];
00600       if(GenerateHash) Data_->LIDHash_LL_->Add(MyGlobalElements[i], i);
00601     }  
00602   IsLongLong = true;
00603   IsInt = false;
00604 }
00605 #endif
00606 
00607 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00608 LightweightMap::LightweightMap(int NumGlobalElements,int NumMyElements, const int * MyGlobalElements, int IndexBase, bool GenerateHash)
00609 {
00610   Construct_int(NumGlobalElements, NumMyElements, MyGlobalElements, IndexBase, GenerateHash);
00611 }
00612 #endif
00613 
00614 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00615 LightweightMap::LightweightMap(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, int IndexBase, bool GenerateHash)
00616 {
00617   Construct_LL(NumGlobalElements, NumMyElements, MyGlobalElements, IndexBase, GenerateHash);
00618 }
00619 
00620 LightweightMap::LightweightMap(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, long long IndexBase, bool GenerateHash)
00621 {
00622   Construct_LL(NumGlobalElements, NumMyElements, MyGlobalElements, IndexBase, GenerateHash);
00623 }
00624 #endif
00625 
00626 //=========================================================================
00627 LightweightMap::LightweightMap(const Epetra_Map & Map)
00628 {
00629   Data_=new LightweightMapData();
00630   Data_->CopyMap_=new Epetra_Map(Map);
00631   IsLongLong = Map.GlobalIndicesLongLong();
00632   IsInt = Map.GlobalIndicesInt();
00633 }
00634 
00635 //=========================================================================
00636 LightweightMap::LightweightMap(const LightweightMap& map)
00637   : Data_(map.Data_)
00638 {
00639   Data_->IncrementReferenceCount();
00640   IsLongLong = map.IsLongLong;
00641   IsInt = map.IsInt;
00642 }
00643 
00644 //=========================================================================
00645 LightweightMap::~LightweightMap(){
00646   CleanupData();
00647 }
00648 
00649 //=========================================================================
00650 LightweightMap & LightweightMap::operator= (const LightweightMap & map)
00651 {
00652   if((this != &map) && (Data_ != map.Data_)) {
00653     CleanupData();
00654     Data_ = map.Data_;
00655     Data_->IncrementReferenceCount();
00656   }
00657   IsLongLong = map.IsLongLong;
00658   IsInt = map.IsInt;
00659   return(*this);
00660 }
00661 
00662 //=========================================================================
00663 void LightweightMap::CleanupData(){
00664   if(Data_){
00665     Data_->DecrementReferenceCount();
00666     if(Data_->ReferenceCount() == 0) {
00667       delete Data_;
00668     }
00669   }
00670 }
00671 
00672 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00673 void LightweightMap::MyGlobalElementsPtr(int *& MyGlobalElementList) const {
00674   MyGlobalElementList = Data_->MyGlobalElements_int_.size() ? &Data_->MyGlobalElements_int_.front() : 0;
00675 }
00676 #endif
00677 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00678 void LightweightMap::MyGlobalElementsPtr(long long *& MyGlobalElementList) const {
00679   MyGlobalElementList = Data_->MyGlobalElements_LL_.size() ? &Data_->MyGlobalElements_LL_.front() : 0;
00680 }
00681 #endif
00682 
00683 //=========================================================================
00684 int LightweightMap::NumMyElements() const {
00685   if(Data_->CopyMap_) return Data_->CopyMap_->NumMyElements();
00686 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00687   if(IsInt)
00688     return Data_->MyGlobalElements_int_.size();
00689   else
00690 #endif
00691 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00692   if(IsLongLong)
00693     return Data_->MyGlobalElements_LL_.size();
00694   else
00695 #endif
00696     throw "EpetraExt::LightweightMap::NumMyElements: Global indices unknowns";
00697 }
00698 
00699 //=========================================================================
00700 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00701 int LightweightMap::LID(int GID) const {
00702   if(Data_->CopyMap_) return Data_->CopyMap_->LID(GID);
00703   if(IsInt)
00704     return Data_->LIDHash_int_->Get(GID);
00705   else if(IsLongLong)
00706     throw "EpetraExt::LightweightMap::LID: Int version called for long long map";
00707   else
00708     throw "EpetraExt::LightweightMap::LID: unknown GID type";
00709 }
00710 #endif
00711 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00712 int LightweightMap::LID(long long GID) const {
00713 
00714   if(Data_->CopyMap_) return Data_->CopyMap_->LID(GID);
00715   if(IsLongLong)
00716     return Data_->LIDHash_LL_->Get(GID);
00717   else if(IsInt)
00718     throw "EpetraExt::LightweightMap::LID: Long long version called for int map";
00719   else
00720     throw "EpetraExt::LightweightMap::LID: unknown GID type";
00721 }
00722 #endif
00723 
00724 //=========================================================================
00725 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00726 int LightweightMap::GID(int LID) const {
00727   if(Data_->CopyMap_) return Data_->CopyMap_->GID(LID);
00728   if(LID < 0 || LID > (int)Data_->MyGlobalElements_int_.size()) return -1;
00729   return Data_->MyGlobalElements_int_[LID];
00730 }
00731 #endif
00732 long long LightweightMap::GID64(int LID) const {
00733   if(Data_->CopyMap_) return Data_->CopyMap_->GID64(LID);
00734 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00735   if(IsInt) {
00736     if(LID < 0 || LID > (int)Data_->MyGlobalElements_int_.size()) return -1;
00737     return Data_->MyGlobalElements_int_[LID];
00738   }
00739   else
00740 #endif
00741 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00742   if(IsLongLong) {
00743     if(LID < 0 || LID > (int)Data_->MyGlobalElements_LL_.size()) return -1;
00744     return Data_->MyGlobalElements_LL_[LID];
00745   }
00746   else
00747 #endif
00748     throw "EpetraExt::LightweightMap::GID64: Global indices unknown.";
00749 }
00750 
00751 //=========================================================================
00752 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00753 int* LightweightMap::MyGlobalElements() const {
00754   if(Data_->CopyMap_) return Data_->CopyMap_->MyGlobalElements();
00755   else if(Data_->MyGlobalElements_int_.size()>0) return const_cast<int*>(&Data_->MyGlobalElements_int_[0]);
00756   else return 0;
00757 }
00758 #endif
00759 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00760 long long* LightweightMap::MyGlobalElements64() const {
00761   if(Data_->CopyMap_) return Data_->CopyMap_->MyGlobalElements64();
00762   else if(Data_->MyGlobalElements_LL_.size()>0) return const_cast<long long*>(&Data_->MyGlobalElements_LL_[0]);
00763   else return 0;
00764 }
00765 #endif
00766 
00767 //=========================================================================
00768 int LightweightMap::MinLID() const {
00769   if(Data_->CopyMap_) return Data_->CopyMap_->MinLID();
00770   else return 0;
00771 }
00772 
00773 //=========================================================================
00774 int LightweightMap::MaxLID() const {
00775   if(Data_->CopyMap_) return Data_->CopyMap_->MaxLID();
00776   else {
00777 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00778   if(IsInt)
00779     return Data_->MyGlobalElements_int_.size() - 1;
00780   else
00781 #endif
00782 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00783   if(IsLongLong)
00784     return Data_->MyGlobalElements_LL_.size() - 1;
00785   else
00786 #endif
00787     throw "EpetraExt::LightweightMap::MaxLID: Global indices unknowns";
00788   }   
00789 }
00790 
00791 
00792 //=========================================================================
00793 //=========================================================================
00794 //=========================================================================
00795 RemoteOnlyImport::RemoteOnlyImport(const Epetra_Import & Importer, LightweightMap & RemoteOnlyTargetMap)
00796 {
00797   int i;
00798 
00799   // Build an "Importer" that only takes the remote parts of the Importer.
00800   SourceMap_=&Importer.SourceMap();
00801   TargetMap_=&RemoteOnlyTargetMap;
00802 
00803   // Pull data from the Importer
00804   NumSend_            = Importer.NumSend();
00805   NumRemoteIDs_       = Importer.NumRemoteIDs();
00806   NumExportIDs_       = Importer.NumExportIDs();
00807   Distor_             = &Importer.Distributor();
00808   int * OldRemoteLIDs = Importer.RemoteLIDs();
00809   int * OldExportLIDs = Importer.ExportLIDs();
00810   int * OldExportPIDs = Importer.ExportPIDs();
00811   
00812   // Sanity Check
00813   if(NumRemoteIDs_ != RemoteOnlyTargetMap.NumMyElements())
00814     throw std::runtime_error("RemoteOnlyImport: Importer doesn't match RemoteOnlyTargetMap for number of remotes.");
00815 
00816   // Copy the ExportIDs_, since they don't change
00817   ExportLIDs_ = new int[NumExportIDs_];
00818   ExportPIDs_ = new int[NumExportIDs_];
00819   for(i=0; i<NumExportIDs_; i++) {
00820     ExportLIDs_[i] = OldExportLIDs[i];
00821     ExportPIDs_[i] = OldExportPIDs[i];
00822   }
00823 
00824   // The RemoteIDs, on the other hand, do change.  So let's do this right.
00825   // Note: We might be able to bypass the LID call by just indexing off the Same and Permute GIDs, but at the moment this
00826   // is fast enough not to worry about it.
00827   RemoteLIDs_ = new int[NumRemoteIDs_];
00828   if(TargetMap_->GlobalIndicesInt()) {
00829     for(i=0; i<NumRemoteIDs_; i++) 
00830       RemoteLIDs_[i] = TargetMap_->LID( (int) Importer.TargetMap().GID64(OldRemoteLIDs[i]));
00831   }
00832   else if(TargetMap_->GlobalIndicesLongLong()) {
00833     for(i=0; i<NumRemoteIDs_; i++) 
00834       RemoteLIDs_[i] = TargetMap_->LID(Importer.TargetMap().GID64(OldRemoteLIDs[i]));
00835   }
00836   else
00837     throw std::runtime_error("RemoteOnlyImport: TargetMap_ index type unknown.");
00838 
00839   // Nowe we make sure these guys are in sorted order.  AztecOO, ML and all that jazz.
00840   for(i=0; i<NumRemoteIDs_-1; i++) 
00841     if(RemoteLIDs_[i] > RemoteLIDs_[i+1])
00842       throw std::runtime_error("RemoteOnlyImport: Importer and RemoteOnlyTargetMap order don't match.");
00843 }
00844 
00845 //=========================================================================
00846 RemoteOnlyImport::~RemoteOnlyImport()
00847 {
00848   delete [] ExportLIDs_;
00849   delete [] ExportPIDs_;
00850   delete [] RemoteLIDs_;
00851   // Don't delete the Distributor, SourceMap_ or TargetMap_ - those were shallow copies
00852 }
00853 
00854 //=========================================================================
00855 //=========================================================================
00856 //=========================================================================
00857 
00858 template <class GO>
00859 void MakeColMapAndReindexSort(int& NumRemoteColGIDs, GO*& RemoteColindices,
00860     std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs);
00861 
00862 template <>
00863 void MakeColMapAndReindexSort<int>(int& NumRemoteColGIDs, int*& RemoteColindices,
00864     std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
00865 {
00866   // Sort External column indices so that all columns coming from a given remote processor are contiguous
00867   int NLists = 2;
00868   int* SortLists[3]; // this array is allocated on the stack, and so we won't need to delete it.
00869   if(NumRemoteColGIDs > 0) {
00870     SortLists[0] = RemoteColindices;    
00871     SortLists[1] = &RemotePermuteIDs[0];
00872     Epetra_Util::Sort(true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, NLists, SortLists);
00873   }
00874 
00875 
00876   bool SortGhostsAssociatedWithEachProcessor_ = false;
00877   if (SortGhostsAssociatedWithEachProcessor_) {
00878     // Sort external column indices so that columns from a given remote processor are not only contiguous
00879     // but also in ascending order. NOTE: I don't know if the number of externals associated
00880     // with a given remote processor is known at this point ... so I count them here.
00881 
00882     NLists=1;
00883     int StartCurrent, StartNext;
00884     StartCurrent = 0; StartNext = 1;
00885     while ( StartNext < NumRemoteColGIDs ) {
00886       if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
00887       else {
00888   SortLists[0] =  &RemotePermuteIDs[StartCurrent];
00889   Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists);
00890         StartCurrent = StartNext; StartNext++;
00891       }
00892     }
00893     SortLists[0] =  &RemotePermuteIDs[StartCurrent];
00894     Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists);
00895   }
00896 }
00897 
00898 template <>
00899 void MakeColMapAndReindexSort<long long>(int& NumRemoteColGIDs, long long*& RemoteColindices,
00900     std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
00901 {
00902   // Sort External column indices so that all columns coming from a given remote processor are contiguous
00903   if(NumRemoteColGIDs > 0) {
00904     long long* SortLists_LL[1] = {RemoteColindices};
00905   int* SortLists_int[1] = {&RemotePermuteIDs[0]};
00906     Epetra_Util::Sort(true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, 1, SortLists_int, 1, SortLists_LL);
00907   }
00908 
00909   bool SortGhostsAssociatedWithEachProcessor_ = false;
00910   if (SortGhostsAssociatedWithEachProcessor_) {
00911     // Sort external column indices so that columns from a given remote processor are not only contiguous
00912     // but also in ascending order. NOTE: I don't know if the number of externals associated
00913     // with a given remote processor is known at this point ... so I count them here.
00914 
00915     int NLists=1;
00916     int StartCurrent, StartNext;
00917     StartCurrent = 0; StartNext = 1;
00918     while ( StartNext < NumRemoteColGIDs ) {
00919       if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
00920       else {
00921         int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
00922         Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
00923         StartCurrent = StartNext; StartNext++;
00924       }
00925     }
00926   int* SortLists[1] =  {&RemotePermuteIDs[StartCurrent]};
00927     Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
00928   }
00929 }
00930 
00931 template <class GO>
00932 int LightweightCrsMatrix::MakeColMapAndReindex(std::vector<int> owningPIDs, std::vector<GO> Gcolind)
00933 {
00934   int i,j;
00935 
00936 #ifdef ENABLE_MMM_TIMINGS
00937   using Teuchos::TimeMonitor;
00938   Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.1")));
00939 #endif
00940 
00941   // Scan all column indices and sort into two groups: 
00942   // Local:  those whose GID matches a GID of the domain map on this processor and
00943   // Remote: All others.
00944   int numDomainElements = DomainMap_.NumMyElements();
00945   std::vector<bool> LocalGIDs(numDomainElements,false);
00946 
00947   bool DoSizes = !DomainMap_.ConstantElementSize(); // If not constant element size, then error
00948   if(DoSizes) EPETRA_CHK_ERR(-1);
00949 
00950   // In principle it is good to have RemoteGIDs and RemoteGIDList be as long as the number of remote GIDs
00951   // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
00952   // and the number of block rows.
00953   int numMyBlockRows;
00954   if(use_lw) numMyBlockRows = RowMapLW_->NumMyElements();
00955   else numMyBlockRows = RowMapEP_->NumMyElements();
00956 
00957   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
00958   Epetra_HashTable<GO> RemoteGIDs(hashsize); 
00959   std::vector<GO>  RemoteGIDList;     RemoteGIDList.reserve(hashsize);
00960   std::vector<int> RemoteOwningPIDs;  RemoteOwningPIDs.reserve(hashsize);
00961     
00962   // In order to do the map reindexing inexpensively, we clobber the GIDs during this pass.  For *local* GID's we clobber them
00963   // with their LID in the domainMap.  For *remote* GIDs, we clobber them with (numDomainElements+NumRemoteColGIDs) before the increment of
00964   // the remote count.  These numberings will be separate because no local LID is greater than numDomainElements. 
00965   int NumLocalColGIDs = 0;
00966   int NumRemoteColGIDs = 0;
00967   for(i = 0; i < numMyBlockRows; i++) {
00968     for(j = rowptr_[i]; j < rowptr_[i+1]; j++) {
00969       GO GID = Gcolind[j];
00970       // Check if GID matches a row GID
00971       int LID = DomainMap_.LID(GID);
00972       if(LID != -1) {
00973   bool alreadyFound = LocalGIDs[LID];
00974   if (!alreadyFound) {
00975           LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
00976           NumLocalColGIDs++;
00977   }
00978   colind_[j] = LID; 
00979       }
00980       else {
00981   GO hash_value=RemoteGIDs.Get(GID);
00982   if(hash_value  == -1) { // This means its a new remote GID
00983     int PID = owningPIDs[j];
00984     if(PID==-1) printf("[%d] ERROR: Remote PID should not be -1\n",DomainMap_.Comm().MyPID());
00985     colind_[j] = numDomainElements + NumRemoteColGIDs;
00986     RemoteGIDs.Add(GID, NumRemoteColGIDs);
00987     RemoteGIDList.push_back(GID);
00988     RemoteOwningPIDs.push_back(PID);
00989     NumRemoteColGIDs++;
00990   }
00991   else
00992     colind_[j] = numDomainElements + hash_value;    
00993       }
00994     }
00995   }
00996 
00997   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
00998   if (DomainMap_.Comm().NumProc()==1) {     
00999     if (NumRemoteColGIDs!=0) {
01000       throw "Some column IDs are not in domainMap.  If matrix is rectangular, you must pass in domainMap to FillComplete";
01001       // Sanity test: When one processor,there can be no remoteGIDs
01002     }
01003     if (NumLocalColGIDs==numDomainElements) {
01004       ColMap_ = DomainMap_;
01005 
01006       // In this case, we just use the domainMap's indices, which is, not coincidently, what we clobbered colind_ with up above anyway. 
01007       // No further reindexing is needed.
01008       return(0); 
01009     }
01010   }
01011       
01012   // Now build integer array containing column GIDs
01013   // Build back end, containing remote GIDs, first
01014   int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
01015   typename Epetra_GIDTypeSerialDenseVector<GO>::impl Colindices;
01016   if(numMyBlockCols > 0) 
01017     Colindices.Size(numMyBlockCols);
01018   GO* RemoteColindices = Colindices.Values() + NumLocalColGIDs; // Points to back end of Colindices
01019 
01020   for(i = 0; i < NumRemoteColGIDs; i++) 
01021     RemoteColindices[i] = RemoteGIDList[i]; 
01022 
01023   // Build permute array for *remote* reindexing.
01024   std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
01025   for(i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
01026 
01027   MakeColMapAndReindexSort<GO>(NumRemoteColGIDs, RemoteColindices, RemotePermuteIDs, RemoteOwningPIDs);
01028 
01029   // Reverse the permutation to get the information we actually care about
01030   std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
01031   for(i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
01032 
01033   // Build permute array for *local* reindexing.
01034   bool use_local_permute=false;
01035   std::vector<int> LocalPermuteIDs(numDomainElements);
01036 
01037   // Now fill front end. Two cases:
01038   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
01039   //     can simply read the domain GIDs into the front part of Colindices, otherwise 
01040   // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
01041   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
01042 
01043   if(NumLocalColGIDs == DomainMap_.NumMyElements()) {
01044     DomainMap_.MyGlobalElements(Colindices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
01045   }
01046   else {
01047     GO* MyGlobalElements = 0;
01048   DomainMap_.MyGlobalElementsPtr(MyGlobalElements);
01049     int NumLocalAgain = 0;
01050     use_local_permute = true;    
01051     for(i = 0; i < numDomainElements; i++) {
01052       if(LocalGIDs[i]) {
01053   LocalPermuteIDs[i] = NumLocalAgain;
01054   Colindices[NumLocalAgain++] = MyGlobalElements[i];
01055       }
01056     }
01057     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
01058   }
01059 
01060 
01061   // Copy the remote PID list correctly
01062   ColMapOwningPIDs_.resize(numMyBlockCols);
01063   ColMapOwningPIDs_.assign(numMyBlockCols,DomainMap_.Comm().MyPID());
01064   for(int i=0;i<NumRemoteColGIDs;i++)
01065     ColMapOwningPIDs_[NumLocalColGIDs+i] = RemoteOwningPIDs[i];
01066 
01067 #ifdef ENABLE_MMM_TIMINGS
01068   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.2")));
01069 #endif
01070 
01071   // Make Column map with same element sizes as Domain map 
01072   LightweightMap temp((GO) -1, numMyBlockCols, Colindices.Values(), (GO) DomainMap_.IndexBase64());
01073   ColMap_ = temp;
01074 
01075 
01076 #ifdef ENABLE_MMM_TIMINGS
01077   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.3")));
01078 #endif
01079 
01080   // Low-cost reindex of the matrix
01081   for(i=0; i<numMyBlockRows; i++){
01082     for(j=rowptr_[i]; j<rowptr_[i+1]; j++){
01083       int ID=colind_[j];
01084       if(ID < numDomainElements){
01085   if(use_local_permute) colind_[j] = LocalPermuteIDs[colind_[j]];
01086   // In the case where use_local_permute==false, we just copy the DomainMap's ordering, which it so happens
01087   // is what we put in colind_ to begin with.
01088       }
01089       else
01090   colind_[j] =  NumLocalColGIDs + ReverseRemotePermuteIDs[colind_[j]-numDomainElements];
01091     }
01092   }
01093 
01094   assert((size_t)ColMap_.NumMyElements() == ColMapOwningPIDs_.size());
01095 
01096   return(0);
01097 }
01098 
01099 
01100 //=========================================================================
01101 // Template params are <PID,GID>
01102 static inline bool lessthan12(std::pair<int,int> i, std::pair<int,int> j){
01103   return ((i.first<j.first) || (i.first==j.first && i.second <j.second));
01104 }
01105 
01106 
01107 template<typename ImportType, typename int_type>
01108 int LightweightCrsMatrix::PackAndPrepareReverseComm(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter,
01109                 std::vector<int> &ReverseSendSizes, std::vector<int_type> &ReverseSendBuffer) {
01110 #ifdef HAVE_MPI
01111   // Buffer pairs are in (PID,GID) order
01112   int i,j,k;
01113   const Epetra_Import *MyImporter= SourceMatrix.Importer();
01114   if(MyImporter == 0) return -1;
01115   const Epetra_MpiComm * MpiComm        = dynamic_cast<const Epetra_MpiComm*>(&SourceMatrix.Comm());
01116   int MyPID                             = MpiComm->MyPID();
01117 
01118   // Things related to messages I and sending in forward mode (RowImporter)
01119   int NumExportIDs                      = RowImporter.NumExportIDs();
01120   int* ExportLIDs                       = RowImporter.ExportLIDs();
01121   int* ExportPIDs                       = RowImporter.ExportPIDs();
01122 
01123   // Things related to messages I am sending in reverse mode (MyImporter)
01124   Epetra_Distributor& Distor            = MyImporter->Distributor();
01125   const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
01126   int NumRecvs                          = MDistor->NumReceives();
01127   int* RemoteLIDs                       = MyImporter->RemoteLIDs();
01128   const int * ProcsFrom                 = MDistor->ProcsFrom();
01129   const int * LengthsFrom               = MDistor->LengthsFrom();
01130 
01131   
01132   // Get the owning pids in a special way, s.t. ProcsFrom[RemotePIDs[i]] is the guy who owns 
01133   // RemoteLIDs[j]....
01134   std::vector<int> RemotePIDOrder(SourceMatrix.NumMyCols(),-1);
01135 
01136   // Now, for each remote ID, record which processor (in ProcsFrom ordering) owns it.
01137   for(i=0,j=0;i<NumRecvs;i++){
01138     for(k=0;k<LengthsFrom[i];k++){
01139       int pid=ProcsFrom[i];
01140       if(pid!=MyPID) RemotePIDOrder[RemoteLIDs[j]]=i;
01141       j++;
01142     }    
01143   }
01144 
01145   // Step One: Start tacking the (GID,PID) pairs on the std sets
01146   std::vector<std::set<std::pair<int,int_type> > > ReversePGIDs(NumRecvs);
01147   int    *rowptr, *colind;
01148   double *vals;
01149   EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
01150   
01151 
01152   // Loop over each exported row and add to the temp list
01153   for(i=0; i < NumExportIDs; i++) {
01154     int lid = ExportLIDs[i];
01155     int exp_pid = ExportPIDs[i];    
01156     for(j=rowptr[lid]; j<rowptr[lid+1]; j++){
01157       int pid_order = RemotePIDOrder[colind[j]];
01158       if(pid_order!=-1) {
01159   int_type gid = (int_type) SourceMatrix.GCID64(colind[j]);
01160   // This GID is getting shipped off somewhere
01161   ReversePGIDs[pid_order].insert(std::pair<int,int_type>(exp_pid,gid));
01162       }
01163     }
01164   }
01165 
01166   // Step 2: Count sizes (one too large to avoid std::vector errors)
01167   ReverseSendSizes.resize(NumRecvs+1);
01168   int totalsize=0;
01169   for(i=0; i<NumRecvs; i++) {
01170     ReverseSendSizes[i] = 2*ReversePGIDs[i].size();
01171     totalsize += ReverseSendSizes[i];
01172   }
01173 
01174   // Step 3: Alloc and fill the send buffer (one too large to avoid std::vector errors)
01175   ReverseSendBuffer.resize(totalsize+1);
01176   for(i=0, j=0; i<NumRecvs; i++) {
01177     for(typename std::set<std::pair<int,int_type> >::iterator it=ReversePGIDs[i].begin(); it!=ReversePGIDs[i].end(); it++) {
01178       ReverseSendBuffer[j]   =  it->first;
01179       ReverseSendBuffer[j+1] =  it->second;
01180       j+=2;
01181     }
01182   }
01183 #endif
01184 
01185   return 0;
01186 }
01187 
01188 //=========================================================================
01189 
01190 template<typename int_type>
01191 void build_type3_exports_sort(std::vector<int_type>& ExportGID3, std::vector<int> &ExportPID3, int total_length3);
01192 
01193 template<>
01194 void build_type3_exports_sort<int>(std::vector<int>& ExportGID3, std::vector<int> &ExportPID3, int total_length3)
01195 {
01196   int * companion = &ExportGID3[0];
01197   Epetra_Util::Sort(true,total_length3,&ExportPID3[0],0,0,1,&companion,0,0);
01198 }
01199 
01200 template<>
01201 void build_type3_exports_sort<long long>(std::vector<long long>& ExportGID3, std::vector<int> &ExportPID3, int total_length3)
01202 {
01203   long long * companion = &ExportGID3[0];
01204   Epetra_Util::Sort(true,total_length3,&ExportPID3[0],0,0,0,0,1,&companion);
01205 }
01206 
01207 template<typename int_type>
01208 int build_type3_exports(int MyPID,int Nrecv, Epetra_BlockMap &DomainMap, std::vector<int> &ReverseRecvSizes, const int_type *ReverseRecvBuffer,  std::vector<int> &ExportLID3, std::vector<int> &ExportPID3){
01209   int i,j;
01210 
01211   // Estimate total length of procs_to for Type 3
01212   int total_length3=0;
01213   for(i=0; i<Nrecv; i++)
01214     total_length3+=ReverseRecvSizes[i]/2;
01215   if(total_length3==0) return 0;
01216 
01217   std::vector<int_type> ExportGID3(total_length3);
01218   ExportLID3.resize(total_length3);
01219   ExportPID3.resize(total_length3);
01220 
01221   // Build a sorted colmap-style list for Type3 (removing any self-sends)    
01222   for(i=0,j=0; i<2*total_length3; i+=2) {
01223     if(ReverseRecvBuffer[i] != MyPID){
01224       ExportPID3[j]=ReverseRecvBuffer[i];
01225       ExportGID3[j]=ReverseRecvBuffer[i+1];
01226       j++;
01227     }    
01228   }
01229   total_length3=j;
01230 
01231   if(total_length3==0) return 0;
01232 
01233   // Sort (ala Epetra_CrsGraph)
01234   build_type3_exports_sort<int_type>(ExportGID3, ExportPID3, total_length3);
01235   int StartCurrent, StartNext;
01236   StartCurrent = 0; StartNext = 1;
01237   while ( StartNext < total_length3 ) {
01238     if(ExportPID3[StartNext] == ExportPID3[StartNext-1]) StartNext++;
01239     else {
01240       Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
01241       StartCurrent = StartNext; StartNext++;
01242     }
01243   }
01244   Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
01245   
01246 
01247   /*  printf("[%d] Type 3 Sorted= ",MyComm.MyPID());
01248   for(i=0; i<total_length3; i++)
01249     printf("(--,%2d,%2d) ",ExportGID3[i],ExportPID3[i]);
01250   printf("\n");
01251   fflush(stdout);*/
01252 
01253 
01254   // Uniq & resize
01255   for(i=1,j=1; i<total_length3; i++){
01256     if(ExportPID3[i]!=ExportPID3[i-1] || ExportGID3[i]!=ExportGID3[i-1]){
01257       ExportPID3[j] = ExportPID3[i];
01258       ExportGID3[j] = ExportGID3[i];
01259       j++;
01260     }
01261   }
01262   ExportPID3.resize(j);
01263   ExportLID3.resize(j);
01264   total_length3=j;
01265 
01266   /*  printf("[%d] Type 3 UNIQ  = ",MyComm.MyPID());
01267   for(i=0; i<total_length3; i++)
01268     printf("(--,%2d,%2d) ",ExportGID3[i],ExportPID3[i]);
01269   printf("\n");
01270   fflush(stdout);*/
01271 
01272 
01273   
01274   // Now index down to LIDs   
01275   for(i=0; i<total_length3; i++) {
01276     ExportLID3[i]=DomainMap.LID(ExportGID3[i]);
01277     if(ExportLID3[i] < 0) throw std::runtime_error("LightweightCrsMatrix:MakeExportLists invalid LID");    
01278   }
01279     
01280   /*  printf("[%d] Type 3 FINAL = ",MyComm.MyPID());
01281   for(i=0; i<total_length3; i++)
01282     printf("(%2d,%2d,%2d) ",ExportLID3[i],ExportGID3[i],ExportPID3[i]);
01283   printf("\n");
01284   fflush(stdout);*/
01285  
01286   return total_length3;
01287 }
01288 
01289 //=========================================================================
01290 template<typename ImportType, typename int_type>
01291 int build_type2_exports(const Epetra_CrsMatrix & SourceMatrix, ImportType & MyImporter, std::vector<int> &ExportLID2, std::vector<int> &ExportPID2){
01292   int total_length2=0;
01293 #ifdef HAVE_MPI
01294   int i,j;
01295   const Epetra_Import *SourceImporter= SourceMatrix.Importer();
01296 
01297   int    *rowptr, *colind;
01298   double *vals;
01299   EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
01300 
01301   // Things related to messages I am sending in forward mode (MyImporter)
01302   int NumExportIDs                      = MyImporter.NumExportIDs();
01303   const int* ExportLIDs                 = MyImporter.ExportLIDs();
01304   const int* ExportPIDs                 = MyImporter.ExportPIDs();
01305   if(NumExportIDs==0) return 0;
01306 
01307   Epetra_Distributor& Distor            = MyImporter.Distributor();
01308   const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
01309 
01310   // Assume I own all the cols, then flag any cols I don't own
01311   // This allows us to avoid LID calls later...
01312   std::vector<bool> IsOwned(SourceMatrix.NumMyCols(),true);
01313   if(SourceImporter) {
01314     const int * RemoteLIDs = SourceImporter->RemoteLIDs();
01315     // Now flag the cols I don't own
01316     for(i=0; i<SourceImporter->NumRemoteIDs(); i++) 
01317       IsOwned[RemoteLIDs[i]]=false;
01318   }
01319 
01320   // Status vector 
01321   std::vector<int> SentTo(SourceMatrix.NumMyCols(),-1);
01322 
01323   // Initial allocation: NumUnknowsSent * Max row size (guaranteed to be too large)
01324   for(i=0,total_length2=0; i<MDistor->NumSends(); i++) total_length2+= MDistor->LengthsTo()[i] * SourceMatrix.MaxNumEntries();
01325   std::vector<int_type> ExportGID2(total_length2);
01326 
01327   ExportLID2.resize(total_length2);
01328   ExportPID2.resize(total_length2);
01329 
01330   int current=0, last_start=0, last_pid=ExportPIDs[0];
01331   for(i=0; i<NumExportIDs; i++){
01332     // For each row I have to send via MyImporter...
01333     int row=ExportLIDs[i];
01334     int pid=ExportPIDs[i];
01335 
01336     if(i!=0 && pid>last_pid) {
01337       // We have a new PID, so lets finish up the current one      
01338       if(current!=last_start){
01339         int *lids = &ExportLID2[last_start];
01340   Epetra_Util::Sort(true,current-last_start,&ExportGID2[last_start],0,0,1,&lids,0,0);
01341   // Note: we don't need to sort the ExportPIDs since they're the same since last_start
01342       }
01343       // Reset the list
01344       last_pid=pid;
01345       last_start=current;         
01346     }
01347     else if(pid < last_pid) {
01348       throw std::runtime_error("build_type2_exports: ExportPIDs are not sorted!");
01349     }
01350 
01351     for(j=rowptr[row]; j<rowptr[row+1]; j++) {
01352       // For each column in that row...      
01353       int col=colind[j];
01354       if(IsOwned[col] && SentTo[col]!=pid){
01355   // We haven't added this guy to the list yet.
01356   if(current>= total_length2) throw std::runtime_error("build_type2_exports: More export ids than I thought!");
01357   SentTo[col]         = pid;
01358   ExportGID2[current] = (int_type) SourceMatrix.GCID64(col);
01359   ExportLID2[current] = SourceMatrix.DomainMap().LID(ExportGID2[current]);
01360   ExportPID2[current] = pid;
01361   current++;
01362       }
01363     }
01364   }//end main loop
01365 
01366   // Final Sort
01367   int *lids = &ExportLID2[last_start];
01368   Epetra_Util::Sort(true,current-last_start,&ExportGID2[last_start],0,0,1,&lids,0,0);
01369   // Note: we don't need to sort the ExportPIDs since they're the same since last_start
01370 
01371   total_length2=current;
01372   ExportLID2.resize(total_length2);
01373   ExportPID2.resize(total_length2);
01374 #endif
01375   return total_length2;
01376 }
01377 
01378 //=========================================================================
01379 template<typename int_type>
01380 void build_type1_exports_sort(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int_type>& ExportGID1, int total_length1);
01381 
01382 template<>
01383 void build_type1_exports_sort<int>(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int>& ExportGID1, int total_length1){
01384   int * companion[2] = {&ExportLID1[0],&ExportGID1[0]};
01385   Epetra_Util::Sort(true,total_length1,&ExportPID1[0],0,0,2,&companion[0],0,0);
01386 }
01387 
01388 template<>
01389 void build_type1_exports_sort<long long>(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<long long>& ExportGID1, int total_length1){
01390   int * companion = &ExportLID1[0];
01391   long long * companion64 = &ExportGID1[0];
01392   Epetra_Util::Sort(true,total_length1,&ExportPID1[0],0,0,1,&companion,1,&companion64);
01393 }
01394 
01395 template<typename int_type>
01396 int build_type1_exports(const Epetra_Import * Importer1, std::vector<int> &ExportLID1, std::vector<int> &ExportPID1){
01397   int i, total_length1=0;
01398   if(!Importer1) return 0;
01399   total_length1 = Importer1->NumSend();
01400   if(total_length1==0) return 0;
01401 
01402   std::vector<int_type> ExportGID1(total_length1);
01403   ExportLID1.resize(total_length1);
01404   ExportPID1.resize(total_length1);
01405   const int * ExportLID1Base = Importer1->ExportLIDs();
01406   const int * ExportPID1Base = Importer1->ExportPIDs();
01407 
01408   for(i=0; i<total_length1; i++){
01409     ExportLID1[i] = ExportLID1Base[i];
01410     ExportPID1[i] = ExportPID1Base[i];
01411     ExportGID1[i] = (int_type) Importer1->SourceMap().GID64(ExportLID1Base[i]);
01412   }
01413 
01414   // Sort (ala Epetra_CrsGraph)
01415   build_type1_exports_sort<int_type>(ExportLID1, ExportPID1, ExportGID1, total_length1);
01416 
01417   int StartCurrent, StartNext;
01418   StartCurrent = 0; StartNext = 1;
01419   while ( StartNext < total_length1 ) {
01420     if(ExportPID1[StartNext] == ExportPID1[StartNext-1]) StartNext++;
01421     else {
01422       int *new_companion = {&ExportLID1[StartCurrent]};
01423       Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
01424       StartCurrent = StartNext; StartNext++;
01425     }
01426   }
01427   int *new_companion = {&ExportLID1[StartCurrent]};
01428   Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
01429   return total_length1;
01430 }
01431 
01432 //=========================================================================
01433 template<typename ImportType, typename int_type>
01434 int LightweightCrsMatrix::MakeExportLists(const Epetra_CrsMatrix & SourceMatrix, ImportType & Importer2,
01435             std::vector<int> &ReverseRecvSizes, const int_type *ReverseRecvBuffer,
01436             std::vector<int> & ExportPIDs, std::vector<int> & ExportLIDs) {
01437 #ifdef HAVE_MPI
01438   int MyPID = SourceMatrix.Comm().MyPID();
01439 
01440   // This could (legitimately) be zero, in which case we don't have any ReverseRecvs either.
01441   const Epetra_Import *Importer1= SourceMatrix.Importer();
01442 
01443   // So for each entry in the DomainMap, I have to answer the question: Do I need to send this to anybody?  And if so, to whom?
01444   //
01445   // This information comes from three sources:
01446   // 1) IDs in my DomainMap that are in someone else's ColMap (e.g. SourceMatrix.Importer()).
01447   // 2) IDs that I own that I sent to another proc in the Forward communication round (e.g. RowImporter).
01448   // 3) IDs that someone else sent on during the Forward round (and they told me about duing the reverse round).
01449   //
01450   // Any of these could legitimately be null.
01451   Epetra_MpiDistributor * Distor1 = (Importer1)?(dynamic_cast<Epetra_MpiDistributor*>(&Importer1->Distributor())):0;
01452 
01453   int Nsend1 = (Distor1)?(Distor1->NumSends()):0; // Also the number of messages we'll need to parse through in build_type3_exports
01454 
01455   std::vector<int> ExportPID3;
01456   std::vector<int> ExportLID3;
01457 
01458   std::vector<int> ExportPID2;
01459   std::vector<int> ExportLID2;
01460 
01461   std::vector<int> ExportPID1;
01462   std::vector<int> ExportLID1;
01463 
01464   // Build (sorted) ExportLID / ExportGID list for Type
01465   int Len1=build_type1_exports<int_type>(Importer1, ExportLID1, ExportPID1);
01466   int Len2=build_type2_exports<ImportType, int_type>(SourceMatrix, Importer2, ExportLID2, ExportPID2);
01467   int Len3=build_type3_exports(MyPID,Nsend1,DomainMap_,ReverseRecvSizes, ReverseRecvBuffer, ExportLID3, ExportPID3);
01468 
01469   // Since everything should now be sorted correctly, we can do a streaming merge of the three Export lists...
01470 #ifdef HAVE_EPETRAEXT_DEBUG
01471   {
01472     int i;
01473     // Now we sanity check the 1 & 2 lists
01474     bool test_passed=true;
01475     for(i=1; i<Len1; i++) {
01476       if(ExportPID1[i] < ExportPID1[i-1] || (ExportPID1[i] == ExportPID1[i-1] && DomainMap_.GID(ExportLID1[i]) < DomainMap_.GID(ExportLID1[i-1])))
01477   test_passed=false;
01478     }
01479     SourceMatrix.Comm().Barrier();
01480     if(!test_passed) {
01481       printf("[%d] Type1 ERRORS  = ",SourceMatrix.Comm().MyPID());
01482       for(int i=0; i<Len1; i++)
01483   printf("(%2d,%2d,%2d) ",ExportLID1[i],DomainMap_.GID(ExportLID1[i]),ExportPID1[i]);
01484       printf("\n");
01485       fflush(stdout);      
01486       throw std::runtime_error("Importer1 fails the sanity test");      
01487     }
01488 
01489     for(i=1; i<Len2; i++) {
01490       if(ExportPID2[i] < ExportPID2[i-1]  || (ExportPID2[i] == ExportPID2[i-1] && DomainMap_.GID(ExportLID2[i]) < DomainMap_.GID(ExportLID2[i-1])))
01491   test_passed=false;
01492     }
01493 
01494     SourceMatrix.Comm().Barrier();
01495     if(!test_passed) {
01496       printf("[%d] Type2 ERRORS  = ",SourceMatrix.Comm().MyPID());
01497       for(int i=0; i<Len2; i++)
01498   printf("(%2d,%2d,%2d) ",ExportLID2[i],DomainMap_.GID(ExportLID2[i]),ExportPID2[i]);
01499       printf("\n");
01500       fflush(stdout);      
01501       throw std::runtime_error("Importer2 fails the sanity test");      
01502     }
01503 
01504     for(i=1; i<Len3; i++) {
01505       if(ExportPID3[i] < ExportPID3[i-1]  || (ExportPID3[i] == ExportPID3[i-1] && DomainMap_.GID(ExportLID3[i]) < DomainMap_.GID(ExportLID3[i-1])))
01506   test_passed=false;
01507     }
01508 
01509     SourceMatrix.Comm().Barrier();
01510     if(!test_passed) {
01511       printf("[%d] Type3 ERRORS  = ",SourceMatrix.Comm().MyPID());
01512       for(int i=0; i<Len3; i++)
01513   printf("(%2d,%2d,%2d) ",ExportLID3[i],DomainMap_.GID(ExportLID3[i]),ExportPID3[i]);
01514       printf("\n");
01515       fflush(stdout);      
01516       throw std::runtime_error("Importer3 fails the sanity test");      
01517     }
01518 
01519 
01520   }
01521 #endif
01522 
01523   if(Importer1 && !Importer1->SourceMap().SameAs(DomainMap_))
01524     throw std::runtime_error("ERROR: Map Mismatch Importer1");
01525 
01526   if(!Importer2.SourceMap().SameAs(SourceMatrix.RowMap()))
01527     throw std::runtime_error("ERROR: Map Mismatch Importer2");
01528 
01529   int_type InfGID = std::numeric_limits<int_type>::max();
01530   int InfPID = INT_MAX;
01531 
01532   int i1=0, i2=0, i3=0, current=0;
01533 
01534   int MyLen=Len1+Len2+Len3;  
01535   ExportLIDs.resize(MyLen);
01536   ExportPIDs.resize(MyLen);
01537 
01538   while(i1 < Len1 || i2 < Len2 || i3 < Len3){
01539     int PID1 = (i1<Len1)?(ExportPID1[i1]):InfPID;
01540     int PID2 = (i2<Len2)?(ExportPID2[i2]):InfPID;
01541     int PID3 = (i3<Len3)?(ExportPID3[i3]):InfPID;
01542 
01543     int_type GID1 = (i1<Len1)?((int_type) DomainMap_.GID64(ExportLID1[i1])):InfGID;
01544     int_type GID2 = (i2<Len2)?((int_type) DomainMap_.GID64(ExportLID2[i2])):InfGID;
01545     int_type GID3 = (i3<Len3)?((int_type) DomainMap_.GID64(ExportLID3[i3])):InfGID;
01546 
01547     int MIN_PID = MIN3(PID1,PID2,PID3);
01548     int_type MIN_GID = MIN3( ((PID1==MIN_PID)?GID1:InfGID), ((PID2==MIN_PID)?GID2:InfGID), ((PID3==MIN_PID)?GID3:InfGID));
01549     bool added_entry=false;
01550 
01551     // Case 1: Add off list 1   
01552     if(PID1 == MIN_PID && GID1 == MIN_GID){
01553       ExportLIDs[current] = ExportLID1[i1];
01554       ExportPIDs[current] = ExportPID1[i1];
01555       current++;
01556       i1++;
01557       added_entry=true;
01558     }
01559 
01560     // Case 2: Add off list 2
01561     if(PID2 == MIN_PID && GID2 == MIN_GID){
01562       if(!added_entry) {
01563   ExportLIDs[current] = ExportLID2[i2];
01564   ExportPIDs[current] = ExportPID2[i2];
01565   current++;
01566   added_entry=true;
01567       }
01568       i2++; 
01569     }
01570 
01571     // Case 3: Add off list 3
01572     if(PID3 == MIN_PID && GID3 == MIN_GID){
01573       if(!added_entry) {
01574   ExportLIDs[current] = ExportLID3[i3];
01575   ExportPIDs[current] = ExportPID3[i3];
01576   current++;
01577       }
01578       i3++; 
01579     }
01580   }// end while
01581   if(current!=MyLen) {
01582     ExportLIDs.resize(current);
01583     ExportPIDs.resize(current);
01584   }
01585 #endif
01586   return 0;
01587 }
01588 
01589 //=========================================================================
01590 template<typename ImportType, typename int_type>
01591 void LightweightCrsMatrix::Construct(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter)
01592 {  
01593   // Do we need to use long long for GCIDs?
01594 
01595 #ifdef ENABLE_MMM_TIMINGS
01596   using Teuchos::TimeMonitor;
01597   Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-1")));
01598 #endif
01599 
01600   // Fused constructor, import & FillComplete
01601   int rv=0;
01602   int N;
01603   if(use_lw) N = RowMapLW_->NumMyElements();
01604   else N = RowMapEP_->NumMyElements();
01605 
01606   int MyPID = SourceMatrix.Comm().MyPID();
01607 
01608 #ifdef HAVE_MPI
01609   std::vector<int> ReverseSendSizes;
01610   std::vector<int_type> ReverseSendBuffer;
01611   std::vector<int> ReverseRecvSizes;
01612   int_type * ReverseRecvBuffer=0;
01613 #endif
01614 
01615   bool communication_needed = RowImporter.SourceMap().DistributedGlobal();
01616 
01617   // The basic algorithm here is:
01618   // 1) Call Distor.Do to handle the import.
01619   // 2) Copy all the Imported and Copy/Permuted data into the raw Epetra_CrsMatrix / Epetra_CrsGraphData pointers, still using GIDs.
01620   // 3) Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids) AND
01621   //    reindexes to LIDs.
01622 
01623   // Sanity Check
01624   if (!SourceMatrix.RowMap().SameAs(RowImporter.SourceMap())) 
01625     throw "LightweightCrsMatrix: Fused copy constructor requires Importer.SourceMap() to match SourceMatrix.RowMap()";
01626 
01627   // Get information from the Importer
01628   int NumSameIDs             = RowImporter.NumSameIDs();
01629   int NumPermuteIDs          = RowImporter.NumPermuteIDs();
01630   int NumRemoteIDs           = RowImporter.NumRemoteIDs();
01631   int NumExportIDs           = RowImporter.NumExportIDs();
01632   int* ExportLIDs            = RowImporter.ExportLIDs();
01633   int* RemoteLIDs            = RowImporter.RemoteLIDs();
01634   int* PermuteToLIDs         = RowImporter.PermuteToLIDs();
01635   int* PermuteFromLIDs       = RowImporter.PermuteFromLIDs();
01636 
01637 #ifdef HAVE_MPI
01638   Epetra_Distributor& Distor            = RowImporter.Distributor();
01639   const Epetra_MpiComm * MpiComm        = dynamic_cast<const Epetra_MpiComm*>(&SourceMatrix.Comm());
01640   const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
01641 #endif
01642 
01643   // Allocate memory
01644   rowptr_.resize(N+1);
01645 
01646 
01647   // Owning PIDs
01648   std::vector<int> SourcePids;
01649   std::vector<int> TargetPids;
01650 
01651 
01652   /***************************************************/
01653   /***** 1) From Epetra_DistObject::DoTransfer() *****/
01654   /***************************************************/
01655   //  rv=SourceMatrix.CheckSizes(SourceMatrix);
01656 
01657   // NTS: Add CheckSizes stuff here.
01658   if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in CheckSizes()";
01659   
01660   // Buffers & Other Relevant Info
01661   char* Exports_  = 0;
01662   char* Imports_  = 0;
01663   int LenImports_ =0;
01664   int LenExports_ = 0;
01665   int *Sizes_     = 0;
01666 
01667   int SizeOfPacket; 
01668   bool VarSizes = false;
01669   if( NumExportIDs > 0) {
01670     Sizes_ = new int[NumExportIDs];
01671   }
01672   
01673 
01674   // Get the owning PIDs
01675   const Epetra_Import *MyImporter= SourceMatrix.Importer();
01676   if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,false);
01677   else {
01678     SourcePids.resize(SourceMatrix.ColMap().NumMyElements());
01679     SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),SourceMatrix.Comm().MyPID());
01680   }
01681   
01682   // Pack & Prepare w/ owning PIDs
01683   rv=Epetra_Import_Util::PackAndPrepareWithOwningPIDs(SourceMatrix,
01684                   NumExportIDs,ExportLIDs,
01685                   LenExports_,Exports_,SizeOfPacket,
01686                   Sizes_,VarSizes,SourcePids);
01687   if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in PackAndPrepare()";
01688 
01689   if (communication_needed) {
01690 #ifdef HAVE_MPI 
01691     // Do the exchange of remote data
01692     int i,curr_pid;
01693     const int * ExportPIDs = RowImporter.ExportPIDs();
01694 
01695     // Use the fact that the export procs are sorted to avoid building a hash table.
01696     // NOTE: The +1's on the message size lists are to avoid std::vector problems if a proc has no sends or recvs. 
01697     std::vector<int> SendSizes(MDistor->NumSends()+1,0);
01698     for(i=0, curr_pid=0; i<NumExportIDs; i++) { 
01699       if(i>0 &&  ExportPIDs[i] > ExportPIDs[i-1]) curr_pid++;
01700       SendSizes[curr_pid] +=Sizes_[i];
01701 
01702       // sanity check
01703       if(i>0 &&  ExportPIDs[i] < ExportPIDs[i-1]) throw "ExportPIDs not sorted";
01704     }
01705 
01706     std::vector<int> RecvSizes(MDistor->NumReceives()+1);
01707     int msg_tag=MpiComm->GetMpiTag();
01708     boundary_exchange_varsize<char>(*MpiComm,MPI_CHAR,MDistor->NumSends(),MDistor->ProcsTo(),&SendSizes[0],Exports_,
01709             MDistor->NumReceives(),MDistor->ProcsFrom(),&RecvSizes[0],Imports_,SizeOfPacket,msg_tag);
01710 
01711     // If the  source matrix doesn't have an importer, then nobody sent data belonging to me in the forward round.
01712     if(SourceMatrix.Importer()) {
01713       Epetra_Import* SourceImporter=const_cast<Epetra_Import*>(SourceMatrix.Importer());
01714       const Epetra_MpiDistributor * MyDistor = dynamic_cast<Epetra_MpiDistributor*>(&SourceImporter->Distributor());
01715 
01716       // Setup the reverse communication
01717       // Note: Buffer pairs are in (PID,GID) order
01718       PackAndPrepareReverseComm<ImportType, int_type>(SourceMatrix,RowImporter,ReverseSendSizes,ReverseSendBuffer);
01719       
01720       // Do the reverse communication
01721       // NOTE: Make the vector one too large to avoid std::vector errors
01722       ReverseRecvSizes.resize(MyDistor->NumSends()+1);
01723       int msg_tag=MpiComm->GetMpiTag();
01724     MPI_Datatype data_type = sizeof(int_type) == 4 ? MPI_INT : MPI_LONG_LONG;
01725       boundary_exchange_varsize<int_type>(*MpiComm,data_type,MyDistor->NumReceives(),MyDistor->ProcsFrom(),&ReverseSendSizes[0],&ReverseSendBuffer[0],
01726               MyDistor->NumSends(),MyDistor->ProcsTo(),&ReverseRecvSizes[0],ReverseRecvBuffer,1,msg_tag);      
01727     }
01728 #endif
01729   }
01730 
01731   if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in Distor.Do";
01732 
01733   /*********************************************************************/
01734   /**** 2) Copy all of the Same/Permute/Remote data into CSR_arrays ****/
01735   /*********************************************************************/
01736 #ifdef ENABLE_MMM_TIMINGS
01737   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-2")));
01738 #endif
01739 
01740   // Count nnz
01741   int mynnz = Epetra_Import_Util::UnpackWithOwningPIDsCount(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_);
01742   // Presume the rowptr_ is the right size
01743   // Allocate colind_ & vals_ arrays
01744   colind_.resize(mynnz);
01745 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01746   if(sizeof(int_type) == sizeof(long long))
01747     colind_LL_.resize(mynnz);
01748 #endif
01749   vals_.resize(mynnz);
01750 
01751   // Reset the Source PIDs (now with -1 rule)
01752   if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,true);
01753   else {
01754     SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),-1);
01755   }
01756    
01757   // Unpack into arrays
01758   int * myrowptr  = rowptr_.size() ? & rowptr_[0] : 0;
01759   double * myvals = vals_.size()   ? & vals_[0] : 0;
01760 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01761   if(sizeof(int_type) == sizeof(int)) {
01762     int * mycolind  = colind_.size() ? & colind_[0] : 0;
01763     Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
01764   }
01765   else
01766 #endif
01767 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01768   if(sizeof(int_type) == sizeof(long long)) {
01769     long long * mycolind  = colind_LL_.size() ? & colind_LL_[0] : 0;
01770     Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
01771   }
01772   else
01773 #endif
01774   throw "EpetraExt::LightweightCrsMatrix::Construct: sizeof(int_type) error.";
01775 
01776   /**************************************************************/
01777   /**** 3) Call Optimized MakeColMap w/ no Directory Lookups ****/
01778   /**************************************************************/
01779 #ifdef ENABLE_MMM_TIMINGS
01780   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-3")));
01781 #endif
01782 
01783   //Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids).
01784   MakeColMapAndReindex<int_type>(TargetPids,getcolind<int_type>());
01785 
01786   /********************************************/
01787   /**** 4) Make Export Lists for Import    ****/
01788   /********************************************/
01789 #ifdef HAVE_MPI
01790   MakeExportLists<ImportType, int_type>(SourceMatrix,RowImporter,ReverseRecvSizes,ReverseRecvBuffer,ExportPIDs_,ExportLIDs_);
01791 #endif 
01792 
01793   /********************************************/
01794   /**** 5) Call sort the entries           ****/
01795   /********************************************/
01796   // NOTE: If we have no entries the &blah[0] will cause the STL to die in debug mode
01797 #ifdef ENABLE_MMM_TIMINGS
01798   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-4")));
01799 #endif
01800   if(N>0) Epetra_Util::SortCrsEntries(N, &rowptr_[0], &colind_[0], &vals_[0]);
01801 
01802   /********************************************/
01803   /**** 6) Cleanup                         ****/
01804   /********************************************/
01805 #ifdef ENABLE_MMM_TIMINGS
01806   MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-5")));
01807 #endif
01808 
01809 #ifdef HAVE_MPI
01810   delete [] ReverseRecvBuffer;
01811 #endif
01812 
01813   delete [] Exports_;
01814   delete [] Imports_;
01815   delete [] Sizes_;
01816 
01817  }// end fused copy constructor
01818 
01819 
01820 
01821 
01822 //=========================================================================
01823 LightweightCrsMatrix::LightweightCrsMatrix(const Epetra_CrsMatrix & SourceMatrix, RemoteOnlyImport & RowImporter):
01824   use_lw(true),
01825   RowMapLW_(0),
01826   RowMapEP_(0),
01827   DomainMap_(SourceMatrix.DomainMap())
01828 { 
01829 #ifdef ENABLE_MMM_TIMINGS
01830   using Teuchos::TimeMonitor;
01831   Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS Total")));
01832 #endif
01833 
01834   RowMapLW_= new LightweightMap(RowImporter.TargetMap());
01835    
01836 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01837   if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
01838     Construct<RemoteOnlyImport, int>(SourceMatrix,RowImporter);
01839   }
01840   else
01841 #endif
01842 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01843   if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
01844     Construct<RemoteOnlyImport, long long>(SourceMatrix,RowImporter);
01845   }
01846   else
01847 #endif
01848     throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
01849 
01850 }
01851 
01852 
01853 //=========================================================================
01854 LightweightCrsMatrix::LightweightCrsMatrix(const Epetra_CrsMatrix & SourceMatrix, Epetra_Import & RowImporter):
01855   use_lw(false),
01856   RowMapLW_(0),
01857   RowMapEP_(0),
01858   DomainMap_(SourceMatrix.DomainMap())
01859 {
01860   RowMapEP_= new Epetra_BlockMap(RowImporter.TargetMap());
01861 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01862   if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
01863     Construct<Epetra_Import, int>(SourceMatrix,RowImporter);
01864   }
01865   else
01866 #endif
01867 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01868   if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
01869     Construct<Epetra_Import, long long>(SourceMatrix,RowImporter);
01870   }
01871   else
01872 #endif
01873     throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
01874 }
01875 
01876 
01877 LightweightCrsMatrix::~LightweightCrsMatrix(){
01878   delete RowMapLW_;
01879   delete RowMapEP_;
01880 }
01881 
01882 
01883 }//namespace EpetraExt
01884 
01885 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01886 template class EpetraExt::CrsWrapper_GraphBuilder<int>;
01887 #endif
01888 
01889 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01890 template class EpetraExt::CrsWrapper_GraphBuilder<long long>;
01891 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines