44 #include <Epetra_Comm.h> 45 #include <Epetra_CrsMatrix.h> 46 #include <Epetra_Import.h> 47 #include <Epetra_Export.h> 48 #include <Epetra_Distributor.h> 49 #include <Epetra_HashTable.h> 50 #include <Epetra_Util.h> 51 #include <Epetra_Import_Util.h> 52 #include <Epetra_GIDTypeSerialDenseVector.h> 54 #include <Teuchos_TimeMonitor.hpp> 59 #include "Epetra_MpiComm.h" 60 #include "Epetra_MpiDistributor.h" 62 #define MIN(x,y) ((x)<(y)?(x):(y)) 63 #define MIN3(x,y,z) ((x)<(y)?(MIN(x,z)):(MIN(y,z))) 69 void debug_print_distor(
const char * label,
const Epetra_Distributor * Distor,
const Epetra_Comm & Comm) {
71 const Epetra_MpiDistributor * MDistor =
dynamic_cast<const Epetra_MpiDistributor*
>(Distor);
72 printf(
"[%d] %s\n",Comm.MyPID(),label);
73 printf(
"[%d] NumSends = %d NumRecvs = %d\n",Comm.MyPID(),MDistor->NumSends(),MDistor->NumReceives());
74 printf(
"[%d] ProcsTo = ",Comm.MyPID());
75 for(
int ii=0; ii<MDistor->NumSends(); ii++)
76 printf(
"%d ",MDistor->ProcsTo()[ii]);
77 printf(
"\n[%d] ProcsFrom = ",Comm.MyPID());
78 for(
int ii=0; ii<MDistor->NumReceives(); ii++)
79 printf(
"%d ",MDistor->ProcsFrom()[ii]);
88 if(!Import1 && !Import2)
return;
89 const Epetra_Comm & Comm = (Import1)? Import1->SourceMap().Comm() : Import2->SourceMap().Comm();
92 if( (!Import1 && Import2) || (Import2 && !Import1) ) {printf(
"[%d] DCI: One Import exists, the other does not\n",PID);
return;}
93 if(!Import1->SourceMap().SameAs(Import2->SourceMap())) {printf(
"[%d] DCI: SourceMaps don't match\n",PID);
return;}
94 if(!Import1->TargetMap().SameAs(Import2->TargetMap())) {printf(
"[%d] DCI: TargetMaps don't match\n",PID);
return;}
96 if(Import1->NumSameIDs() != Import2->NumSameIDs()) {printf(
"[%d] DCI NumSameIDs() mismatch %d vs. %d\n",PID,Import1->NumSameIDs(),Import2->NumSameIDs());flag=
false;}
98 if(Import1->NumPermuteIDs() != Import2->NumPermuteIDs()) {printf(
"[%d] DCI NumPermuteIDs() mismatch %d vs. %d\n",PID,Import1->NumPermuteIDs(),Import2->NumPermuteIDs()); flag=
false;}
100 if(Import1->NumExportIDs() != Import2->NumExportIDs()) {printf(
"[%d] DCI NumExportIDs() mismatch %d vs. %d\n",PID,Import1->NumExportIDs(),Import2->NumExportIDs()); flag=
false;}
102 if(Import1->NumRemoteIDs() != Import2->NumRemoteIDs()) {printf(
"[%d] DCI NumRemoteIDs() mismatch %d vs. %d\n",PID,Import1->NumRemoteIDs(),Import2->NumRemoteIDs()); flag=
false;}
104 if(Import1->NumSend() != Import2->NumSend()) {printf(
"[%d] DCI NumSend() mismatch %d vs. %d\n",PID,Import1->NumSend(),Import2->NumSend()); flag=
false;}
106 if(Import1->NumRecv() != Import2->NumRecv()) {printf(
"[%d] DCI NumRecv() mismatch %d vs. %d\n",PID,Import1->NumRecv(),Import2->NumRecv()); flag=
false;}
109 if(flag) printf(
"[%d] DCI Importers compare OK\n",PID);
111 Import1->SourceMap().Comm().Barrier();
112 Import1->SourceMap().Comm().Barrier();
113 Import1->SourceMap().Comm().Barrier();
121 : numRows(0), numEntriesPerRow(NULL), indices(NULL), values(NULL),
122 remote(NULL), numRemote(0), importColMap(NULL), rowMap(NULL), colMap(NULL),
123 domainMap(NULL), importMatrix(NULL), origMatrix(NULL)
150 std::cout <<
"proc " << M.
rowMap->Comm().MyPID()<<std::endl;
151 std::cout <<
"numRows: " << M.
numRows<<std::endl;
152 for(
int i=0; i<M.
numRows; ++i) {
155 std::cout <<
" *"<<M.
rowMap->GID64(i)<<
" " 159 std::cout <<
" "<<M.
rowMap->GID64(i)<<
" " 168 : ecrsmat_(epetracrsmatrix)
187 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 191 return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
197 return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
201 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 205 return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
211 return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
218 template<
typename int_type>
224 int num_rows = emap.NumMyElements();
226 emap.MyGlobalElementsPtr(rows);
228 for(
int i=0; i<num_rows; ++i) {
229 graph_[rows[i]] =
new std::set<int_type>;
233 template<
typename int_type>
236 typename std::map<int_type,std::set<int_type>*>::iterator
237 iter = graph_.begin(), iter_end = graph_.end();
238 for(; iter!=iter_end; ++iter) {
245 template<
typename int_type>
251 template<
typename int_type>
255 typename std::map<int_type,std::set<int_type>*>::iterator
256 iter = graph_.find(GlobalRow);
258 if (iter == graph_.end())
return(-1);
260 std::set<int_type>& cols = *(iter->second);
262 for(
int i=0; i<NumEntries; ++i) {
263 cols.insert(Indices[i]);
266 int row_length = cols.size();
267 if (row_length > max_row_length_) max_row_length_ = row_length;
272 template<
typename int_type>
276 return InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
279 template<
typename int_type>
280 std::map<int_type,std::set<int_type>*>&
286 template<
typename int_type>
291 if (max_row_length < 1)
return;
293 std::vector<int_type> indices(max_row_length);
294 int_type* indices_ptr = &indices[0];
295 std::vector<double> zeros(max_row_length, 0.0);
296 double* zeros_ptr = &zeros[0];
298 std::map<int_type,std::set<int_type>*>& graph = graphbuilder.
get_graph();
300 typename std::map<int_type,std::set<int_type>*>::iterator
301 iter = graph.begin(), iter_end = graph.end();
303 for(; iter!=iter_end; ++iter) {
304 int_type row = iter->first;
305 std::set<int_type>& cols = *(iter->second);
306 int num_entries = cols.size();
308 typename std::set<int_type>::iterator
309 col_iter = cols.begin(), col_end = cols.end();
310 for(
int j=0; col_iter!=col_end; ++col_iter, ++j) {
311 indices_ptr[j] = *col_iter;
314 C.InsertGlobalValues(row, num_entries, zeros_ptr, indices_ptr);
318 template<
typename int_type>
320 const std::vector<int_type>& proc_col_ranges,
321 std::vector<int_type>& send_rows,
322 std::vector<int>& rows_per_send_proc)
324 const Epetra_Map& rowmap = mtx.RowMap();
325 int numrows = mtx.NumMyRows();
326 const Epetra_CrsGraph& graph = mtx.Graph();
328 int* col_indices = NULL;
329 int_type* Tcol_indices = NULL;
330 int num_col_ranges = proc_col_ranges.size()/2;
331 rows_per_send_proc.resize(num_col_ranges);
333 for(
int nc=0; nc<num_col_ranges; ++nc) {
334 int_type first_col = proc_col_ranges[nc*2];
335 int_type last_col = proc_col_ranges[nc*2+1];
336 int num_send_rows = 0;
337 for(
int i=0; i<numrows; ++i) {
338 int_type grow = (int_type) rowmap.GID64(i);
340 const Epetra_Map& colmap = mtx.ColMap();
341 graph.ExtractMyRowView(i, rowlen, col_indices);
342 for(
int j=0; j<rowlen; ++j) {
343 int_type col = (int_type) colmap.GID64(col_indices[j]);
344 if (first_col <= col && last_col >= col) {
346 send_rows.push_back(grow);
352 graph.ExtractGlobalRowView(grow, rowlen, Tcol_indices);
353 for(
int j=0; j<rowlen; ++j) {
354 if (first_col <= Tcol_indices[j] && last_col >= Tcol_indices[j]) {
356 send_rows.push_back(grow);
362 rows_per_send_proc[nc] = num_send_rows;
366 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 368 const std::vector<int>& proc_col_ranges,
369 std::vector<int>& send_rows,
370 std::vector<int>& rows_per_send_proc)
372 if(mtx.RowMatrixRowMap().GlobalIndicesInt()) {
373 Tpack_outgoing_rows<int>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
376 throw "EpetraExt::pack_outgoing_rows: Global indices not int";
381 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 383 const std::vector<long long>& proc_col_ranges,
384 std::vector<long long>& send_rows,
385 std::vector<int>& rows_per_send_proc)
387 if(mtx.RowMatrixRowMap().GlobalIndicesLongLong()) {
388 Tpack_outgoing_rows<long long>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
391 throw "EpetraExt::pack_outgoing_rows: Global indices not long long";
396 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 400 if(emap.GlobalIndicesInt())
401 return std::make_pair(emap.MinMyGID(),emap.MaxMyGID());
403 throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_Map";
407 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 411 if(emap.GlobalIndicesLongLong())
412 return std::make_pair(emap.MinMyGID64(),emap.MaxMyGID64());
414 throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_Map";
418 template<
typename int_type>
421 std::pair<int_type,int_type> col_range;
423 col_range = get_col_range<int_type>(mtx.ColMap());
426 const Epetra_Map& row_map = mtx.RowMap();
427 col_range.first = row_map.MaxMyGID64();
428 col_range.second = row_map.MinMyGID64();
430 int_type* col_indices = NULL;
431 const Epetra_CrsGraph& graph = mtx.Graph();
432 for(
int i=0; i<row_map.NumMyElements(); ++i) {
433 graph.ExtractGlobalRowView((int_type) row_map.GID64(i), rowlen, col_indices);
434 for(
int j=0; j<rowlen; ++j) {
435 if (col_indices[j] < col_range.first) col_range.first = col_indices[j];
436 if (col_indices[j] > col_range.second) col_range.second = col_indices[j];
444 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 448 if(mtx.RowMatrixRowMap().GlobalIndicesInt())
449 return Tget_col_range<int>(mtx);
451 throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_CrsMatrix";
455 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 459 if(mtx.RowMatrixRowMap().GlobalIndicesLongLong())
460 return Tget_col_range<long long>(mtx);
462 throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_CrsMatrix";
471 template <
typename MyType>
472 void boundary_exchange(
const Epetra_MpiComm Comm, MPI_Datatype DataType,
473 int NumSends,
const int * SendProcs,
const int * SendSizes, MyType* SendBuffer,
474 int NumRecvs,
const int * RecvProcs,
const int * RecvSizes, MyType* RecvBuffer,
int SizeOfPacket,
int msg_tag)
477 MPI_Comm comm = Comm.Comm();
478 std::vector<MPI_Request> requests(NumRecvs);
479 std::vector<MPI_Status> status(NumRecvs);
481 int i,num_waits=0,MyPID=Comm.MyPID();
482 int start, self_recv_len=-1,self_recv_start=-1, self_send_start=-1;
485 int mysendsize=1, myrecvsize=1;
489 for(i=0; i<NumRecvs; i++){
490 if(RecvSizes) myrecvsize=RecvSizes[i]*SizeOfPacket;
491 if(RecvProcs[i] != MyPID) {
492 MPI_Irecv(RecvBuffer + start, myrecvsize, DataType, RecvProcs[i], msg_tag, comm, &requests[num_waits]);
496 self_recv_len = myrecvsize;
497 self_recv_start=start;
504 for(i=0; i<NumSends; i++){
505 if(SendSizes) mysendsize=SendSizes[i]*SizeOfPacket;
506 if(SendProcs[i] != MyPID)
507 MPI_Send(SendBuffer + start, mysendsize,DataType,SendProcs[i],msg_tag,comm);
509 self_send_start=start;
514 if(self_recv_len != -1)
515 memcpy(RecvBuffer+self_recv_start,SendBuffer+self_send_start,self_recv_len*
sizeof(MyType)*SizeOfPacket);
519 MPI_Waitall(num_waits, &requests[0],&status[0]);
525 template <
typename MyType>
526 void boundary_exchange_varsize(
const Epetra_MpiComm Comm, MPI_Datatype DataType,
527 int NumSends,
const int * SendProcs,
const int * SendSizes, MyType* SendBuffer,
528 int NumRecvs,
const int * RecvProcs,
int * RecvSizes, MyType*& RecvBuffer,
int SizeOfPacket,
int msg_tag)
534 boundary_exchange<int>(Comm,MPI_INT,NumSends,SendProcs,(
int*)0,const_cast<int*>(SendSizes),NumRecvs,RecvProcs,(
int*)0,RecvSizes,1,msg_tag);
537 for(i=0; i<NumRecvs; i++) rbuffersize+=RecvSizes[i]*SizeOfPacket;
538 RecvBuffer =
new MyType[rbuffersize];
541 boundary_exchange<MyType>(Comm,DataType,NumSends,SendProcs,SendSizes,SendBuffer,NumRecvs,RecvProcs,RecvSizes,RecvBuffer,SizeOfPacket,msg_tag+100);
552 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
555 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
563 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 566 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 576 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 583 if(GenerateHash)
Data_->
LIDHash_int_ =
new Epetra_HashTable<int>(numMyElements + 1 );
584 for(
int i=0; i < numMyElements; ++i ) {
593 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 600 if(GenerateHash)
Data_->
LIDHash_LL_ =
new Epetra_HashTable<long long>(numMyElements + 1 );
601 for(
int i=0; i < numMyElements; ++i ) {
610 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 613 Construct_int(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
617 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 620 Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
625 Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
635 IsInt = Map.GlobalIndicesInt();
642 Data_->IncrementReferenceCount();
658 Data_->IncrementReferenceCount();
668 Data_->DecrementReferenceCount();
669 if(
Data_->ReferenceCount() == 0) {
675 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 680 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 689 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 694 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 699 throw "EpetraExt::LightweightMap::NumMyElements: Global indices unknowns";
703 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 709 throw "EpetraExt::LightweightMap::LID: Int version called for long long map";
711 throw "EpetraExt::LightweightMap::LID: unknown GID type";
714 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 721 throw "EpetraExt::LightweightMap::LID: Long long version called for int map";
723 throw "EpetraExt::LightweightMap::LID: unknown GID type";
728 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 737 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 744 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 751 throw "EpetraExt::LightweightMap::GID64: Global indices unknown.";
755 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 762 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 780 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 785 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 790 throw "EpetraExt::LightweightMap::MaxLID: Global indices unknowns";
810 Distor_ = &Importer.Distributor();
811 int * OldRemoteLIDs = Importer.RemoteLIDs();
812 int * OldExportLIDs = Importer.ExportLIDs();
813 int * OldExportPIDs = Importer.ExportPIDs();
817 throw std::runtime_error(
"RemoteOnlyImport: Importer doesn't match RemoteOnlyTargetMap for number of remotes.");
840 throw std::runtime_error(
"RemoteOnlyImport: TargetMap_ index type unknown.");
845 throw std::runtime_error(
"RemoteOnlyImport: Importer and RemoteOnlyTargetMap order don't match.");
863 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs,
bool SortGhostsAssociatedWithEachProcessor_);
867 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs,
bool SortGhostsAssociatedWithEachProcessor_)
872 if(NumRemoteColGIDs > 0) {
873 SortLists[0] = RemoteColindices;
874 SortLists[1] = &RemotePermuteIDs[0];
875 Epetra_Util::Sort(
true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, NLists, SortLists);
880 if (SortGhostsAssociatedWithEachProcessor_) {
886 int StartCurrent, StartNext;
887 StartCurrent = 0; StartNext = 1;
888 while ( StartNext < NumRemoteColGIDs ) {
889 if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
891 SortLists[0] = &RemotePermuteIDs[StartCurrent];
892 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists);
893 StartCurrent = StartNext; StartNext++;
896 SortLists[0] = &RemotePermuteIDs[StartCurrent];
897 Epetra_Util::Sort(
true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists);
903 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs,
bool SortGhostsAssociatedWithEachProcessor_)
906 if(NumRemoteColGIDs > 0) {
907 long long* SortLists_LL[1] = {RemoteColindices};
908 int* SortLists_int[1] = {&RemotePermuteIDs[0]};
909 Epetra_Util::Sort(
true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, 1, SortLists_int, 1, SortLists_LL);
913 if (SortGhostsAssociatedWithEachProcessor_) {
919 int StartCurrent, StartNext;
920 StartCurrent = 0; StartNext = 1;
921 while ( StartNext < NumRemoteColGIDs ) {
922 if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
924 int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
925 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
926 StartCurrent = StartNext; StartNext++;
929 int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
930 Epetra_Util::Sort(
true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
937 #ifdef ENABLE_MMM_TIMINGS 939 if(label) tpref = std::string(label);
940 using Teuchos::TimeMonitor;
941 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS-3.1"))));
947 int numDomainElements =
DomainMap_.NumMyElements();
948 std::vector<bool> LocalGIDs(numDomainElements,
false);
950 bool DoSizes = !
DomainMap_.ConstantElementSize();
951 if(DoSizes) EPETRA_CHK_ERR(-1);
958 else numMyBlockRows =
RowMapEP_->NumMyElements();
960 int hashsize = numMyBlockRows;
if (hashsize < 100) hashsize = 100;
961 Epetra_HashTable<GO> RemoteGIDs(hashsize);
962 std::vector<GO> RemoteGIDList; RemoteGIDList.reserve(hashsize);
963 std::vector<int> RemoteOwningPIDs; RemoteOwningPIDs.reserve(hashsize);
968 int NumLocalColGIDs = 0;
969 int NumRemoteColGIDs = 0;
970 for(
int i = 0; i < numMyBlockRows; i++) {
976 bool alreadyFound = LocalGIDs[LID];
978 LocalGIDs[LID] =
true;
984 GO hash_value=RemoteGIDs.Get(GID);
985 if(hash_value == -1) {
986 int PID = owningPIDs[j];
987 if(PID==-1) printf(
"[%d] ERROR: Remote PID should not be -1\n",
DomainMap_.Comm().MyPID());
988 colind_[j] = numDomainElements + NumRemoteColGIDs;
989 RemoteGIDs.Add(GID, NumRemoteColGIDs);
990 RemoteGIDList.push_back(GID);
991 RemoteOwningPIDs.push_back(PID);
995 colind_[j] = numDomainElements + hash_value;
1002 if (NumRemoteColGIDs!=0) {
1003 throw "Some column IDs are not in domainMap. If matrix is rectangular, you must pass in domainMap to FillComplete";
1006 if (NumLocalColGIDs==numDomainElements) {
1017 int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1018 typename Epetra_GIDTypeSerialDenseVector<GO>::impl Colindices;
1019 if(numMyBlockCols > 0)
1020 Colindices.Size(numMyBlockCols);
1021 GO* RemoteColindices = Colindices.Values() + NumLocalColGIDs;
1023 for(
int i = 0; i < NumRemoteColGIDs; i++)
1024 RemoteColindices[i] = RemoteGIDList[i];
1027 std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
1028 for(
int i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
1030 MakeColMapAndReindexSort<GO>(NumRemoteColGIDs, RemoteColindices, RemotePermuteIDs, RemoteOwningPIDs,SortGhosts);
1033 std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
1034 for(
int i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
1037 bool use_local_permute=
false;
1038 std::vector<int> LocalPermuteIDs(numDomainElements);
1046 if(NumLocalColGIDs ==
DomainMap_.NumMyElements()) {
1047 DomainMap_.MyGlobalElements(Colindices.Values());
1050 GO* MyGlobalElements = 0;
1051 DomainMap_.MyGlobalElementsPtr(MyGlobalElements);
1052 int NumLocalAgain = 0;
1053 use_local_permute =
true;
1054 for(
int i = 0; i < numDomainElements; i++) {
1056 LocalPermuteIDs[i] = NumLocalAgain;
1057 Colindices[NumLocalAgain++] = MyGlobalElements[i];
1060 assert(NumLocalAgain==NumLocalColGIDs);
1067 for(
int i=0;i<NumRemoteColGIDs;i++)
1070 #ifdef ENABLE_MMM_TIMINGS 1071 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS-3.2"))));
1079 #ifdef ENABLE_MMM_TIMINGS 1080 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS-3.3"))));
1084 for(
int i=0; i<numMyBlockRows; i++){
1087 if(ID < numDomainElements){
1093 colind_[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[
colind_[j]-numDomainElements];
1107 static inline bool lessthan12(std::pair<int,int> i, std::pair<int,int> j){
1108 return ((i.first<j.first) || (i.first==j.first && i.second <j.second));
1113 template<
typename ImportType,
typename int_type>
1115 std::vector<int> &ReverseSendSizes, std::vector<int_type> &ReverseSendBuffer) {
1119 const Epetra_Import *MyImporter= SourceMatrix.Importer();
1120 if(MyImporter == 0)
return -1;
1121 const Epetra_MpiComm * MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&SourceMatrix.Comm());
1122 int MyPID = MpiComm->MyPID();
1125 int NumExportIDs = RowImporter.NumExportIDs();
1126 int* ExportLIDs = RowImporter.ExportLIDs();
1127 int* ExportPIDs = RowImporter.ExportPIDs();
1130 Epetra_Distributor& Distor = MyImporter->Distributor();
1131 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1132 int NumRecvs = MDistor->NumReceives();
1133 int* RemoteLIDs = MyImporter->RemoteLIDs();
1134 const int * ProcsFrom = MDistor->ProcsFrom();
1135 const int * LengthsFrom = MDistor->LengthsFrom();
1140 std::vector<int> RemotePIDOrder(SourceMatrix.NumMyCols(),-1);
1143 for(i=0,j=0;i<NumRecvs;i++){
1144 for(k=0;k<LengthsFrom[i];k++){
1145 int pid=ProcsFrom[i];
1146 if(pid!=MyPID) RemotePIDOrder[RemoteLIDs[j]]=i;
1152 std::vector<std::set<std::pair<int,int_type> > > ReversePGIDs(NumRecvs);
1153 int *rowptr, *colind;
1155 EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1159 for(i=0; i < NumExportIDs; i++) {
1160 int lid = ExportLIDs[i];
1161 int exp_pid = ExportPIDs[i];
1162 for(j=rowptr[lid]; j<rowptr[lid+1]; j++){
1163 int pid_order = RemotePIDOrder[colind[j]];
1165 int_type gid = (int_type) SourceMatrix.GCID64(colind[j]);
1167 ReversePGIDs[pid_order].insert(std::pair<int,int_type>(exp_pid,gid));
1173 ReverseSendSizes.resize(NumRecvs+1);
1175 for(i=0; i<NumRecvs; i++) {
1176 ReverseSendSizes[i] = 2*ReversePGIDs[i].size();
1177 totalsize += ReverseSendSizes[i];
1181 ReverseSendBuffer.resize(totalsize+1);
1182 for(i=0, j=0; i<NumRecvs; i++) {
1183 for(
typename std::set<std::pair<int,int_type> >::iterator it=ReversePGIDs[i].begin(); it!=ReversePGIDs[i].end(); it++) {
1184 ReverseSendBuffer[j] = it->first;
1185 ReverseSendBuffer[j+1] = it->second;
1196 template<
typename int_type>
1202 int * companion = &ExportGID3[0];
1203 Epetra_Util::Sort(
true,total_length3,&ExportPID3[0],0,0,1,&companion,0,0);
1209 long long * companion = &ExportGID3[0];
1210 Epetra_Util::Sort(
true,total_length3,&ExportPID3[0],0,0,0,0,1,&companion);
1213 template<
typename int_type>
1214 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){
1218 int total_length3=0;
1219 for(i=0; i<Nrecv; i++)
1220 total_length3+=ReverseRecvSizes[i]/2;
1221 if(total_length3==0)
return 0;
1223 std::vector<int_type> ExportGID3(total_length3);
1224 ExportLID3.resize(total_length3);
1225 ExportPID3.resize(total_length3);
1228 for(i=0,j=0; i<2*total_length3; i+=2) {
1229 if(ReverseRecvBuffer[i] != MyPID){
1230 ExportPID3[j]=ReverseRecvBuffer[i];
1231 ExportGID3[j]=ReverseRecvBuffer[i+1];
1237 if(total_length3==0)
return 0;
1240 build_type3_exports_sort<int_type>(ExportGID3, ExportPID3, total_length3);
1241 int StartCurrent, StartNext;
1242 StartCurrent = 0; StartNext = 1;
1243 while ( StartNext < total_length3 ) {
1244 if(ExportPID3[StartNext] == ExportPID3[StartNext-1]) StartNext++;
1246 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1247 StartCurrent = StartNext; StartNext++;
1250 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1261 for(i=1,j=1; i<total_length3; i++){
1262 if(ExportPID3[i]!=ExportPID3[i-1] || ExportGID3[i]!=ExportGID3[i-1]){
1263 ExportPID3[j] = ExportPID3[i];
1264 ExportGID3[j] = ExportGID3[i];
1268 ExportPID3.resize(j);
1269 ExportLID3.resize(j);
1281 for(i=0; i<total_length3; i++) {
1282 ExportLID3[i]=DomainMap.LID(ExportGID3[i]);
1283 if(ExportLID3[i] < 0)
throw std::runtime_error(
"LightweightCrsMatrix:MakeExportLists invalid LID");
1292 return total_length3;
1296 template<
typename ImportType,
typename int_type>
1297 int build_type2_exports(
const Epetra_CrsMatrix & SourceMatrix, ImportType & MyImporter, std::vector<int> &ExportLID2, std::vector<int> &ExportPID2){
1298 int total_length2=0;
1301 const Epetra_Import *SourceImporter= SourceMatrix.Importer();
1303 int *rowptr, *colind;
1305 EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1308 int NumExportIDs = MyImporter.NumExportIDs();
1309 const int* ExportLIDs = MyImporter.ExportLIDs();
1310 const int* ExportPIDs = MyImporter.ExportPIDs();
1311 if(NumExportIDs==0)
return 0;
1313 Epetra_Distributor& Distor = MyImporter.Distributor();
1314 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1318 std::vector<bool> IsOwned(SourceMatrix.NumMyCols(),
true);
1319 if(SourceImporter) {
1320 const int * RemoteLIDs = SourceImporter->RemoteLIDs();
1322 for(i=0; i<SourceImporter->NumRemoteIDs(); i++)
1323 IsOwned[RemoteLIDs[i]]=
false;
1327 std::vector<int> SentTo(SourceMatrix.NumMyCols(),-1);
1330 for(i=0,total_length2=0; i<MDistor->NumSends(); i++) total_length2+= MDistor->LengthsTo()[i] * SourceMatrix.MaxNumEntries();
1331 std::vector<int_type> ExportGID2(total_length2);
1333 ExportLID2.resize(total_length2);
1334 ExportPID2.resize(total_length2);
1336 int current=0, last_start=0, last_pid=ExportPIDs[0];
1337 for(i=0; i<NumExportIDs; i++){
1339 int row=ExportLIDs[i];
1340 int pid=ExportPIDs[i];
1342 if(i!=0 && pid>last_pid) {
1344 if(current!=last_start){
1345 int *lids = &ExportLID2[last_start];
1346 Epetra_Util::Sort(
true,current-last_start,&ExportGID2[last_start],0,0,1,&lids,0,0);
1353 else if(pid < last_pid) {
1354 throw std::runtime_error(
"build_type2_exports: ExportPIDs are not sorted!");
1357 for(j=rowptr[row]; j<rowptr[row+1]; j++) {
1360 if(IsOwned[col] && SentTo[col]!=pid){
1362 if(current>= total_length2)
throw std::runtime_error(
"build_type2_exports: More export ids than I thought!");
1364 ExportGID2[current] = (int_type) SourceMatrix.GCID64(col);
1365 ExportLID2[current] = SourceMatrix.DomainMap().LID(ExportGID2[current]);
1366 ExportPID2[current] = pid;
1373 int *lids = ExportLID2.size() > (std::size_t) last_start ? &ExportLID2[last_start] : 0;
1374 Epetra_Util::Sort(
true,current-last_start,ExportGID2.size() > (std::size_t) last_start ? &ExportGID2[last_start] : 0,0,0,1,&lids,0,0);
1377 total_length2=current;
1378 ExportLID2.resize(total_length2);
1379 ExportPID2.resize(total_length2);
1381 return total_length2;
1385 template<
typename int_type>
1386 void build_type1_exports_sort(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int_type>& ExportGID1,
int total_length1);
1390 int * companion[2] = {&ExportLID1[0],&ExportGID1[0]};
1391 Epetra_Util::Sort(
true,total_length1,&ExportPID1[0],0,0,2,&companion[0],0,0);
1396 int * companion = &ExportLID1[0];
1397 long long * companion64 = &ExportGID1[0];
1398 Epetra_Util::Sort(
true,total_length1,&ExportPID1[0],0,0,1,&companion,1,&companion64);
1401 template<
typename int_type>
1402 int build_type1_exports(
const Epetra_Import * Importer1, std::vector<int> &ExportLID1, std::vector<int> &ExportPID1){
1403 int i, total_length1=0;
1404 if(!Importer1)
return 0;
1405 total_length1 = Importer1->NumSend();
1406 if(total_length1==0)
return 0;
1408 std::vector<int_type> ExportGID1(total_length1);
1409 ExportLID1.resize(total_length1);
1410 ExportPID1.resize(total_length1);
1411 const int * ExportLID1Base = Importer1->ExportLIDs();
1412 const int * ExportPID1Base = Importer1->ExportPIDs();
1414 for(i=0; i<total_length1; i++){
1415 ExportLID1[i] = ExportLID1Base[i];
1416 ExportPID1[i] = ExportPID1Base[i];
1417 ExportGID1[i] = (int_type) Importer1->SourceMap().GID64(ExportLID1Base[i]);
1421 build_type1_exports_sort<int_type>(ExportLID1, ExportPID1, ExportGID1, total_length1);
1423 int StartCurrent, StartNext;
1424 StartCurrent = 0; StartNext = 1;
1425 while ( StartNext < total_length1 ) {
1426 if(ExportPID1[StartNext] == ExportPID1[StartNext-1]) StartNext++;
1428 int *new_companion = {&ExportLID1[StartCurrent]};
1429 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1430 StartCurrent = StartNext; StartNext++;
1433 int *new_companion = {&ExportLID1[StartCurrent]};
1434 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1435 return total_length1;
1439 template<
typename ImportType,
typename int_type>
1441 std::vector<int> &ReverseRecvSizes,
const int_type *ReverseRecvBuffer,
1442 std::vector<int> & ExportPIDs, std::vector<int> & ExportLIDs) {
1444 int MyPID = SourceMatrix.Comm().MyPID();
1447 const Epetra_Import *Importer1= SourceMatrix.Importer();
1457 Epetra_MpiDistributor * Distor1 = (Importer1)?(dynamic_cast<Epetra_MpiDistributor*>(&Importer1->Distributor())):0;
1459 int Nsend1 = (Distor1)?(Distor1->NumSends()):0;
1461 std::vector<int> ExportPID3;
1462 std::vector<int> ExportLID3;
1464 std::vector<int> ExportPID2;
1465 std::vector<int> ExportLID2;
1467 std::vector<int> ExportPID1;
1468 std::vector<int> ExportLID1;
1471 int Len1=build_type1_exports<int_type>(Importer1, ExportLID1, ExportPID1);
1472 int Len2=build_type2_exports<ImportType, int_type>(SourceMatrix, Importer2, ExportLID2, ExportPID2);
1476 #ifdef HAVE_EPETRAEXT_DEBUG 1480 bool test_passed=
true;
1481 for(i=1; i<Len1; i++) {
1482 if(ExportPID1[i] < ExportPID1[i-1] || (ExportPID1[i] == ExportPID1[i-1] &&
DomainMap_.GID(ExportLID1[i]) <
DomainMap_.GID(ExportLID1[i-1])))
1485 SourceMatrix.Comm().Barrier();
1487 printf(
"[%d] Type1 ERRORS = ",SourceMatrix.Comm().MyPID());
1488 for(
int i=0; i<Len1; i++)
1489 printf(
"(%2d,%2d,%2d) ",ExportLID1[i],
DomainMap_.GID(ExportLID1[i]),ExportPID1[i]);
1492 throw std::runtime_error(
"Importer1 fails the sanity test");
1495 for(i=1; i<Len2; i++) {
1496 if(ExportPID2[i] < ExportPID2[i-1] || (ExportPID2[i] == ExportPID2[i-1] &&
DomainMap_.GID(ExportLID2[i]) <
DomainMap_.GID(ExportLID2[i-1])))
1500 SourceMatrix.Comm().Barrier();
1502 printf(
"[%d] Type2 ERRORS = ",SourceMatrix.Comm().MyPID());
1503 for(
int i=0; i<Len2; i++)
1504 printf(
"(%2d,%2d,%2d) ",ExportLID2[i],
DomainMap_.GID(ExportLID2[i]),ExportPID2[i]);
1507 throw std::runtime_error(
"Importer2 fails the sanity test");
1510 for(i=1; i<Len3; i++) {
1511 if(ExportPID3[i] < ExportPID3[i-1] || (ExportPID3[i] == ExportPID3[i-1] &&
DomainMap_.GID(ExportLID3[i]) <
DomainMap_.GID(ExportLID3[i-1])))
1515 SourceMatrix.Comm().Barrier();
1517 printf(
"[%d] Type3 ERRORS = ",SourceMatrix.Comm().MyPID());
1518 for(
int i=0; i<Len3; i++)
1519 printf(
"(%2d,%2d,%2d) ",ExportLID3[i],
DomainMap_.GID(ExportLID3[i]),ExportPID3[i]);
1522 throw std::runtime_error(
"Importer3 fails the sanity test");
1529 if(Importer1 && !Importer1->SourceMap().SameAs(
DomainMap_))
1530 throw std::runtime_error(
"ERROR: Map Mismatch Importer1");
1532 if(!Importer2.SourceMap().SameAs(SourceMatrix.RowMap()))
1533 throw std::runtime_error(
"ERROR: Map Mismatch Importer2");
1535 int_type InfGID = std::numeric_limits<int_type>::max();
1536 int InfPID = INT_MAX;
1538 int i1=0, i2=0, i3=0, current=0;
1540 int MyLen=Len1+Len2+Len3;
1541 ExportLIDs.resize(MyLen);
1542 ExportPIDs.resize(MyLen);
1544 while(i1 < Len1 || i2 < Len2 || i3 < Len3){
1545 int PID1 = (i1<Len1)?(ExportPID1[i1]):InfPID;
1546 int PID2 = (i2<Len2)?(ExportPID2[i2]):InfPID;
1547 int PID3 = (i3<Len3)?(ExportPID3[i3]):InfPID;
1549 int_type GID1 = (i1<Len1)?((int_type)
DomainMap_.GID64(ExportLID1[i1])):InfGID;
1550 int_type GID2 = (i2<Len2)?((int_type)
DomainMap_.GID64(ExportLID2[i2])):InfGID;
1551 int_type GID3 = (i3<Len3)?((int_type)
DomainMap_.GID64(ExportLID3[i3])):InfGID;
1553 int MIN_PID =
MIN3(PID1,PID2,PID3);
1554 int_type MIN_GID =
MIN3( ((PID1==MIN_PID)?GID1:InfGID), ((PID2==MIN_PID)?GID2:InfGID), ((PID3==MIN_PID)?GID3:InfGID));
1555 bool added_entry=
false;
1558 if(PID1 == MIN_PID && GID1 == MIN_GID){
1559 ExportLIDs[current] = ExportLID1[i1];
1560 ExportPIDs[current] = ExportPID1[i1];
1567 if(PID2 == MIN_PID && GID2 == MIN_GID){
1569 ExportLIDs[current] = ExportLID2[i2];
1570 ExportPIDs[current] = ExportPID2[i2];
1578 if(PID3 == MIN_PID && GID3 == MIN_GID){
1580 ExportLIDs[current] = ExportLID3[i3];
1581 ExportPIDs[current] = ExportPID3[i3];
1587 if(current!=MyLen) {
1588 ExportLIDs.resize(current);
1589 ExportPIDs.resize(current);
1596 template<
typename ImportType,
typename int_type>
1601 #ifdef ENABLE_MMM_TIMINGS 1603 if(label) tpref = std::string(label);
1604 using Teuchos::TimeMonitor;
1605 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1"))));
1614 int MyPID = SourceMatrix.Comm().MyPID();
1617 std::vector<int> ReverseSendSizes;
1618 std::vector<int_type> ReverseSendBuffer;
1619 std::vector<int> ReverseRecvSizes;
1620 int_type * ReverseRecvBuffer=0;
1623 bool communication_needed = RowImporter.SourceMap().DistributedGlobal();
1632 if (!SourceMatrix.RowMap().SameAs(RowImporter.SourceMap()))
1633 throw "LightweightCrsMatrix: Fused copy constructor requires Importer.SourceMap() to match SourceMatrix.RowMap()";
1636 int NumSameIDs = RowImporter.NumSameIDs();
1637 int NumPermuteIDs = RowImporter.NumPermuteIDs();
1638 int NumRemoteIDs = RowImporter.NumRemoteIDs();
1639 int NumExportIDs = RowImporter.NumExportIDs();
1640 int* ExportLIDs = RowImporter.ExportLIDs();
1641 int* RemoteLIDs = RowImporter.RemoteLIDs();
1642 int* PermuteToLIDs = RowImporter.PermuteToLIDs();
1643 int* PermuteFromLIDs = RowImporter.PermuteFromLIDs();
1646 Epetra_Distributor& Distor = RowImporter.Distributor();
1647 const Epetra_MpiComm * MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&SourceMatrix.Comm());
1648 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1656 std::vector<int> SourcePids;
1657 std::vector<int> TargetPids;
1666 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in CheckSizes()";
1672 int LenExports_ = 0;
1676 bool VarSizes =
false;
1677 if( NumExportIDs > 0) {
1678 Sizes_ =
new int[NumExportIDs];
1683 const Epetra_Import *MyImporter= SourceMatrix.Importer();
1684 if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,
false);
1686 SourcePids.resize(SourceMatrix.ColMap().NumMyElements());
1687 SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),SourceMatrix.Comm().MyPID());
1690 #ifdef ENABLE_MMM_TIMINGS 1691 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1.1 Forward Pack"))));
1695 rv=Epetra_Import_Util::PackAndPrepareWithOwningPIDs(SourceMatrix,
1696 NumExportIDs,ExportLIDs,
1697 LenExports_,Exports_,SizeOfPacket,
1698 Sizes_,VarSizes,SourcePids);
1699 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in PackAndPrepare()";
1702 #ifdef ENABLE_MMM_TIMINGS 1703 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1.2 Reverse"))));
1706 if (communication_needed) {
1709 const int * ExportPIDs = RowImporter.ExportPIDs();
1713 std::vector<int> SendSizes(MDistor->NumSends()+1,0);
1714 for(
int i=0, curr_pid=0; i<NumExportIDs; i++) {
1715 if(i>0 && ExportPIDs[i] > ExportPIDs[i-1]) curr_pid++;
1716 SendSizes[curr_pid] +=Sizes_[i];
1719 if(i>0 && ExportPIDs[i] < ExportPIDs[i-1])
throw "ExportPIDs not sorted";
1722 #ifdef ENABLE_MMM_TIMINGS 1723 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1.2 Forward Send"))));
1726 std::vector<int> RecvSizes(MDistor->NumReceives()+1);
1727 int msg_tag=MpiComm->GetMpiTag();
1728 boundary_exchange_varsize<char>(*MpiComm,MPI_CHAR,MDistor->NumSends(),MDistor->ProcsTo(),SendSizes.size() ? &SendSizes[0] : 0,Exports_,
1729 MDistor->NumReceives(),MDistor->ProcsFrom(),RecvSizes.size() ? &RecvSizes[0] : 0,Imports_,SizeOfPacket,msg_tag);
1732 if(SourceMatrix.Importer()) {
1733 Epetra_Import* SourceImporter=
const_cast<Epetra_Import*
>(SourceMatrix.Importer());
1734 const Epetra_MpiDistributor * MyDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&SourceImporter->Distributor());
1737 #ifdef ENABLE_MMM_TIMINGS 1738 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1.3 Reverse Pack"))));
1743 PackAndPrepareReverseComm<ImportType, int_type>(SourceMatrix,RowImporter,ReverseSendSizes,ReverseSendBuffer);
1745 #ifdef ENABLE_MMM_TIMINGS 1746 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-1.4 Reverse Send"))));
1751 ReverseRecvSizes.resize(MyDistor->NumSends()+1);
1752 const int msg_tag2 = MpiComm->GetMpiTag ();
1753 MPI_Datatype data_type =
sizeof(int_type) == 4 ? MPI_INT : MPI_LONG_LONG;
1754 boundary_exchange_varsize<int_type> (*MpiComm, data_type, MyDistor->NumReceives (),
1755 MyDistor->ProcsFrom (),
1756 ReverseSendSizes.size() ? &ReverseSendSizes[0] : 0,
1757 ReverseSendBuffer.size() ? &ReverseSendBuffer[0] : 0,
1758 MyDistor->NumSends (), MyDistor->ProcsTo (),
1759 ReverseRecvSizes.size() ? &ReverseRecvSizes[0] : 0,
1760 ReverseRecvBuffer, 1, msg_tag2);
1765 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in Distor.Do";
1770 #ifdef ENABLE_MMM_TIMINGS 1771 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-2"))));
1775 int mynnz = Epetra_Import_Util::UnpackWithOwningPIDsCount(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_);
1779 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1780 if(
sizeof(int_type) ==
sizeof(
long long))
1783 vals_.resize(mynnz);
1786 if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,
true);
1788 SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),-1);
1793 double * myvals =
vals_.size() ? &
vals_[0] : 0;
1794 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1795 if(
sizeof(int_type) ==
sizeof(
int)) {
1797 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,
N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1801 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1802 if(
sizeof(int_type) ==
sizeof(
long long)) {
1804 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,
N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1808 throw "EpetraExt::LightweightCrsMatrix::Construct: sizeof(int_type) error.";
1813 #ifdef ENABLE_MMM_TIMINGS 1814 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-3"))));
1818 MakeColMapAndReindex<int_type>(TargetPids,getcolind<int_type>(),SortGhosts);
1824 MakeExportLists<ImportType, int_type>(SourceMatrix,RowImporter,ReverseRecvSizes,ReverseRecvBuffer,
ExportPIDs_,
ExportLIDs_);
1831 #ifdef ENABLE_MMM_TIMINGS 1832 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-4"))));
1839 #ifdef ENABLE_MMM_TIMINGS 1840 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS C-5"))));
1844 delete [] ReverseRecvBuffer;
1861 DomainMap_(SourceMatrix.DomainMap())
1863 #ifdef ENABLE_MMM_TIMINGS 1865 if(label) tpref = std::string(label);
1866 using Teuchos::TimeMonitor;
1867 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string(
"EpetraExt: LWCRS Total"))));
1872 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1873 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1874 Construct<RemoteOnlyImport, int>(SourceMatrix,RowImporter,SortGhosts,label);
1878 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1879 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1880 Construct<RemoteOnlyImport, long long>(SourceMatrix,RowImporter,SortGhosts,label);
1884 throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1894 DomainMap_(SourceMatrix.DomainMap())
1896 RowMapEP_=
new Epetra_BlockMap(RowImporter.TargetMap());
1897 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1898 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1899 Construct<Epetra_Import, int>(SourceMatrix,RowImporter);
1903 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1904 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1905 Construct<Epetra_Import, long long>(SourceMatrix,RowImporter);
1909 throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1922 template <
class TransferType>
1924 if(!Transfer)
return;
1926 const Epetra_MpiDistributor & Distor =
dynamic_cast<Epetra_MpiDistributor&
>(Transfer->Distributor());
1927 const Epetra_Comm & Comm = Transfer->SourceMap().Comm();
1929 int rows_send = Transfer->NumExportIDs();
1930 int rows_recv = Transfer->NumRemoteIDs();
1932 int round1_send = Transfer->NumExportIDs() *
sizeof(int);
1933 int round1_recv = Transfer->NumRemoteIDs() *
sizeof(int);
1934 int num_send_neighbors = Distor.NumSends();
1935 int num_recv_neighbors = Distor.NumReceives();
1936 int round2_send, round2_recv;
1937 Distor.GetLastDoStatistics(round2_send,round2_recv);
1939 int myPID = Comm.MyPID();
1940 int NumProcs = Comm.NumProc();
1947 int lstats[8] = {num_send_neighbors,num_recv_neighbors,rows_send,rows_recv,round1_send,round1_recv,round2_send,round2_recv};
1948 int gstats_min[8], gstats_max[8];
1950 double lstats_avg[8], gstats_avg[8];
1951 for(
int i=0; i<8; i++)
1952 lstats_avg[i] = ((
double)lstats[i])/NumProcs;
1954 Comm.MinAll(lstats,gstats_min,8);
1955 Comm.MaxAll(lstats,gstats_max,8);
1956 Comm.SumAll(lstats_avg,gstats_avg,8);
1959 printf(
"%s Send Statistics[min/avg/max]: neigh=%d/%4.1f/%d rows=%d/%4.1f/%d round1=%d/%4.1f/%d round2=%d/%4.1f/%d\n", label.c_str(),
1960 (int)gstats_min[0],gstats_avg[0],(
int)gstats_max[0], (int)gstats_min[2],gstats_avg[2],(
int)gstats_max[2],
1961 (int)gstats_min[4],gstats_avg[4],(
int)gstats_max[4], (int)gstats_min[6],gstats_avg[6],(
int)gstats_max[6]);
1962 printf(
"%s Recv Statistics[min/avg/max]: neigh=%d/%4.1f/%d rows=%d/%4.1f/%d round1=%d/%4.1f/%d round2=%d/%4.1f/%d\n", label.c_str(),
1963 (int)gstats_min[1],gstats_avg[1],(
int)gstats_max[1], (int)gstats_min[3],gstats_avg[3],(
int)gstats_max[3],
1964 (int)gstats_min[5],gstats_avg[5],(
int)gstats_max[5], (int)gstats_min[7],gstats_avg[7],(
int)gstats_max[7]);
1984 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1988 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
void Construct_int(int NumGlobalElements, int NumMyElements, const int *MyGlobalElements, long long IndexBase, bool GenerateHash=true)
LightweightCrsMatrix * importMatrix
const LightweightMap * TargetMap_
RemoteOnlyImport(const Epetra_Import &Importer, LightweightMap &RemoteOnlyTargetMap)
Epetra_HashTable< long long > * LIDHash_LL_
Epetra_BlockMap * RowMapEP_
std::vector< int > targetMapToOrigRow
virtual ~CrsWrapper_GraphBuilder()
void debug_print_distor(const char *label, const Epetra_Distributor *Distor, const Epetra_Comm &Comm)
const LightweightMap & TargetMap() const
std::pair< int_type, int_type > Tget_col_range(const Epetra_CrsMatrix &mtx)
int MakeColMapAndReindex(std::vector< int > owningPIDs, std::vector< GO > Gcolind, bool SortGhosts=false, const char *label=0)
LightweightMap * RowMapLW_
void pack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int > &proc_col_ranges, std::vector< int > &send_rows, std::vector< int > &rows_per_send_proc)
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
int PackAndPrepareReverseComm(const Epetra_CrsMatrix &SourceMatrix, ImportType &RowImporter, std::vector< int > &ReverseSendSizes, std::vector< int_type > &ReverseSendBuffer)
Epetra_CrsMatrix & ecrsmat_
void debug_compare_import(const Epetra_Import *Import1, const Epetra_Import *Import2)
const Epetra_Map * rowMap
std::vector< int > colind_
void printMultiplicationStatistics(Epetra_Import *Transfer, const std::string &label)
const Epetra_BlockMap * SourceMap_
std::vector< int > ExportPIDs_
const Epetra_CrsMatrix * origMatrix
void MakeColMapAndReindexSort(int &NumRemoteColGIDs, GO *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs, bool SortGhostsAssociatedWithEachProcessor_)
long long GID64(int LID) const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
void build_type1_exports_sort(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< int_type > &ExportGID1, int total_length1)
void build_type1_exports_sort< int >(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< int > &ExportGID1, int total_length1)
int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
std::pair< long long, long long > get_col_range< long long >(const Epetra_Map &emap)
std::vector< int > targetMapToImportRow
long long * MyGlobalElements64() const
LightweightMap & operator=(const LightweightMap &map)
LightweightMapData * Data_
void build_type3_exports_sort< int >(std::vector< int > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
void build_type1_exports_sort< long long >(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< long long > &ExportGID1, int total_length1)
virtual ~CrsWrapper_Epetra_CrsMatrix()
void build_type3_exports_sort(std::vector< int_type > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
int * MyGlobalElements() const
int build_type2_exports(const Epetra_CrsMatrix &SourceMatrix, ImportType &MyImporter, std::vector< int > &ExportLID2, std::vector< int > &ExportPID2)
void TPrintMultiplicationStatistics(TransferType *Transfer, const std::string &label)
int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
bool GlobalIndicesInt() const
std::vector< double > vals_
CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix &epetracrsmatrix)
std::vector< int > MyGlobalElements_int_
std::map< int_type, std::set< int_type > * > graph_
void MakeColMapAndReindexSort< long long >(int &NumRemoteColGIDs, long long *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs, bool SortGhostsAssociatedWithEachProcessor_)
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)
void MakeColMapAndReindexSort< int >(int &NumRemoteColGIDs, int *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs, bool SortGhostsAssociatedWithEachProcessor_)
bool GlobalIndicesLongLong() const
void Construct_LL(long long NumGlobalElements, int NumMyElements, const long long *MyGlobalElements, long long IndexBase, bool GenerateHash=true)
void Tpack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int_type > &proc_col_ranges, std::vector< int_type > &send_rows, std::vector< int > &rows_per_send_proc)
const Epetra_Map & RowMap() const
CrsWrapper_GraphBuilder(const Epetra_Map &emap)
virtual ~CrsMatrixStruct()
Epetra_HashTable< int > * LIDHash_int_
std::vector< long long > MyGlobalElements_LL_
int build_type1_exports(const Epetra_Import *Importer1, std::vector< int > &ExportLID1, std::vector< int > &ExportPID1)
LightweightCrsMatrix(const Epetra_CrsMatrix &A, RemoteOnlyImport &RowImporter, bool SortGhosts=false, const char *label=0)
std::map< int_type, std::set< int_type > * > & get_graph()
void build_type3_exports_sort< long long >(std::vector< long long > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
std::vector< int > ColMapOwningPIDs_
std::vector< int > ExportLIDs_
std::vector< int > rowptr_
const Epetra_Map * colMap
const Epetra_BlockMap * importColMap
void insert_matrix_locations(CrsWrapper_GraphBuilder< int_type > &graphbuilder, Epetra_CrsMatrix &C)
int dumpCrsMatrixStruct(const CrsMatrixStruct &M)
int NumMyElements() const
int InsertGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
void Construct(const Epetra_CrsMatrix &A, ImportType &RowImporter, bool SortGhosts=false, const char *label=0)
std::vector< long long > colind_LL_
std::pair< int, int > get_col_range< int >(const Epetra_Map &emap)
int MakeExportLists(const Epetra_CrsMatrix &SourceMatrix, ImportType &RowImporter, std::vector< int > &ReverseRecvSizes, const int_type *ReverseRecvBuffer, std::vector< int > &ExportPIDs, std::vector< int > &ExportLIDs)
Epetra_Distributor * Distor_