61 template<
typename int_type>
70 std::vector<int>& pids)
74 int TotalSendLength = 0;
76 if( NumExportIDs>0 ) IntSizes =
new int[NumExportIDs];
78 int SizeofIntType =
sizeof(int_type);
80 for(
int i = 0; i < NumExportIDs; ++i) {
82 A.NumMyRowEntries( ExportLIDs[i], NumEntries );
85 Sizes[i] = NumEntries;
87 IntSizes[i] = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
88 TotalSendLength += (Sizes[i]+IntSizes[i]);
91 double * DoubleExports = 0;
92 SizeOfPacket = (int)
sizeof(
double);
95 if( TotalSendLength*SizeOfPacket > LenExports ) {
96 if( LenExports > 0 )
delete [] Exports;
97 LenExports = TotalSendLength*SizeOfPacket;
98 DoubleExports =
new double[TotalSendLength];
99 for(
int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
100 Exports = (
char *) DoubleExports;
105 double * valptr, * dintptr;
117 if (NumExportIDs > 0) {
122 int maxNumEntries =
A.MaxNumEntries();
123 std::vector<int> MyIndices(maxNumEntries);
126 dintptr = (
double *) Exports;
127 valptr = dintptr + IntSizes[0];
130 intptr = (int_type *) dintptr;
131 for (
int i = 0; i < NumExportIDs; ++i) {
132 FromRow = (int_type) rowMap.
GID64(ExportLIDs[i]);
135 Indices = intptr + 2;
136 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, MyIndices.size() ? &MyIndices[0] : 0));
137 for (
int j = 0; j < NumEntries; ++j) {
138 Indices[2*j] = (int_type) colMap.
GID64(MyIndices[j]);
139 Indices[2*j+1] = pids[MyIndices[j]];
141 intptr[1] = NumEntries;
142 if (i < (NumExportIDs-1)) {
143 dintptr += (IntSizes[i]+Sizes[i]);
144 valptr = dintptr + IntSizes[i+1];
147 intptr = (int_type *) dintptr;
151 for (
int i = 0; i < NumExportIDs; ++i) {
152 Sizes[i] += IntSizes[i];
156 if( IntSizes )
delete [] IntSizes;
172 std::vector<int>& SourcePids)
175 return TPackAndPrepareWithOwningPIDs<int>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
178 return TPackAndPrepareWithOwningPIDs<long long>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
181 throw std::runtime_error(
"UnpackWithOwningPIDsCount: Unable to determine source global index type");
187 template<
typename int_type>
191 const int * RemoteLIDs,
193 const int *PermuteToLIDs,
194 const int *PermuteFromLIDs,
199 int SizeofIntType = (int)
sizeof(int_type);
202 for(i=0; i<NumSameIDs; i++)
206 for(i=0; i<NumPermuteIDs; i++)
210 if(NumRemoteIDs > 0) {
211 double * dintptr = (
double *) Imports;
213 int_type * intptr = (int_type *) dintptr;
214 int NumEntries = (int) intptr[1];
215 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
216 for(i=0; i<NumRemoteIDs; i++) {
219 if( i < (NumRemoteIDs-1) ) {
220 dintptr += IntSize + NumEntries;
221 intptr = (int_type *) dintptr;
222 NumEntries = (int) intptr[1];
223 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
237 const int * RemoteLIDs,
239 const int *PermuteToLIDs,
240 const int *PermuteFromLIDs,
245 return TUnpackWithOwningPIDsCount<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
246 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
249 return TUnpackWithOwningPIDsCount<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
250 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
253 throw std::runtime_error(
"UnpackWithOwningPIDsCount: Unable to determine source global index type");
260 template<
typename int_type>
264 const int * RemoteLIDs,
266 const int *PermuteToLIDs,
267 const int *PermuteFromLIDs,
271 int TargetNumNonzeros,
274 int_type * CSR_colind,
276 const std::vector<int> &SourcePids,
277 std::vector<int> &TargetPids)
290 int N = TargetNumRows;
291 int mynnz = TargetNumNonzeros;
293 int MyPID = MyTargetPID;
295 int SizeofIntType =
sizeof(int_type);
298 for(i=0; i<N+1; i++) CSR_rowptr[i]=0;
301 for(i=0; i<NumSameIDs; i++)
305 for(i=0; i<NumPermuteIDs; i++)
306 CSR_rowptr[PermuteToLIDs[i]] = SourceMatrix.
NumMyEntries(PermuteFromLIDs[i]);
309 if(NumRemoteIDs > 0) {
310 double * dintptr = (
double *) Imports;
311 int_type * intptr = (int_type *) dintptr;
312 int NumEntries = (int) intptr[1];
313 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
314 for(i=0; i<NumRemoteIDs; i++) {
315 CSR_rowptr[RemoteLIDs[i]] += NumEntries;
317 if( i < (NumRemoteIDs-1) ) {
318 dintptr += IntSize + NumEntries;
319 intptr = (int_type *) dintptr;
320 NumEntries = (int) intptr[1];
321 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
327 std::vector<int> NewStartRow(N+1);
330 int last_len = CSR_rowptr[0];
332 for(i=1; i<N+1; i++){
333 int new_len = CSR_rowptr[i];
334 CSR_rowptr[i] = last_len + CSR_rowptr[i-1];
335 NewStartRow[i] = CSR_rowptr[i];
340 if(TargetPids.size()!=(size_t)mynnz) TargetPids.resize(mynnz);
341 TargetPids.assign(mynnz,-1);
344 int * Source_rowptr, * Source_colind;
345 double * Source_vals;
347 if(rv)
throw std::runtime_error(
"UnpackAndCombineIntoCrsArrays: failed in ExtractCrsDataPointers");
350 for(i=0; i<NumSameIDs; i++) {
351 int FromRow = Source_rowptr[i];
352 int ToRow = CSR_rowptr[i];
353 NewStartRow[i] += Source_rowptr[i+1]-Source_rowptr[i];
355 for(j=Source_rowptr[i]; j<Source_rowptr[i+1]; j++) {
356 CSR_vals[ToRow + j - FromRow] = Source_vals[j];
357 CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.
GCID64(Source_colind[j]);
358 TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
363 for(i=0; i<NumPermuteIDs; i++) {
364 int FromLID = PermuteFromLIDs[i];
365 int FromRow = Source_rowptr[FromLID];
366 int ToRow = CSR_rowptr[PermuteToLIDs[i]];
368 NewStartRow[PermuteToLIDs[i]] += Source_rowptr[FromLID+1]-Source_rowptr[FromLID];
370 for(j=Source_rowptr[FromLID]; j<Source_rowptr[FromLID+1]; j++) {
371 CSR_vals[ToRow + j - FromRow] = Source_vals[j];
372 CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.
GCID64(Source_colind[j]);
373 TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
378 if(NumRemoteIDs > 0) {
379 double * dintptr = (
double *) Imports;
380 int_type * intptr = (int_type *) dintptr;
381 int NumEntries = (int) intptr[1];
382 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
383 double* valptr = dintptr + IntSize;
385 for (i=0; i<NumRemoteIDs; i++) {
386 int ToLID = RemoteLIDs[i];
387 int StartRow = NewStartRow[ToLID];
388 NewStartRow[ToLID]+=NumEntries;
390 double * values = valptr;
391 int_type * Indices = intptr + 2;
392 for(j=0; j<NumEntries; j++){
393 CSR_vals[StartRow + j] = values[j];
394 CSR_colind[StartRow + j] = Indices[2*j];
395 TargetPids[StartRow + j] = (Indices[2*j+1] != MyPID) ? Indices[2*j+1] : -1;
397 if( i < (NumRemoteIDs-1) ) {
398 dintptr += IntSize + NumEntries;
399 intptr = (int_type *) dintptr;
400 NumEntries = (int) intptr[1];
401 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
402 valptr = dintptr + IntSize;
417 const int * RemoteLIDs,
419 const int *PermuteToLIDs,
420 const int *PermuteFromLIDs,
424 int TargetNumNonzeros,
429 const std::vector<int> &SourcePids,
430 std::vector<int> &TargetPids)
433 return TUnpackAndCombineIntoCrsArrays<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
434 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
435 CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
438 throw std::runtime_error(
"UnpackAndCombineIntoCrsArrays: int type not matched.");
446 const int * RemoteLIDs,
448 const int *PermuteToLIDs,
449 const int *PermuteFromLIDs,
453 int TargetNumNonzeros,
456 long long * CSR_colind,
458 const std::vector<int> &SourcePids,
459 std::vector<int> &TargetPids)
462 return TUnpackAndCombineIntoCrsArrays<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
463 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
464 CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
467 throw std::runtime_error(
"UnpackAndCombineIntoCrsArrays: int type not matched.");
473 template<
typename int_type,
class MapType1,
class MapType2>
474 int TLowCommunicationMakeColMapAndReindex(
int N,
const int * rowptr,
int * colind_LID,
const int_type *colind_GID,
const Epetra_Map& domainMap,
const int * owningPIDs,
bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, MapType1 & NewColMap)
483 else throw std::runtime_error(
"LowCommunicationMakeColMapAndReindex: cannot detect int type.");
489 bool * LocalGIDs = 0;
490 if (numDomainElements>0) LocalGIDs =
new bool[numDomainElements];
491 for (i=0; i<numDomainElements; i++) LocalGIDs[i] =
false;
494 if(DoSizes)
throw std::runtime_error(
"LowCommunicationMakeColMapAndReindex: cannot handle non-constant sized domainMap.");
500 const int numMyBlockRows = N;
501 int hashsize = numMyBlockRows;
if (hashsize < 100) hashsize = 100;
503 std::vector<int_type> RemoteGIDList; RemoteGIDList.reserve(hashsize);
504 std::vector<int> PIDList; PIDList.reserve(hashsize);
513 int NumLocalColGIDs = 0;
514 int NumRemoteColGIDs = 0;
515 for(i = 0; i < numMyBlockRows; i++) {
516 for(j = rowptr[i]; j < rowptr[i+1]; j++) {
517 int_type GID = colind_GID[j];
519 int LID = domainMap.
LID(GID);
521 bool alreadyFound = LocalGIDs[LID];
523 LocalGIDs[LID] =
true;
529 int_type hash_value=RemoteGIDs.
Get(GID);
530 if(hash_value == -1) {
531 int PID = owningPIDs[j];
532 if(PID==-1)
throw std::runtime_error(
"LowCommunicationMakeColMapAndReindex: Cannot figure out if PID is owned.");
533 colind_LID[j] = numDomainElements + NumRemoteColGIDs;
534 RemoteGIDs.
Add(GID, NumRemoteColGIDs);
535 RemoteGIDList.push_back(GID);
536 PIDList.push_back(PID);
540 colind_LID[j] = numDomainElements + hash_value;
548 if (NumRemoteColGIDs!=0) {
549 throw std::runtime_error(
"Some column IDs are not in domainMap. If matrix is rectangular, you must pass in a domainMap");
552 if (NumLocalColGIDs==numDomainElements) {
553 if (LocalGIDs!=0)
delete [] LocalGIDs;
556 NewColMap = domainMap;
563 int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
564 std::vector<int_type> ColIndices;
565 int_type * RemoteColIndices=0;
566 if(numMyBlockCols > 0) {
567 ColIndices.resize(numMyBlockCols);
568 if(NumLocalColGIDs!=numMyBlockCols) RemoteColIndices = &ColIndices[NumLocalColGIDs];
569 else RemoteColIndices=0;
572 for(i = 0; i < NumRemoteColGIDs; i++)
573 RemoteColIndices[i] = RemoteGIDList[i];
576 std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
577 for(i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
582 int * IntSortLists[2];
583 long long * LLSortLists[2];
584 int * RemotePermuteIDs_ptr = RemotePermuteIDs.size() ? &RemotePermuteIDs[0] : 0;
587 IntSortLists[0] = (
int*) RemoteColIndices;
588 IntSortLists[1] = RemotePermuteIDs_ptr;
593 LLSortLists[0] = (
long long*) RemoteColIndices;
594 IntSortLists[0] = RemotePermuteIDs_ptr;
595 NumListsInt = NumListsLL = 1;
598 int * PIDList_ptr = PIDList.size() ? &PIDList[0] : 0;
599 Epetra_Util::Sort(
true, NumRemoteColGIDs, PIDList_ptr, 0, 0, NumListsInt, IntSortLists,NumListsLL,LLSortLists);
602 PIDList.resize(NumRemoteColGIDs);
603 RemotePIDs = PIDList;
605 if (SortGhostsAssociatedWithEachProcessor) {
611 int StartCurrent, StartNext;
612 StartCurrent = 0; StartNext = 1;
613 while ( StartNext < NumRemoteColGIDs ) {
614 if (PIDList[StartNext]==PIDList[StartNext-1]) StartNext++;
616 IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
617 Epetra_Util::Sort(
true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,1,IntSortLists,0,0);
618 StartCurrent = StartNext; StartNext++;
621 IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
622 Epetra_Util::Sort(
true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, 1,IntSortLists,0,0);
626 std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
627 for(i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
630 bool use_local_permute=
false;
631 std::vector<int> LocalPermuteIDs(numDomainElements);
640 if(NumLocalColGIDs > 0) {
645 int_type* MyGlobalElements = 0;
648 int NumLocalAgain = 0;
649 use_local_permute =
true;
650 for(i = 0; i < numDomainElements; i++) {
652 LocalPermuteIDs[i] = NumLocalAgain;
653 ColIndices[NumLocalAgain++] = MyGlobalElements[i];
656 assert(NumLocalAgain==NumLocalColGIDs);
660 if (LocalGIDs!=0)
delete [] LocalGIDs;
663 int_type * ColIndices_ptr = ColIndices.size() ? &ColIndices[0] : 0;
664 MapType2 temp((int_type)(-1), numMyBlockCols, ColIndices_ptr, (int_type)domainMap.
IndexBase64(), domainMap.
Comm());
668 for(i=0; i<numMyBlockRows; i++){
669 for(j=rowptr[i]; j<rowptr[i+1]; j++){
670 int ID=colind_LID[j];
671 if(ID < numDomainElements){
672 if(use_local_permute) colind_LID[j] = LocalPermuteIDs[colind_LID[j]];
677 colind_LID[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_LID[j]-numDomainElements];
687 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 692 return TLowCommunicationMakeColMapAndReindex<int,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind,colind,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
694 throw std::runtime_error(
"LowCommunicationMakeColMapAndReindex: GID type mismatch.");
701 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 706 return TLowCommunicationMakeColMapAndReindex<long long,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind_LID,colind_GID,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
708 throw std::runtime_error(
"LowCommunicationMakeColMapAndReindex: GID type mismatch.");
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
int TLowCommunicationMakeColMapAndReindex(int N, const int *rowptr, int *colind_LID, const int_type *colind_GID, const Epetra_Map &domainMap, const int *owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector< int > &RemotePIDs, MapType1 &NewColMap)
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Returns internal data pointers associated with Crs matrix format.
bool ConstantElementSize() const
Returns true if map has constant element size.
value_type Get(const long long key)
#define EPETRA_CHK_ERR(a)
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int PackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix &SourceMatrix, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, std::vector< int > &SourcePids)
PackAndPrepareWithOwningPIDs.
int TUnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports, int TargetNumRows, int TargetNumNonzeros, int MyTargetPID, int *CSR_rowptr, int_type *CSR_colind, double *CSR_vals, const std::vector< int > &SourcePids, std::vector< int > &TargetPids)
int UnpackWithOwningPIDsCount(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports)
UnpackWithOwningPIDsCount.
int UnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports, int TargetNumRows, int TargetNumNonzeros, int MyTargetPID, int *CSR_rowptr, int *CSR_colind, double *CSR_values, const std::vector< int > &SourcePids, std::vector< int > &TargetPids)
UnpackAndCombineIntoCrsArrays.
Epetra_Import_Util: The Epetra ImportUtil Wrapper Namespace.
long long GID64(int LID) const
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
int NumMyElements() const
Number of elements on the calling processor.
int TUnpackWithOwningPIDsCount(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int LowCommunicationMakeColMapAndReindex(int N, const int *rowptr, int *colind, const Epetra_Map &domainMap, const int *owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector< int > &RemotePIDs, Epetra_BlockMap &NewColMap)
LowCommunicationMakeColMapAndReindex.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
long long IndexBase64() const
virtual int NumProc() const =0
Returns total number of processes.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
static void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_Util Sort Routine (Shell sort)
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
int TPackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, std::vector< int > &pids)
long long GCID64(int LCID_in) const
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
void Add(const long long key, const value_type value)