44 #include <Epetra_Comm.h> 45 #include <Epetra_CrsMatrix.h> 46 #include <Epetra_Import.h> 47 #include <Epetra_Distributor.h> 48 #include <Epetra_HashTable.h> 49 #include <Epetra_Util.h> 50 #include <Epetra_Import_Util.h> 51 #include <Epetra_GIDTypeSerialDenseVector.h> 53 #include <Teuchos_TimeMonitor.hpp> 57 #include "Epetra_MpiComm.h" 58 #include "Epetra_MpiDistributor.h" 60 #define MIN(x,y) ((x)<(y)?(x):(y)) 61 #define MIN3(x,y,z) ((x)<(y)?(MIN(x,z)):(MIN(y,z))) 67 void debug_print_distor(
const char * label,
const Epetra_Distributor * Distor,
const Epetra_Comm & Comm) {
69 const Epetra_MpiDistributor * MDistor =
dynamic_cast<const Epetra_MpiDistributor*
>(Distor);
70 printf(
"[%d] %s\n",Comm.MyPID(),label);
71 printf(
"[%d] NumSends = %d NumRecvs = %d\n",Comm.MyPID(),MDistor->NumSends(),MDistor->NumReceives());
72 printf(
"[%d] ProcsTo = ",Comm.MyPID());
73 for(
int ii=0; ii<MDistor->NumSends(); ii++)
74 printf(
"%d ",MDistor->ProcsTo()[ii]);
75 printf(
"\n[%d] ProcsFrom = ",Comm.MyPID());
76 for(
int ii=0; ii<MDistor->NumReceives(); ii++)
77 printf(
"%d ",MDistor->ProcsFrom()[ii]);
86 if(!Import1 && !Import2)
return;
87 const Epetra_Comm & Comm = (Import1)? Import1->SourceMap().Comm() : Import2->SourceMap().Comm();
90 if( (!Import1 && Import2) || (Import2 && !Import1) ) {printf(
"[%d] DCI: One Import exists, the other does not\n",PID);
return;}
91 if(!Import1->SourceMap().SameAs(Import2->SourceMap())) {printf(
"[%d] DCI: SourceMaps don't match\n",PID);
return;}
92 if(!Import1->TargetMap().SameAs(Import2->TargetMap())) {printf(
"[%d] DCI: TargetMaps don't match\n",PID);
return;}
94 if(Import1->NumSameIDs() != Import2->NumSameIDs()) {printf(
"[%d] DCI NumSameIDs() mismatch %d vs. %d\n",PID,Import1->NumSameIDs(),Import2->NumSameIDs());flag=
false;}
96 if(Import1->NumPermuteIDs() != Import2->NumPermuteIDs()) {printf(
"[%d] DCI NumPermuteIDs() mismatch %d vs. %d\n",PID,Import1->NumPermuteIDs(),Import2->NumPermuteIDs()); flag=
false;}
98 if(Import1->NumExportIDs() != Import2->NumExportIDs()) {printf(
"[%d] DCI NumExportIDs() mismatch %d vs. %d\n",PID,Import1->NumExportIDs(),Import2->NumExportIDs()); flag=
false;}
100 if(Import1->NumRemoteIDs() != Import2->NumRemoteIDs()) {printf(
"[%d] DCI NumRemoteIDs() mismatch %d vs. %d\n",PID,Import1->NumRemoteIDs(),Import2->NumRemoteIDs()); flag=
false;}
102 if(Import1->NumSend() != Import2->NumSend()) {printf(
"[%d] DCI NumSend() mismatch %d vs. %d\n",PID,Import1->NumSend(),Import2->NumSend()); flag=
false;}
104 if(Import1->NumRecv() != Import2->NumRecv()) {printf(
"[%d] DCI NumRecv() mismatch %d vs. %d\n",PID,Import1->NumRecv(),Import2->NumRecv()); flag=
false;}
107 if(flag) printf(
"[%d] DCI Importers compare OK\n",PID);
109 Import1->SourceMap().Comm().Barrier();
110 Import1->SourceMap().Comm().Barrier();
111 Import1->SourceMap().Comm().Barrier();
119 : numRows(0), numEntriesPerRow(NULL), indices(NULL), values(NULL),
120 remote(NULL), numRemote(0), importColMap(NULL), rowMap(NULL), colMap(NULL),
121 domainMap(NULL), importMatrix(NULL), origMatrix(NULL)
148 std::cout <<
"proc " << M.
rowMap->Comm().MyPID()<<std::endl;
149 std::cout <<
"numRows: " << M.
numRows<<std::endl;
150 for(
int i=0; i<M.
numRows; ++i) {
153 std::cout <<
" *"<<M.
rowMap->GID64(i)<<
" " 157 std::cout <<
" "<<M.
rowMap->GID64(i)<<
" " 166 : ecrsmat_(epetracrsmatrix)
177 return ecrsmat_.RowMap();
182 return ecrsmat_.Filled();
185 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 189 return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
195 return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
199 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 203 return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
209 return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
216 template<
typename int_type>
222 int num_rows = emap.NumMyElements();
224 emap.MyGlobalElementsPtr(rows);
226 for(
int i=0; i<num_rows; ++i) {
227 graph_[rows[i]] =
new std::set<int_type>;
231 template<
typename int_type>
234 typename std::map<int_type,std::set<int_type>*>::iterator
235 iter = graph_.begin(), iter_end = graph_.end();
236 for(; iter!=iter_end; ++iter) {
243 template<
typename int_type>
249 template<
typename int_type>
253 typename std::map<int_type,std::set<int_type>*>::iterator
254 iter = graph_.find(GlobalRow);
256 if (iter == graph_.end())
return(-1);
258 std::set<int_type>& cols = *(iter->second);
260 for(
int i=0; i<NumEntries; ++i) {
261 cols.insert(Indices[i]);
264 int row_length = cols.size();
265 if (row_length > max_row_length_) max_row_length_ = row_length;
270 template<
typename int_type>
277 template<
typename int_type>
278 std::map<int_type,std::set<int_type>*>&
284 template<
typename int_type>
289 if (max_row_length < 1)
return;
291 std::vector<int_type> indices(max_row_length);
292 int_type* indices_ptr = &indices[0];
293 std::vector<double> zeros(max_row_length, 0.0);
294 double* zeros_ptr = &zeros[0];
296 std::map<int_type,std::set<int_type>*>& graph = graphbuilder.
get_graph();
298 typename std::map<int_type,std::set<int_type>*>::iterator
299 iter = graph.begin(), iter_end = graph.end();
301 for(; iter!=iter_end; ++iter) {
302 int_type row = iter->first;
303 std::set<int_type>& cols = *(iter->second);
304 int num_entries = cols.size();
306 typename std::set<int_type>::iterator
307 col_iter = cols.begin(), col_end = cols.end();
308 for(
int j=0; col_iter!=col_end; ++col_iter, ++j) {
309 indices_ptr[j] = *col_iter;
312 C.InsertGlobalValues(row, num_entries, zeros_ptr, indices_ptr);
316 template<
typename int_type>
318 const std::vector<int_type>& proc_col_ranges,
319 std::vector<int_type>& send_rows,
320 std::vector<int>& rows_per_send_proc)
322 const Epetra_Map& rowmap = mtx.RowMap();
323 int numrows = mtx.NumMyRows();
324 const Epetra_CrsGraph& graph = mtx.Graph();
326 int* col_indices = NULL;
327 int_type* Tcol_indices = NULL;
328 int num_col_ranges = proc_col_ranges.size()/2;
329 rows_per_send_proc.resize(num_col_ranges);
331 for(
int nc=0; nc<num_col_ranges; ++nc) {
332 int_type first_col = proc_col_ranges[nc*2];
333 int_type last_col = proc_col_ranges[nc*2+1];
334 int num_send_rows = 0;
335 for(
int i=0; i<numrows; ++i) {
336 int_type grow = (int_type) rowmap.GID64(i);
338 const Epetra_Map& colmap = mtx.ColMap();
339 graph.ExtractMyRowView(i, rowlen, col_indices);
340 for(
int j=0; j<rowlen; ++j) {
341 int_type col = (int_type) colmap.GID64(col_indices[j]);
342 if (first_col <= col && last_col >= col) {
344 send_rows.push_back(grow);
350 graph.ExtractGlobalRowView(grow, rowlen, Tcol_indices);
351 for(
int j=0; j<rowlen; ++j) {
352 if (first_col <= Tcol_indices[j] && last_col >= Tcol_indices[j]) {
354 send_rows.push_back(grow);
360 rows_per_send_proc[nc] = num_send_rows;
364 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 366 const std::vector<int>& proc_col_ranges,
367 std::vector<int>& send_rows,
368 std::vector<int>& rows_per_send_proc)
370 if(mtx.RowMatrixRowMap().GlobalIndicesInt()) {
371 Tpack_outgoing_rows<int>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
374 throw "EpetraExt::pack_outgoing_rows: Global indices not int";
379 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 381 const std::vector<long long>& proc_col_ranges,
382 std::vector<long long>& send_rows,
383 std::vector<int>& rows_per_send_proc)
385 if(mtx.RowMatrixRowMap().GlobalIndicesLongLong()) {
386 Tpack_outgoing_rows<long long>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
389 throw "EpetraExt::pack_outgoing_rows: Global indices not long long";
394 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 398 if(emap.GlobalIndicesInt())
399 return std::make_pair(emap.MinMyGID(),emap.MaxMyGID());
401 throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_Map";
405 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 409 if(emap.GlobalIndicesLongLong())
410 return std::make_pair(emap.MinMyGID64(),emap.MaxMyGID64());
412 throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_Map";
416 template<
typename int_type>
419 std::pair<int_type,int_type> col_range;
421 col_range = get_col_range<int_type>(mtx.ColMap());
424 const Epetra_Map& row_map = mtx.RowMap();
425 col_range.first = row_map.MaxMyGID64();
426 col_range.second = row_map.MinMyGID64();
428 int_type* col_indices = NULL;
429 const Epetra_CrsGraph& graph = mtx.Graph();
430 for(
int i=0; i<row_map.NumMyElements(); ++i) {
431 graph.ExtractGlobalRowView((int_type) row_map.GID64(i), rowlen, col_indices);
432 for(
int j=0; j<rowlen; ++j) {
433 if (col_indices[j] < col_range.first) col_range.first = col_indices[j];
434 if (col_indices[j] > col_range.second) col_range.second = col_indices[j];
442 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 446 if(mtx.RowMatrixRowMap().GlobalIndicesInt())
447 return Tget_col_range<int>(mtx);
449 throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_CrsMatrix";
453 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 457 if(mtx.RowMatrixRowMap().GlobalIndicesLongLong())
458 return Tget_col_range<long long>(mtx);
460 throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_CrsMatrix";
469 template <
typename MyType>
470 void boundary_exchange(
const Epetra_MpiComm Comm, MPI_Datatype DataType,
471 int NumSends,
const int * SendProcs,
const int * SendSizes, MyType* SendBuffer,
472 int NumRecvs,
const int * RecvProcs,
const int * RecvSizes, MyType* RecvBuffer,
int SizeOfPacket,
int msg_tag)
475 MPI_Comm comm = Comm.Comm();
476 std::vector<MPI_Request> requests(NumRecvs);
477 std::vector<MPI_Status> status(NumRecvs);
479 int i,num_waits=0,MyPID=Comm.MyPID();
480 int start, self_recv_len=-1,self_recv_start=-1, self_send_start=-1;
483 int mysendsize=1, myrecvsize=1;
487 for(i=0; i<NumRecvs; i++){
488 if(RecvSizes) myrecvsize=RecvSizes[i]*SizeOfPacket;
489 if(RecvProcs[i] != MyPID) {
490 MPI_Irecv(RecvBuffer + start, myrecvsize, DataType, RecvProcs[i], msg_tag, comm, &requests[num_waits]);
494 self_recv_len = myrecvsize;
495 self_recv_start=start;
502 for(i=0; i<NumSends; i++){
503 if(SendSizes) mysendsize=SendSizes[i]*SizeOfPacket;
504 if(SendProcs[i] != MyPID)
505 MPI_Send(SendBuffer + start, mysendsize,DataType,SendProcs[i],msg_tag,comm);
507 self_send_start=start;
512 if(self_recv_len != -1)
513 memcpy(RecvBuffer+self_recv_start,SendBuffer+self_send_start,self_recv_len*
sizeof(MyType)*SizeOfPacket);
517 MPI_Waitall(num_waits, &requests[0],&status[0]);
523 template <
typename MyType>
524 void boundary_exchange_varsize(
const Epetra_MpiComm Comm, MPI_Datatype DataType,
525 int NumSends,
const int * SendProcs,
const int * SendSizes, MyType* SendBuffer,
526 int NumRecvs,
const int * RecvProcs,
int * RecvSizes, MyType*& RecvBuffer,
int SizeOfPacket,
int msg_tag)
532 boundary_exchange<int>(Comm,MPI_INT,NumSends,SendProcs,(
int*)0,const_cast<int*>(SendSizes),NumRecvs,RecvProcs,(
int*)0,RecvSizes,1,msg_tag);
535 for(i=0; i<NumRecvs; i++) rbuffersize+=RecvSizes[i]*SizeOfPacket;
536 RecvBuffer =
new MyType[rbuffersize];
539 boundary_exchange<MyType>(Comm,DataType,NumSends,SendProcs,SendSizes,SendBuffer,NumRecvs,RecvProcs,RecvSizes,RecvBuffer,SizeOfPacket,msg_tag+100);
550 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
553 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
561 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 564 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 574 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 575 void LightweightMap::Construct_int(
int ,
int numMyElements,
const int * myGlobalElements,
long long ,
bool GenerateHash)
581 if(GenerateHash) Data_->
LIDHash_int_ =
new Epetra_HashTable<int>(numMyElements + 1 );
582 for(
int i=0; i < numMyElements; ++i ) {
584 if(GenerateHash) Data_->
LIDHash_int_->Add(myGlobalElements[i], i);
591 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 592 void LightweightMap::Construct_LL(
long long ,
int numMyElements,
const long long * myGlobalElements,
long long ,
bool GenerateHash)
598 if(GenerateHash) Data_->
LIDHash_LL_ =
new Epetra_HashTable<long long>(numMyElements + 1 );
599 for(
int i=0; i < numMyElements; ++i ) {
601 if(GenerateHash) Data_->
LIDHash_LL_->Add(myGlobalElements[i], i);
608 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 611 Construct_int(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
615 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 618 Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
623 Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
631 Data_->
CopyMap_=
new Epetra_Map(Map);
632 IsLongLong = Map.GlobalIndicesLongLong();
633 IsInt = Map.GlobalIndicesInt();
640 Data_->IncrementReferenceCount();
641 IsLongLong = map.IsLongLong;
653 if((
this != &map) && (Data_ != map.Data_)) {
656 Data_->IncrementReferenceCount();
658 IsLongLong = map.IsLongLong;
664 void LightweightMap::CleanupData(){
666 Data_->DecrementReferenceCount();
667 if(Data_->ReferenceCount() == 0) {
673 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 678 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 687 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 692 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 697 throw "EpetraExt::LightweightMap::NumMyElements: Global indices unknowns";
701 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 707 throw "EpetraExt::LightweightMap::LID: Int version called for long long map";
709 throw "EpetraExt::LightweightMap::LID: unknown GID type";
712 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 719 throw "EpetraExt::LightweightMap::LID: Long long version called for int map";
721 throw "EpetraExt::LightweightMap::LID: unknown GID type";
726 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 735 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 742 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 749 throw "EpetraExt::LightweightMap::GID64: Global indices unknown.";
753 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 760 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 778 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 783 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 788 throw "EpetraExt::LightweightMap::MaxLID: Global indices unknowns";
801 SourceMap_=&Importer.SourceMap();
802 TargetMap_=&RemoteOnlyTargetMap;
805 NumSend_ = Importer.NumSend();
806 NumRemoteIDs_ = Importer.NumRemoteIDs();
807 NumExportIDs_ = Importer.NumExportIDs();
808 Distor_ = &Importer.Distributor();
809 int * OldRemoteLIDs = Importer.RemoteLIDs();
810 int * OldExportLIDs = Importer.ExportLIDs();
811 int * OldExportPIDs = Importer.ExportPIDs();
815 throw std::runtime_error(
"RemoteOnlyImport: Importer doesn't match RemoteOnlyTargetMap for number of remotes.");
818 ExportLIDs_ =
new int[NumExportIDs_];
819 ExportPIDs_ =
new int[NumExportIDs_];
820 for(i=0; i<NumExportIDs_; i++) {
821 ExportLIDs_[i] = OldExportLIDs[i];
822 ExportPIDs_[i] = OldExportPIDs[i];
828 RemoteLIDs_ =
new int[NumRemoteIDs_];
829 if(TargetMap_->GlobalIndicesInt()) {
830 for(i=0; i<NumRemoteIDs_; i++)
831 RemoteLIDs_[i] = TargetMap_->
LID( (
int) Importer.TargetMap().GID64(OldRemoteLIDs[i]));
833 else if(TargetMap_->GlobalIndicesLongLong()) {
834 for(i=0; i<NumRemoteIDs_; i++)
835 RemoteLIDs_[i] = TargetMap_->LID(Importer.TargetMap().GID64(OldRemoteLIDs[i]));
838 throw std::runtime_error(
"RemoteOnlyImport: TargetMap_ index type unknown.");
841 for(i=0; i<NumRemoteIDs_-1; i++)
842 if(RemoteLIDs_[i] > RemoteLIDs_[i+1])
843 throw std::runtime_error(
"RemoteOnlyImport: Importer and RemoteOnlyTargetMap order don't match.");
849 delete [] ExportLIDs_;
850 delete [] ExportPIDs_;
851 delete [] RemoteLIDs_;
861 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs);
865 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
870 if(NumRemoteColGIDs > 0) {
871 SortLists[0] = RemoteColindices;
872 SortLists[1] = &RemotePermuteIDs[0];
873 Epetra_Util::Sort(
true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, NLists, SortLists);
877 bool SortGhostsAssociatedWithEachProcessor_ =
false;
878 if (SortGhostsAssociatedWithEachProcessor_) {
884 int StartCurrent, StartNext;
885 StartCurrent = 0; StartNext = 1;
886 while ( StartNext < NumRemoteColGIDs ) {
887 if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
889 SortLists[0] = &RemotePermuteIDs[StartCurrent];
890 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists);
891 StartCurrent = StartNext; StartNext++;
894 SortLists[0] = &RemotePermuteIDs[StartCurrent];
895 Epetra_Util::Sort(
true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists);
901 std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
904 if(NumRemoteColGIDs > 0) {
905 long long* SortLists_LL[1] = {RemoteColindices};
906 int* SortLists_int[1] = {&RemotePermuteIDs[0]};
907 Epetra_Util::Sort(
true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, 1, SortLists_int, 1, SortLists_LL);
910 bool SortGhostsAssociatedWithEachProcessor_ =
false;
911 if (SortGhostsAssociatedWithEachProcessor_) {
917 int StartCurrent, StartNext;
918 StartCurrent = 0; StartNext = 1;
919 while ( StartNext < NumRemoteColGIDs ) {
920 if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
922 int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
923 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
924 StartCurrent = StartNext; StartNext++;
927 int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
928 Epetra_Util::Sort(
true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
933 int LightweightCrsMatrix::MakeColMapAndReindex(std::vector<int> owningPIDs, std::vector<GO> Gcolind)
935 #ifdef ENABLE_MMM_TIMINGS 936 using Teuchos::TimeMonitor;
937 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS-3.1")));
943 int numDomainElements = DomainMap_.NumMyElements();
944 std::vector<bool> LocalGIDs(numDomainElements,
false);
946 bool DoSizes = !DomainMap_.ConstantElementSize();
947 if(DoSizes) EPETRA_CHK_ERR(-1);
953 if(use_lw) numMyBlockRows = RowMapLW_->NumMyElements();
954 else numMyBlockRows = RowMapEP_->NumMyElements();
956 int hashsize = numMyBlockRows;
if (hashsize < 100) hashsize = 100;
957 Epetra_HashTable<GO> RemoteGIDs(hashsize);
958 std::vector<GO> RemoteGIDList; RemoteGIDList.reserve(hashsize);
959 std::vector<int> RemoteOwningPIDs; RemoteOwningPIDs.reserve(hashsize);
964 int NumLocalColGIDs = 0;
965 int NumRemoteColGIDs = 0;
966 for(
int i = 0; i < numMyBlockRows; i++) {
967 for(
int j = rowptr_[i]; j < rowptr_[i+1]; j++) {
970 int LID = DomainMap_.LID(GID);
972 bool alreadyFound = LocalGIDs[
LID];
974 LocalGIDs[
LID] =
true;
980 GO hash_value=RemoteGIDs.Get(GID);
981 if(hash_value == -1) {
982 int PID = owningPIDs[j];
983 if(PID==-1) printf(
"[%d] ERROR: Remote PID should not be -1\n",DomainMap_.Comm().MyPID());
984 colind_[j] = numDomainElements + NumRemoteColGIDs;
985 RemoteGIDs.Add(GID, NumRemoteColGIDs);
986 RemoteGIDList.push_back(GID);
987 RemoteOwningPIDs.push_back(PID);
991 colind_[j] = numDomainElements + hash_value;
997 if (DomainMap_.Comm().NumProc()==1) {
998 if (NumRemoteColGIDs!=0) {
999 throw "Some column IDs are not in domainMap. If matrix is rectangular, you must pass in domainMap to FillComplete";
1002 if (NumLocalColGIDs==numDomainElements) {
1003 ColMap_ = DomainMap_;
1013 int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1014 typename Epetra_GIDTypeSerialDenseVector<GO>::impl Colindices;
1015 if(numMyBlockCols > 0)
1016 Colindices.Size(numMyBlockCols);
1017 GO* RemoteColindices = Colindices.Values() + NumLocalColGIDs;
1019 for(
int i = 0; i < NumRemoteColGIDs; i++)
1020 RemoteColindices[i] = RemoteGIDList[i];
1023 std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
1024 for(
int i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
1026 MakeColMapAndReindexSort<GO>(NumRemoteColGIDs, RemoteColindices, RemotePermuteIDs, RemoteOwningPIDs);
1029 std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
1030 for(
int i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
1033 bool use_local_permute=
false;
1034 std::vector<int> LocalPermuteIDs(numDomainElements);
1042 if(NumLocalColGIDs == DomainMap_.NumMyElements()) {
1043 DomainMap_.MyGlobalElements(Colindices.Values());
1047 DomainMap_.MyGlobalElementsPtr(MyGlobalElements);
1048 int NumLocalAgain = 0;
1049 use_local_permute =
true;
1050 for(
int i = 0; i < numDomainElements; i++) {
1052 LocalPermuteIDs[i] = NumLocalAgain;
1053 Colindices[NumLocalAgain++] = MyGlobalElements[i];
1056 assert(NumLocalAgain==NumLocalColGIDs);
1061 ColMapOwningPIDs_.resize(numMyBlockCols);
1062 ColMapOwningPIDs_.assign(numMyBlockCols,DomainMap_.Comm().MyPID());
1063 for(
int i=0;i<NumRemoteColGIDs;i++)
1064 ColMapOwningPIDs_[NumLocalColGIDs+i] = RemoteOwningPIDs[i];
1066 #ifdef ENABLE_MMM_TIMINGS 1067 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS-3.2")));
1071 LightweightMap temp((GO) -1, numMyBlockCols, Colindices.Values(), (GO) DomainMap_.IndexBase64());
1075 #ifdef ENABLE_MMM_TIMINGS 1076 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS-3.3")));
1080 for(
int i=0; i<numMyBlockRows; i++){
1081 for(
int j=rowptr_[i]; j<rowptr_[i+1]; j++){
1083 if(ID < numDomainElements){
1084 if(use_local_permute) colind_[j] = LocalPermuteIDs[colind_[j]];
1089 colind_[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_[j]-numDomainElements];
1093 assert((
size_t)ColMap_.NumMyElements() == ColMapOwningPIDs_.size());
1103 static inline bool lessthan12(std::pair<int,int> i, std::pair<int,int> j){
1104 return ((i.first<j.first) || (i.first==j.first && i.second <j.second));
1109 template<
typename ImportType,
typename int_type>
1110 int LightweightCrsMatrix::PackAndPrepareReverseComm(
const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter,
1111 std::vector<int> &ReverseSendSizes, std::vector<int_type> &ReverseSendBuffer) {
1115 const Epetra_Import *MyImporter= SourceMatrix.Importer();
1116 if(MyImporter == 0)
return -1;
1117 const Epetra_MpiComm * MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&SourceMatrix.Comm());
1118 int MyPID = MpiComm->MyPID();
1121 int NumExportIDs = RowImporter.NumExportIDs();
1122 int* ExportLIDs = RowImporter.ExportLIDs();
1123 int* ExportPIDs = RowImporter.ExportPIDs();
1126 Epetra_Distributor& Distor = MyImporter->Distributor();
1127 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1128 int NumRecvs = MDistor->NumReceives();
1129 int* RemoteLIDs = MyImporter->RemoteLIDs();
1130 const int * ProcsFrom = MDistor->ProcsFrom();
1131 const int * LengthsFrom = MDistor->LengthsFrom();
1136 std::vector<int> RemotePIDOrder(SourceMatrix.NumMyCols(),-1);
1139 for(i=0,j=0;i<NumRecvs;i++){
1140 for(k=0;k<LengthsFrom[i];k++){
1141 int pid=ProcsFrom[i];
1142 if(pid!=MyPID) RemotePIDOrder[RemoteLIDs[j]]=i;
1148 std::vector<std::set<std::pair<int,int_type> > > ReversePGIDs(NumRecvs);
1149 int *rowptr, *colind;
1151 EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1155 for(i=0; i < NumExportIDs; i++) {
1156 int lid = ExportLIDs[i];
1157 int exp_pid = ExportPIDs[i];
1158 for(j=rowptr[lid]; j<rowptr[lid+1]; j++){
1159 int pid_order = RemotePIDOrder[colind[j]];
1161 int_type gid = (int_type) SourceMatrix.GCID64(colind[j]);
1163 ReversePGIDs[pid_order].insert(std::pair<int,int_type>(exp_pid,gid));
1169 ReverseSendSizes.resize(NumRecvs+1);
1171 for(i=0; i<NumRecvs; i++) {
1172 ReverseSendSizes[i] = 2*ReversePGIDs[i].size();
1173 totalsize += ReverseSendSizes[i];
1177 ReverseSendBuffer.resize(totalsize+1);
1178 for(i=0, j=0; i<NumRecvs; i++) {
1179 for(
typename std::set<std::pair<int,int_type> >::iterator it=ReversePGIDs[i].begin(); it!=ReversePGIDs[i].end(); it++) {
1180 ReverseSendBuffer[j] = it->first;
1181 ReverseSendBuffer[j+1] = it->second;
1192 template<
typename int_type>
1198 int * companion = &ExportGID3[0];
1199 Epetra_Util::Sort(
true,total_length3,&ExportPID3[0],0,0,1,&companion,0,0);
1205 long long * companion = &ExportGID3[0];
1206 Epetra_Util::Sort(
true,total_length3,&ExportPID3[0],0,0,0,0,1,&companion);
1209 template<
typename int_type>
1210 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){
1214 int total_length3=0;
1215 for(i=0; i<Nrecv; i++)
1216 total_length3+=ReverseRecvSizes[i]/2;
1217 if(total_length3==0)
return 0;
1219 std::vector<int_type> ExportGID3(total_length3);
1220 ExportLID3.resize(total_length3);
1221 ExportPID3.resize(total_length3);
1224 for(i=0,j=0; i<2*total_length3; i+=2) {
1225 if(ReverseRecvBuffer[i] != MyPID){
1226 ExportPID3[j]=ReverseRecvBuffer[i];
1227 ExportGID3[j]=ReverseRecvBuffer[i+1];
1233 if(total_length3==0)
return 0;
1236 build_type3_exports_sort<int_type>(ExportGID3, ExportPID3, total_length3);
1237 int StartCurrent, StartNext;
1238 StartCurrent = 0; StartNext = 1;
1239 while ( StartNext < total_length3 ) {
1240 if(ExportPID3[StartNext] == ExportPID3[StartNext-1]) StartNext++;
1242 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1243 StartCurrent = StartNext; StartNext++;
1246 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1257 for(i=1,j=1; i<total_length3; i++){
1258 if(ExportPID3[i]!=ExportPID3[i-1] || ExportGID3[i]!=ExportGID3[i-1]){
1259 ExportPID3[j] = ExportPID3[i];
1260 ExportGID3[j] = ExportGID3[i];
1264 ExportPID3.resize(j);
1265 ExportLID3.resize(j);
1277 for(i=0; i<total_length3; i++) {
1278 ExportLID3[i]=DomainMap.LID(ExportGID3[i]);
1279 if(ExportLID3[i] < 0)
throw std::runtime_error(
"LightweightCrsMatrix:MakeExportLists invalid LID");
1288 return total_length3;
1292 template<
typename ImportType,
typename int_type>
1293 int build_type2_exports(
const Epetra_CrsMatrix & SourceMatrix, ImportType & MyImporter, std::vector<int> &ExportLID2, std::vector<int> &ExportPID2){
1294 int total_length2=0;
1297 const Epetra_Import *SourceImporter= SourceMatrix.Importer();
1299 int *rowptr, *colind;
1301 EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1304 int NumExportIDs = MyImporter.NumExportIDs();
1305 const int* ExportLIDs = MyImporter.ExportLIDs();
1306 const int* ExportPIDs = MyImporter.ExportPIDs();
1307 if(NumExportIDs==0)
return 0;
1309 Epetra_Distributor& Distor = MyImporter.Distributor();
1310 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1314 std::vector<bool> IsOwned(SourceMatrix.NumMyCols(),
true);
1315 if(SourceImporter) {
1316 const int * RemoteLIDs = SourceImporter->RemoteLIDs();
1318 for(i=0; i<SourceImporter->NumRemoteIDs(); i++)
1319 IsOwned[RemoteLIDs[i]]=
false;
1323 std::vector<int> SentTo(SourceMatrix.NumMyCols(),-1);
1326 for(i=0,total_length2=0; i<MDistor->NumSends(); i++) total_length2+= MDistor->LengthsTo()[i] * SourceMatrix.MaxNumEntries();
1327 std::vector<int_type> ExportGID2(total_length2);
1329 ExportLID2.resize(total_length2);
1330 ExportPID2.resize(total_length2);
1332 int current=0, last_start=0, last_pid=ExportPIDs[0];
1333 for(i=0; i<NumExportIDs; i++){
1335 int row=ExportLIDs[i];
1336 int pid=ExportPIDs[i];
1338 if(i!=0 && pid>last_pid) {
1340 if(current!=last_start){
1341 int *lids = &ExportLID2[last_start];
1342 Epetra_Util::Sort(
true,current-last_start,&ExportGID2[last_start],0,0,1,&lids,0,0);
1349 else if(pid < last_pid) {
1350 throw std::runtime_error(
"build_type2_exports: ExportPIDs are not sorted!");
1353 for(j=rowptr[row]; j<rowptr[row+1]; j++) {
1356 if(IsOwned[col] && SentTo[col]!=pid){
1358 if(current>= total_length2)
throw std::runtime_error(
"build_type2_exports: More export ids than I thought!");
1360 ExportGID2[current] = (int_type) SourceMatrix.GCID64(col);
1361 ExportLID2[current] = SourceMatrix.DomainMap().LID(ExportGID2[current]);
1362 ExportPID2[current] = pid;
1369 int *lids = ExportLID2.size() > (std::size_t) last_start ? &ExportLID2[last_start] : 0;
1370 Epetra_Util::Sort(
true,current-last_start,ExportGID2.size() > (std::size_t) last_start ? &ExportGID2[last_start] : 0,0,0,1,&lids,0,0);
1373 total_length2=current;
1374 ExportLID2.resize(total_length2);
1375 ExportPID2.resize(total_length2);
1377 return total_length2;
1381 template<
typename int_type>
1382 void build_type1_exports_sort(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int_type>& ExportGID1,
int total_length1);
1386 int * companion[2] = {&ExportLID1[0],&ExportGID1[0]};
1387 Epetra_Util::Sort(
true,total_length1,&ExportPID1[0],0,0,2,&companion[0],0,0);
1392 int * companion = &ExportLID1[0];
1393 long long * companion64 = &ExportGID1[0];
1394 Epetra_Util::Sort(
true,total_length1,&ExportPID1[0],0,0,1,&companion,1,&companion64);
1397 template<
typename int_type>
1398 int build_type1_exports(
const Epetra_Import * Importer1, std::vector<int> &ExportLID1, std::vector<int> &ExportPID1){
1399 int i, total_length1=0;
1400 if(!Importer1)
return 0;
1401 total_length1 = Importer1->NumSend();
1402 if(total_length1==0)
return 0;
1404 std::vector<int_type> ExportGID1(total_length1);
1405 ExportLID1.resize(total_length1);
1406 ExportPID1.resize(total_length1);
1407 const int * ExportLID1Base = Importer1->ExportLIDs();
1408 const int * ExportPID1Base = Importer1->ExportPIDs();
1410 for(i=0; i<total_length1; i++){
1411 ExportLID1[i] = ExportLID1Base[i];
1412 ExportPID1[i] = ExportPID1Base[i];
1413 ExportGID1[i] = (int_type) Importer1->SourceMap().GID64(ExportLID1Base[i]);
1417 build_type1_exports_sort<int_type>(ExportLID1, ExportPID1, ExportGID1, total_length1);
1419 int StartCurrent, StartNext;
1420 StartCurrent = 0; StartNext = 1;
1421 while ( StartNext < total_length1 ) {
1422 if(ExportPID1[StartNext] == ExportPID1[StartNext-1]) StartNext++;
1424 int *new_companion = {&ExportLID1[StartCurrent]};
1425 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1426 StartCurrent = StartNext; StartNext++;
1429 int *new_companion = {&ExportLID1[StartCurrent]};
1430 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1431 return total_length1;
1435 template<
typename ImportType,
typename int_type>
1436 int LightweightCrsMatrix::MakeExportLists(
const Epetra_CrsMatrix & SourceMatrix, ImportType & Importer2,
1437 std::vector<int> &ReverseRecvSizes,
const int_type *ReverseRecvBuffer,
1438 std::vector<int> & ExportPIDs, std::vector<int> & ExportLIDs) {
1440 int MyPID = SourceMatrix.Comm().MyPID();
1443 const Epetra_Import *Importer1= SourceMatrix.Importer();
1453 Epetra_MpiDistributor * Distor1 = (Importer1)?(dynamic_cast<Epetra_MpiDistributor*>(&Importer1->Distributor())):0;
1455 int Nsend1 = (Distor1)?(Distor1->NumSends()):0;
1457 std::vector<int> ExportPID3;
1458 std::vector<int> ExportLID3;
1460 std::vector<int> ExportPID2;
1461 std::vector<int> ExportLID2;
1463 std::vector<int> ExportPID1;
1464 std::vector<int> ExportLID1;
1467 int Len1=build_type1_exports<int_type>(Importer1, ExportLID1, ExportPID1);
1468 int Len2=build_type2_exports<ImportType, int_type>(SourceMatrix, Importer2, ExportLID2, ExportPID2);
1469 int Len3=
build_type3_exports(MyPID,Nsend1,DomainMap_,ReverseRecvSizes, ReverseRecvBuffer, ExportLID3, ExportPID3);
1472 #ifdef HAVE_EPETRAEXT_DEBUG 1476 bool test_passed=
true;
1477 for(i=1; i<Len1; i++) {
1478 if(ExportPID1[i] < ExportPID1[i-1] || (ExportPID1[i] == ExportPID1[i-1] && DomainMap_.GID(ExportLID1[i]) < DomainMap_.GID(ExportLID1[i-1])))
1481 SourceMatrix.Comm().Barrier();
1483 printf(
"[%d] Type1 ERRORS = ",SourceMatrix.Comm().MyPID());
1484 for(
int i=0; i<Len1; i++)
1485 printf(
"(%2d,%2d,%2d) ",ExportLID1[i],DomainMap_.GID(ExportLID1[i]),ExportPID1[i]);
1488 throw std::runtime_error(
"Importer1 fails the sanity test");
1491 for(i=1; i<Len2; i++) {
1492 if(ExportPID2[i] < ExportPID2[i-1] || (ExportPID2[i] == ExportPID2[i-1] && DomainMap_.GID(ExportLID2[i]) < DomainMap_.GID(ExportLID2[i-1])))
1496 SourceMatrix.Comm().Barrier();
1498 printf(
"[%d] Type2 ERRORS = ",SourceMatrix.Comm().MyPID());
1499 for(
int i=0; i<Len2; i++)
1500 printf(
"(%2d,%2d,%2d) ",ExportLID2[i],DomainMap_.GID(ExportLID2[i]),ExportPID2[i]);
1503 throw std::runtime_error(
"Importer2 fails the sanity test");
1506 for(i=1; i<Len3; i++) {
1507 if(ExportPID3[i] < ExportPID3[i-1] || (ExportPID3[i] == ExportPID3[i-1] && DomainMap_.GID(ExportLID3[i]) < DomainMap_.GID(ExportLID3[i-1])))
1511 SourceMatrix.Comm().Barrier();
1513 printf(
"[%d] Type3 ERRORS = ",SourceMatrix.Comm().MyPID());
1514 for(
int i=0; i<Len3; i++)
1515 printf(
"(%2d,%2d,%2d) ",ExportLID3[i],DomainMap_.GID(ExportLID3[i]),ExportPID3[i]);
1518 throw std::runtime_error(
"Importer3 fails the sanity test");
1525 if(Importer1 && !Importer1->SourceMap().SameAs(DomainMap_))
1526 throw std::runtime_error(
"ERROR: Map Mismatch Importer1");
1528 if(!Importer2.SourceMap().SameAs(SourceMatrix.RowMap()))
1529 throw std::runtime_error(
"ERROR: Map Mismatch Importer2");
1531 int_type InfGID = std::numeric_limits<int_type>::max();
1532 int InfPID = INT_MAX;
1534 int i1=0, i2=0, i3=0, current=0;
1536 int MyLen=Len1+Len2+Len3;
1537 ExportLIDs.resize(MyLen);
1538 ExportPIDs.resize(MyLen);
1540 while(i1 < Len1 || i2 < Len2 || i3 < Len3){
1541 int PID1 = (i1<Len1)?(ExportPID1[i1]):InfPID;
1542 int PID2 = (i2<Len2)?(ExportPID2[i2]):InfPID;
1543 int PID3 = (i3<Len3)?(ExportPID3[i3]):InfPID;
1545 int_type GID1 = (i1<Len1)?((int_type) DomainMap_.GID64(ExportLID1[i1])):InfGID;
1546 int_type GID2 = (i2<Len2)?((int_type) DomainMap_.GID64(ExportLID2[i2])):InfGID;
1547 int_type GID3 = (i3<Len3)?((int_type) DomainMap_.GID64(ExportLID3[i3])):InfGID;
1549 int MIN_PID =
MIN3(PID1,PID2,PID3);
1550 int_type MIN_GID =
MIN3( ((PID1==MIN_PID)?GID1:InfGID), ((PID2==MIN_PID)?GID2:InfGID), ((PID3==MIN_PID)?GID3:InfGID));
1551 bool added_entry=
false;
1554 if(PID1 == MIN_PID && GID1 == MIN_GID){
1555 ExportLIDs[current] = ExportLID1[i1];
1556 ExportPIDs[current] = ExportPID1[i1];
1563 if(PID2 == MIN_PID && GID2 == MIN_GID){
1565 ExportLIDs[current] = ExportLID2[i2];
1566 ExportPIDs[current] = ExportPID2[i2];
1574 if(PID3 == MIN_PID && GID3 == MIN_GID){
1576 ExportLIDs[current] = ExportLID3[i3];
1577 ExportPIDs[current] = ExportPID3[i3];
1583 if(current!=MyLen) {
1584 ExportLIDs.resize(current);
1585 ExportPIDs.resize(current);
1592 template<
typename ImportType,
typename int_type>
1593 void LightweightCrsMatrix::Construct(
const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter)
1597 #ifdef ENABLE_MMM_TIMINGS 1598 using Teuchos::TimeMonitor;
1599 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS C-1")));
1605 if(use_lw) N = RowMapLW_->NumMyElements();
1606 else N = RowMapEP_->NumMyElements();
1608 int MyPID = SourceMatrix.Comm().MyPID();
1611 std::vector<int> ReverseSendSizes;
1612 std::vector<int_type> ReverseSendBuffer;
1613 std::vector<int> ReverseRecvSizes;
1614 int_type * ReverseRecvBuffer=0;
1617 bool communication_needed = RowImporter.SourceMap().DistributedGlobal();
1626 if (!SourceMatrix.RowMap().SameAs(RowImporter.SourceMap()))
1627 throw "LightweightCrsMatrix: Fused copy constructor requires Importer.SourceMap() to match SourceMatrix.RowMap()";
1630 int NumSameIDs = RowImporter.NumSameIDs();
1631 int NumPermuteIDs = RowImporter.NumPermuteIDs();
1632 int NumRemoteIDs = RowImporter.NumRemoteIDs();
1633 int NumExportIDs = RowImporter.NumExportIDs();
1634 int* ExportLIDs = RowImporter.ExportLIDs();
1635 int* RemoteLIDs = RowImporter.RemoteLIDs();
1636 int* PermuteToLIDs = RowImporter.PermuteToLIDs();
1637 int* PermuteFromLIDs = RowImporter.PermuteFromLIDs();
1640 Epetra_Distributor& Distor = RowImporter.Distributor();
1641 const Epetra_MpiComm * MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&SourceMatrix.Comm());
1642 const Epetra_MpiDistributor * MDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&Distor);
1646 rowptr_.resize(N+1);
1650 std::vector<int> SourcePids;
1651 std::vector<int> TargetPids;
1660 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in CheckSizes()";
1666 int LenExports_ = 0;
1670 bool VarSizes =
false;
1671 if( NumExportIDs > 0) {
1672 Sizes_ =
new int[NumExportIDs];
1677 const Epetra_Import *MyImporter= SourceMatrix.Importer();
1678 if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,
false);
1680 SourcePids.resize(SourceMatrix.ColMap().NumMyElements());
1681 SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),SourceMatrix.Comm().MyPID());
1685 rv=Epetra_Import_Util::PackAndPrepareWithOwningPIDs(SourceMatrix,
1686 NumExportIDs,ExportLIDs,
1687 LenExports_,Exports_,SizeOfPacket,
1688 Sizes_,VarSizes,SourcePids);
1689 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in PackAndPrepare()";
1691 if (communication_needed) {
1694 const int * ExportPIDs = RowImporter.ExportPIDs();
1698 std::vector<int> SendSizes(MDistor->NumSends()+1,0);
1699 for(
int i=0, curr_pid=0; i<NumExportIDs; i++) {
1700 if(i>0 && ExportPIDs[i] > ExportPIDs[i-1]) curr_pid++;
1701 SendSizes[curr_pid] +=Sizes_[i];
1704 if(i>0 && ExportPIDs[i] < ExportPIDs[i-1])
throw "ExportPIDs not sorted";
1707 std::vector<int> RecvSizes(MDistor->NumReceives()+1);
1708 int msg_tag=MpiComm->GetMpiTag();
1709 boundary_exchange_varsize<char>(*MpiComm,MPI_CHAR,MDistor->NumSends(),MDistor->ProcsTo(),SendSizes.size() ? &SendSizes[0] : 0,Exports_,
1710 MDistor->NumReceives(),MDistor->ProcsFrom(),RecvSizes.size() ? &RecvSizes[0] : 0,Imports_,SizeOfPacket,msg_tag);
1713 if(SourceMatrix.Importer()) {
1714 Epetra_Import* SourceImporter=
const_cast<Epetra_Import*
>(SourceMatrix.Importer());
1715 const Epetra_MpiDistributor * MyDistor =
dynamic_cast<Epetra_MpiDistributor*
>(&SourceImporter->Distributor());
1719 PackAndPrepareReverseComm<ImportType, int_type>(SourceMatrix,RowImporter,ReverseSendSizes,ReverseSendBuffer);
1723 ReverseRecvSizes.resize(MyDistor->NumSends()+1);
1724 const int msg_tag2 = MpiComm->GetMpiTag ();
1725 MPI_Datatype data_type =
sizeof(int_type) == 4 ? MPI_INT : MPI_LONG_LONG;
1726 boundary_exchange_varsize<int_type> (*MpiComm, data_type, MyDistor->NumReceives (),
1727 MyDistor->ProcsFrom (),
1728 ReverseSendSizes.size() ? &ReverseSendSizes[0] : 0,
1729 ReverseSendBuffer.size() ? &ReverseSendBuffer[0] : 0,
1730 MyDistor->NumSends (), MyDistor->ProcsTo (),
1731 ReverseRecvSizes.size() ? &ReverseRecvSizes[0] : 0,
1732 ReverseRecvBuffer, 1, msg_tag2);
1737 if(rv)
throw "LightweightCrsMatrix: Fused copy constructor failed in Distor.Do";
1742 #ifdef ENABLE_MMM_TIMINGS 1743 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS C-2")));
1747 int mynnz = Epetra_Import_Util::UnpackWithOwningPIDsCount(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_);
1750 colind_.resize(mynnz);
1751 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1752 if(
sizeof(int_type) ==
sizeof(
long long))
1753 colind_LL_.resize(mynnz);
1755 vals_.resize(mynnz);
1758 if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,
true);
1760 SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),-1);
1764 int * myrowptr = rowptr_.size() ? & rowptr_[0] : 0;
1765 double * myvals = vals_.size() ? & vals_[0] : 0;
1766 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1767 if(
sizeof(int_type) ==
sizeof(
int)) {
1768 int * mycolind = colind_.size() ? & colind_[0] : 0;
1769 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1773 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1774 if(
sizeof(int_type) ==
sizeof(
long long)) {
1775 long long * mycolind = colind_LL_.size() ? & colind_LL_[0] : 0;
1776 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1780 throw "EpetraExt::LightweightCrsMatrix::Construct: sizeof(int_type) error.";
1785 #ifdef ENABLE_MMM_TIMINGS 1786 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS C-3")));
1790 MakeColMapAndReindex<int_type>(TargetPids,getcolind<int_type>());
1796 MakeExportLists<ImportType, int_type>(SourceMatrix,RowImporter,ReverseRecvSizes,ReverseRecvBuffer,ExportPIDs_,ExportLIDs_);
1803 #ifdef ENABLE_MMM_TIMINGS 1804 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS C-4")));
1806 if(N>0) Epetra_Util::SortCrsEntries(N, &rowptr_[0], colind_.size() ? &colind_[0] : 0, vals_.size() ? &vals_[0] : 0);
1811 #ifdef ENABLE_MMM_TIMINGS 1812 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS C-5")));
1816 delete [] ReverseRecvBuffer;
1833 DomainMap_(SourceMatrix.DomainMap())
1835 #ifdef ENABLE_MMM_TIMINGS 1836 using Teuchos::TimeMonitor;
1837 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: LWCRS Total")));
1842 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1843 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1844 Construct<RemoteOnlyImport, int>(SourceMatrix,RowImporter);
1848 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1849 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1850 Construct<RemoteOnlyImport, long long>(SourceMatrix,RowImporter);
1854 throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1866 RowMapEP_=
new Epetra_BlockMap(RowImporter.TargetMap());
1867 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1868 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1869 Construct<Epetra_Import, int>(SourceMatrix,RowImporter);
1873 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1874 if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1875 Construct<Epetra_Import, long long>(SourceMatrix,RowImporter);
1879 throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1891 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1895 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
LightweightCrsMatrix * importMatrix
RemoteOnlyImport(const Epetra_Import &Importer, LightweightMap &RemoteOnlyTargetMap)
Epetra_HashTable< long long > * LIDHash_LL_
Epetra_BlockMap * RowMapEP_
std::vector< int > targetMapToOrigRow
long long * MyGlobalElements64() const
virtual ~CrsWrapper_GraphBuilder()
void debug_print_distor(const char *label, const Epetra_Distributor *Distor, const Epetra_Comm &Comm)
std::pair< int_type, int_type > Tget_col_range(const Epetra_CrsMatrix &mtx)
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
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 MakeColMapAndReindexSort< long long >(int &NumRemoteColGIDs, long long *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
void debug_compare_import(const Epetra_Import *Import1, const Epetra_Import *Import2)
const Epetra_Map * rowMap
long long GID64(int LID) const
void MakeColMapAndReindexSort< int >(int &NumRemoteColGIDs, int *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
const Epetra_CrsMatrix * origMatrix
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
int NumMyElements() const
LightweightMap & operator=(const LightweightMap &map)
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 build_type2_exports(const Epetra_CrsMatrix &SourceMatrix, ImportType &MyImporter, std::vector< int > &ExportLID2, std::vector< int > &ExportPID2)
int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
const Epetra_Map & RowMap() const
LightweightCrsMatrix(const Epetra_CrsMatrix &A, RemoteOnlyImport &RowImporter)
CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix &epetracrsmatrix)
std::vector< int > MyGlobalElements_int_
int * MyGlobalElements() const
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)
const LightweightMap & TargetMap() const
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)
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)
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)
const Epetra_Map * colMap
const Epetra_BlockMap * importColMap
void insert_matrix_locations(CrsWrapper_GraphBuilder< int_type > &graphbuilder, Epetra_CrsMatrix &C)
void MakeColMapAndReindexSort(int &NumRemoteColGIDs, GO *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
int dumpCrsMatrixStruct(const CrsMatrixStruct &M)
int InsertGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
std::pair< int, int > get_col_range< int >(const Epetra_Map &emap)