Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_Import_Util.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
44#include "Epetra_ConfigDefs.h"
45#include "Epetra_Import_Util.h"
46#include "Epetra_Util.h"
47#include "Epetra_Comm.h"
48#include "Epetra_BlockMap.h"
49#include "Epetra_Map.h"
50#include "Epetra_Import.h"
51#include "Epetra_CrsMatrix.h"
52#include "Epetra_HashTable.h"
53#include "Epetra_Util.h"
54
55#include <stdexcept>
56
57
59
60// =========================================================================
61// =========================================================================
62template<typename int_type>
64 int NumExportIDs,
65 int * ExportLIDs,
66 int & LenExports,
67 char *& Exports,
68 int & SizeOfPacket,
69 int * Sizes,
70 bool & VarSizes,
71 std::vector<int>& pids)
72{
73 VarSizes = true; //enable variable block size data comm
74
75 int TotalSendLength = 0;
76 int * IntSizes = 0;
77 if( NumExportIDs>0 ) IntSizes = new int[NumExportIDs];
78
79 int SizeofIntType = sizeof(int_type);
80
81 for(int i = 0; i < NumExportIDs; ++i) {
82 int NumEntries;
83 A.NumMyRowEntries( ExportLIDs[i], NumEntries );
84 // Will have NumEntries doubles, 2*NumEntries +2 ints pack them interleaved Sizes[i] = NumEntries;
85 // NTS: We add the owning PID as the SECOND int of the pair for each entry
86 Sizes[i] = NumEntries;
87 // NOTE: Mixing and matching Int Types would be more efficient, BUT what about variable alignment?
88 IntSizes[i] = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
89 TotalSendLength += (Sizes[i]+IntSizes[i]);
90 }
91
92 double * DoubleExports = 0;
93 SizeOfPacket = (int)sizeof(double);
94
95 //setup buffer locally for memory management by this object
96 if( TotalSendLength*SizeOfPacket > LenExports ) {
97 if( LenExports > 0 ) delete [] Exports;
98 LenExports = TotalSendLength*SizeOfPacket;
99 DoubleExports = new double[TotalSendLength];
100 for( int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
101 Exports = (char *) DoubleExports;
102 }
103
104 int NumEntries;
105 double * values;
106 double * valptr, * dintptr;
107
108 // Each segment of Exports will be filled by a packed row of information for each row as follows:
109 // 1st int: GRID of row where GRID is the global row ID for the source matrix
110 // next int: NumEntries, Number of indices in row
111 // next 2*NumEntries: The actual indices and owning [1] PID each for the row in (GID,PID) pairs with the GID first.
112
113 // [1] Owning is defined in the sense of "Who owns the GID in the DomainMap," aka, who sends the GID in the importer
114
115 const Epetra_Map & rowMap = A.RowMap();
116 const Epetra_Map & colMap = A.ColMap();
117
118 if (NumExportIDs > 0) {
119 int_type * Indices;
120 int_type FromRow;
121 int_type * intptr;
122
123 int maxNumEntries = A.MaxNumEntries();
124 std::vector<int> MyIndices(maxNumEntries);
125 // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
126 // prohibition of unaligned type-punning access.
127 dintptr = (double *) Exports;
128 valptr = dintptr + IntSizes[0];
129 // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
130 // prohibition of unaligned type-punning access.
131 intptr = (int_type *) dintptr;
132 for (int i = 0; i < NumExportIDs; ++i) {
133 FromRow = (int_type) rowMap.GID64(ExportLIDs[i]);
134 intptr[0] = FromRow;
135 values = valptr;
136 Indices = intptr + 2;
137 EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, Epetra_Util_data_ptr(MyIndices)));
138 for (int j = 0; j < NumEntries; ++j) {
139 Indices[2*j] = (int_type) colMap.GID64(MyIndices[j]); // convert to GIDs
140 Indices[2*j+1] = pids[MyIndices[j]]; // PID owning the entry.
141 }
142 intptr[1] = NumEntries; // Load second slot of segment
143 if (i < (NumExportIDs-1)) {
144 dintptr += (IntSizes[i]+Sizes[i]);
145 valptr = dintptr + IntSizes[i+1];
146 // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
147 // prohibition of unaligned type-punning access.
148 intptr = (int_type *) dintptr;
149 }
150 }
151
152 for (int i = 0; i < NumExportIDs; ++i) {
153 Sizes[i] += IntSizes[i];
154 }
155 }
156
157 if( IntSizes ) delete [] IntSizes;
158
159 return 0;
160}
161
162
163// =========================================================================
164// =========================================================================
166 int NumExportIDs,
167 int * ExportLIDs,
168 int & LenExports,
169 char *& Exports,
170 int & SizeOfPacket,
171 int * Sizes,
172 bool & VarSizes,
173 std::vector<int>& SourcePids)
174{
175 if(SourceMatrix.RowMap().GlobalIndicesInt()) {
176 return TPackAndPrepareWithOwningPIDs<int>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
177 }
178 else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
179 return TPackAndPrepareWithOwningPIDs<long long>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
180 }
181 else
182 throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
183}
184
185
186// =========================================================================
187// =========================================================================
188template<typename int_type>
190 int NumSameIDs,
191 int NumRemoteIDs,
192 const int * /* RemoteLIDs */,
193 int NumPermuteIDs,
194 const int */* PermuteToLIDs */,
195 const int *PermuteFromLIDs,
196 int /* LenImports */,
197 char* Imports)
198{
199 int i,nnz=0;
200 int SizeofIntType = (int)sizeof(int_type);
201
202 // SameIDs: Always first, always in the same place
203 for(i=0; i<NumSameIDs; i++)
204 nnz+=SourceMatrix.NumMyEntries(i);
205
206 // PermuteIDs: Still local, but reordered
207 for(i=0; i<NumPermuteIDs; i++)
208 nnz += SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
209
210 // RemoteIDs: RemoteLIDs tells us the ID, we need to look up the length the hard way. See UnpackAndCombine for where this code came from
211 if(NumRemoteIDs > 0) {
212 double * dintptr = (double *) Imports;
213 // General version
214 int_type * intptr = (int_type *) dintptr;
215 int NumEntries = (int) intptr[1];
216 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
217 for(i=0; i<NumRemoteIDs; i++) {
218 nnz += NumEntries;
219
220 if( i < (NumRemoteIDs-1) ) {
221 dintptr += IntSize + NumEntries;
222 intptr = (int_type *) dintptr;
223 NumEntries = (int) intptr[1];
224 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
225 }
226 }
227 }
228
229 return nnz;
230}
231
232
233// =========================================================================
234// =========================================================================
236 int NumSameIDs,
237 int NumRemoteIDs,
238 const int * RemoteLIDs,
239 int NumPermuteIDs,
240 const int *PermuteToLIDs,
241 const int *PermuteFromLIDs,
242 int LenImports,
243 char* Imports)
244{
245 if(SourceMatrix.RowMap().GlobalIndicesInt()) {
246 return TUnpackWithOwningPIDsCount<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
247 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
248 }
249 else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
250 return TUnpackWithOwningPIDsCount<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
251 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
252 }
253 else
254 throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
255
256}
257
258
259// =========================================================================
260// =========================================================================
261template<typename int_type>
263 int NumSameIDs,
264 int NumRemoteIDs,
265 const int * RemoteLIDs,
266 int NumPermuteIDs,
267 const int *PermuteToLIDs,
268 const int *PermuteFromLIDs,
269 int /* LenImports */,
270 char* Imports,
271 int TargetNumRows,
272 int TargetNumNonzeros,
273 int MyTargetPID,
274 int * CSR_rowptr,
275 int_type * CSR_colind,
276 double * CSR_vals,
277 const std::vector<int> &SourcePids,
278 std::vector<int> &TargetPids)
279{
280 // What we really need to know is where in the CSR arrays each row should start (aka the rowptr).
281 // We do that by (a) having each row record it's size in the rowptr (b) doing a cumulative sum to get the rowptr values correct and
282 // (c) Having each row copied into the right colind / values locations.
283
284 // From Epetra_CrsMatrix UnpackAndCombine():
285 // Each segment of Exports will be filled by a packed row of information for each row as follows:
286 // 1st int: GRID of row where GRID is the global row ID for the source matrix
287 // next int: NumEntries, Number of indices in row.
288 // next NumEntries: The actual indices for the row.
289
290 int i,j,rv;
291 int N = TargetNumRows;
292 int mynnz = TargetNumNonzeros;
293 // In the case of reduced communicators, the SourceMatrix won't have the right "MyPID", so thus we have to supply it.
294 int MyPID = MyTargetPID;
295
296 int SizeofIntType = sizeof(int_type);
297
298 // Zero the rowptr
299 for(i=0; i<N+1; i++) CSR_rowptr[i]=0;
300
301 // SameIDs: Always first, always in the same place
302 for(i=0; i<NumSameIDs; i++)
303 CSR_rowptr[i]=SourceMatrix.NumMyEntries(i);
304
305 // PermuteIDs: Still local, but reordered
306 for(i=0; i<NumPermuteIDs; i++)
307 CSR_rowptr[PermuteToLIDs[i]] = SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
308
309 // RemoteIDs: RemoteLIDs tells us the ID, we need to look up the length the hard way. See UnpackAndCombine for where this code came from
310 if(NumRemoteIDs > 0) {
311 double * dintptr = (double *) Imports;
312 int_type * intptr = (int_type *) dintptr;
313 int NumEntries = (int) intptr[1];
314 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
315 for(i=0; i<NumRemoteIDs; i++) {
316 CSR_rowptr[RemoteLIDs[i]] += NumEntries;
317
318 if( i < (NumRemoteIDs-1) ) {
319 dintptr += IntSize + NumEntries;
320 intptr = (int_type *) dintptr;
321 NumEntries = (int) intptr[1];
322 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
323 }
324 }
325 }
326
327 // If multiple procs contribute to a row;
328 std::vector<int> NewStartRow(N+1);
329
330 // Turn row length into a real CSR_rowptr
331 int last_len = CSR_rowptr[0];
332 CSR_rowptr[0] = 0;
333 for(i=1; i<N+1; i++){
334 int new_len = CSR_rowptr[i];
335 CSR_rowptr[i] = last_len + CSR_rowptr[i-1];
336 NewStartRow[i] = CSR_rowptr[i];
337 last_len = new_len;
338 }
339
340 // Preseed TargetPids with -1 for local
341 if(TargetPids.size()!=(size_t)mynnz) TargetPids.resize(mynnz);
342 TargetPids.assign(mynnz,-1);
343
344 // Grab pointers for SourceMatrix
345 int * Source_rowptr, * Source_colind;
346 double * Source_vals;
347 rv=SourceMatrix.ExtractCrsDataPointers(Source_rowptr,Source_colind,Source_vals);
348 if(rv) throw std::runtime_error("UnpackAndCombineIntoCrsArrays: failed in ExtractCrsDataPointers");
349
350 // SameIDs: Copy the data over
351 for(i=0; i<NumSameIDs; i++) {
352 int FromRow = Source_rowptr[i];
353 int ToRow = CSR_rowptr[i];
354 NewStartRow[i] += Source_rowptr[i+1]-Source_rowptr[i];
355
356 for(j=Source_rowptr[i]; j<Source_rowptr[i+1]; j++) {
357 CSR_vals[ToRow + j - FromRow] = Source_vals[j];
358 CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.GCID64(Source_colind[j]);
359 TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
360 }
361 }
362
363 // PermuteIDs: Copy the data over
364 for(i=0; i<NumPermuteIDs; i++) {
365 int FromLID = PermuteFromLIDs[i];
366 int FromRow = Source_rowptr[FromLID];
367 int ToRow = CSR_rowptr[PermuteToLIDs[i]];
368
369 NewStartRow[PermuteToLIDs[i]] += Source_rowptr[FromLID+1]-Source_rowptr[FromLID];
370
371 for(j=Source_rowptr[FromLID]; j<Source_rowptr[FromLID+1]; j++) {
372 CSR_vals[ToRow + j - FromRow] = Source_vals[j];
373 CSR_colind[ToRow + j - FromRow] = (int_type) SourceMatrix.GCID64(Source_colind[j]);
374 TargetPids[ToRow + j - FromRow] = (SourcePids[Source_colind[j]] != MyPID) ? SourcePids[Source_colind[j]] : -1;
375 }
376 }
377
378 // RemoteIDs: Loop structure following UnpackAndCombine
379 if(NumRemoteIDs > 0) {
380 double * dintptr = (double *) Imports;
381 int_type * intptr = (int_type *) dintptr;
382 int NumEntries = (int) intptr[1];
383 int IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
384 double* valptr = dintptr + IntSize;
385
386 for (i=0; i<NumRemoteIDs; i++) {
387 int ToLID = RemoteLIDs[i];
388 int StartRow = NewStartRow[ToLID];
389 NewStartRow[ToLID]+=NumEntries;
390
391 double * values = valptr;
392 int_type * Indices = intptr + 2;
393 for(j=0; j<NumEntries; j++){
394 CSR_vals[StartRow + j] = values[j];
395 CSR_colind[StartRow + j] = Indices[2*j];
396 TargetPids[StartRow + j] = (Indices[2*j+1] != MyPID) ? Indices[2*j+1] : -1;
397 }
398 if( i < (NumRemoteIDs-1) ) {
399 dintptr += IntSize + NumEntries;
400 intptr = (int_type *) dintptr;
401 NumEntries = (int) intptr[1];
402 IntSize = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
403 valptr = dintptr + IntSize;
404 }
405 }
406 }
407
408 return 0;
409}
410
411
412
413// =========================================================================
414// =========================================================================
416 int NumSameIDs,
417 int NumRemoteIDs,
418 const int * RemoteLIDs,
419 int NumPermuteIDs,
420 const int *PermuteToLIDs,
421 const int *PermuteFromLIDs,
422 int LenImports,
423 char* Imports,
424 int TargetNumRows,
425 int TargetNumNonzeros,
426 int MyTargetPID,
427 int * CSR_rowptr,
428 int * CSR_colind,
429 double * CSR_values,
430 const std::vector<int> &SourcePids,
431 std::vector<int> &TargetPids)
432{
433 if(SourceMatrix.RowMap().GlobalIndicesInt()) {
434 return TUnpackAndCombineIntoCrsArrays<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
435 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
436 CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
437 }
438 else
439 throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
440}
441
442// =========================================================================
443// =========================================================================
445 int NumSameIDs,
446 int NumRemoteIDs,
447 const int * RemoteLIDs,
448 int NumPermuteIDs,
449 const int *PermuteToLIDs,
450 const int *PermuteFromLIDs,
451 int LenImports,
452 char* Imports,
453 int TargetNumRows,
454 int TargetNumNonzeros,
455 int MyTargetPID,
456 int * CSR_rowptr,
457 long long * CSR_colind,
458 double * CSR_values,
459 const std::vector<int> &SourcePids,
460 std::vector<int> &TargetPids)
461{
462 if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
463 return TUnpackAndCombineIntoCrsArrays<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
464 PermuteToLIDs,PermuteFromLIDs,LenImports,Imports,TargetNumRows,TargetNumNonzeros,MyTargetPID,
465 CSR_rowptr,CSR_colind,CSR_values,SourcePids,TargetPids);
466 }
467 else
468 throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
469}
470
471
472// =========================================================================
473// =========================================================================
474 template<typename int_type, class MapType1, class MapType2>
475 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)
476 {
477 int i,j;
478
479
480 // Sanity checks
481 bool UseLL;
482 if(domainMap.GlobalIndicesLongLong()) UseLL=true;
483 else if(domainMap.GlobalIndicesInt()) UseLL=false;
484 else throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot detect int type.");
485
486 // Scan all column indices and sort into two groups:
487 // Local: those whose GID matches a GID of the domain map on this processor and
488 // Remote: All others.
489 int numDomainElements = domainMap.NumMyElements();
490 bool * LocalGIDs = 0;
491 if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
492 for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
493
494 bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then error
495 if(DoSizes) throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot handle non-constant sized domainMap.");
496
497
498 // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
499 // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
500 // and the number of block rows.
501 const int numMyBlockRows = N;
502 int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
503 Epetra_HashTable<int_type> RemoteGIDs(hashsize);
504 std::vector<int_type> RemoteGIDList; RemoteGIDList.reserve(hashsize);
505 std::vector<int> PIDList; PIDList.reserve(hashsize);
506
507 // Here we start using the *int* colind array. If int_type==int this clobbers the GIDs, if
508 // int_type==long long, then this is the first use of the colind array.
509 // For *local* GID's set colind with with their LID in the domainMap. For *remote* GIDs,
510 // we set colind with (numDomainElements+NumRemoteColGIDs) before the increment of
511 // the remote count. These numberings will be separate because no local LID is greater
512 // than numDomainElements.
513
514 int NumLocalColGIDs = 0;
515 int NumRemoteColGIDs = 0;
516 for(i = 0; i < numMyBlockRows; i++) {
517 for(j = rowptr[i]; j < rowptr[i+1]; j++) {
518 int_type GID = colind_GID[j];
519 // Check if GID matches a row GID
520 int LID = domainMap.LID(GID);
521 if(LID != -1) {
522 bool alreadyFound = LocalGIDs[LID];
523 if (!alreadyFound) {
524 LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
525 NumLocalColGIDs++;
526 }
527 colind_LID[j] = LID;
528 }
529 else {
530 int_type hash_value=RemoteGIDs.Get(GID);
531 if(hash_value == -1) { // This means its a new remote GID
532 int PID = owningPIDs[j];
533 if(PID==-1) throw std::runtime_error("LowCommunicationMakeColMapAndReindex: Cannot figure out if PID is owned.");
534 colind_LID[j] = numDomainElements + NumRemoteColGIDs;
535 RemoteGIDs.Add(GID, NumRemoteColGIDs);
536 RemoteGIDList.push_back(GID);
537 PIDList.push_back(PID);
538 NumRemoteColGIDs++;
539 }
540 else
541 colind_LID[j] = numDomainElements + hash_value;
542 }
543 }
544 }
545
546 // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
547 if (domainMap.Comm().NumProc()==1) {
548
549 if (NumRemoteColGIDs!=0) {
550 throw std::runtime_error("Some column IDs are not in domainMap. If matrix is rectangular, you must pass in a domainMap");
551 // Sanity test: When one processor,there can be no remoteGIDs
552 }
553 if (NumLocalColGIDs==numDomainElements) {
554 if (LocalGIDs!=0) delete [] LocalGIDs;
555 // In this case, we just use the domainMap's indices, which is, not coincidently, what we clobbered colind with up above anyway.
556 // No further reindexing is needed.
557 NewColMap = domainMap;
558 return 0;
559 }
560 }
561
562 // Now build integer array containing column GIDs
563 // Build back end, containing remote GIDs, first
564 int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
565 std::vector<int_type> ColIndices;
566 int_type * RemoteColIndices=0;
567 if(numMyBlockCols > 0) {
568 ColIndices.resize(numMyBlockCols);
569 if(NumLocalColGIDs!=numMyBlockCols) RemoteColIndices = &ColIndices[NumLocalColGIDs]; // Points to back end of ColIndices
570 else RemoteColIndices=0;
571 }
572
573 for(i = 0; i < NumRemoteColGIDs; i++)
574 RemoteColIndices[i] = RemoteGIDList[i];
575
576 // Build permute array for *remote* reindexing.
577 std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
578 for(i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
579
580 // Sort External column indices so that all columns coming from a given remote processor are contiguous
581 int NumListsInt=0;
582 int NumListsLL =0;
583 int * IntSortLists[2];
584 long long * LLSortLists[2];
585 int * RemotePermuteIDs_ptr = Epetra_Util_data_ptr(RemotePermuteIDs);
586 if(!UseLL) {
587 // int version
588 IntSortLists[0] = (int*) RemoteColIndices;
589 IntSortLists[1] = RemotePermuteIDs_ptr;
590 NumListsInt=2;
591 }
592 else {
593 //LL version
594 LLSortLists[0] = (long long*) RemoteColIndices;
595 IntSortLists[0] = RemotePermuteIDs_ptr;
596 NumListsInt = NumListsLL = 1;
597 }
598
599 int * PIDList_ptr = Epetra_Util_data_ptr(PIDList);
600 Epetra_Util::Sort(true, NumRemoteColGIDs, PIDList_ptr, 0, 0, NumListsInt, IntSortLists,NumListsLL,LLSortLists);
601
602 // Stash the RemotePIDs
603 PIDList.resize(NumRemoteColGIDs);
604 RemotePIDs = PIDList;
605
606 if (SortGhostsAssociatedWithEachProcessor) {
607 // Sort external column indices so that columns from a given remote processor are not only contiguous
608 // but also in ascending order. NOTE: I don't know if the number of externals associated
609 // with a given remote processor is known at this point ... so I count them here.
610
611 // NTS: Only sort the RemoteColIndices this time...
612 int StartCurrent, StartNext;
613 StartCurrent = 0; StartNext = 1;
614 while ( StartNext < NumRemoteColGIDs ) {
615 if (PIDList[StartNext]==PIDList[StartNext-1]) StartNext++;
616 else {
617 IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
618 Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,1,IntSortLists,0,0);
619 StartCurrent = StartNext; StartNext++;
620 }
621 }
622 IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
623 Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, 1,IntSortLists,0,0);
624 }
625
626 // Reverse the permutation to get the information we actually care about
627 std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
628 for(i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
629
630 // Build permute array for *local* reindexing.
631 bool use_local_permute=false;
632 std::vector<int> LocalPermuteIDs(numDomainElements);
633
634 // Now fill front end. Two cases:
635 // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
636 // can simply read the domain GIDs into the front part of ColIndices, otherwise
637 // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
638 // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
639
640 if(NumLocalColGIDs == domainMap.NumMyElements()) {
641 if(NumLocalColGIDs > 0) {
642 domainMap.MyGlobalElements(Epetra_Util_data_ptr(ColIndices)); // Load Global Indices into first numMyBlockCols elements column GID list
643 }
644 }
645 else {
646 int_type* MyGlobalElements = 0;
647 domainMap.MyGlobalElementsPtr(MyGlobalElements);
648
649 int NumLocalAgain = 0;
650 use_local_permute = true;
651 for(i = 0; i < numDomainElements; i++) {
652 if(LocalGIDs[i]) {
653 LocalPermuteIDs[i] = NumLocalAgain;
654 ColIndices[NumLocalAgain++] = MyGlobalElements[i];
655 }
656 }
657 assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
658 }
659
660 // Done with this array
661 if (LocalGIDs!=0) delete [] LocalGIDs;
662
663 // Make Column map with same element sizes as Domain map
664 int_type * ColIndices_ptr = Epetra_Util_data_ptr(ColIndices);
665 MapType2 temp((int_type)(-1), numMyBlockCols, ColIndices_ptr, (int_type)domainMap.IndexBase64(), domainMap.Comm());
666 NewColMap = temp;
667
668 // Low-cost reindex of the matrix
669 for(i=0; i<numMyBlockRows; i++){
670 for(j=rowptr[i]; j<rowptr[i+1]; j++){
671 int ID=colind_LID[j];
672 if(ID < numDomainElements){
673 if(use_local_permute) colind_LID[j] = LocalPermuteIDs[colind_LID[j]];
674 // In the case where use_local_permute==false, we just copy the DomainMap's ordering, which it so happens
675 // is what we put in colind to begin with.
676 }
677 else
678 colind_LID[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_LID[j]-numDomainElements];
679 }
680 }
681
682 return 0;
683}
684
685
686// =========================================================================
687// =========================================================================
688#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
689int LowCommunicationMakeColMapAndReindex(int N, const int * rowptr, int * colind, const Epetra_Map& domainMap, const int * owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, Epetra_BlockMap & NewColMap) {
690
691
692 if(domainMap.GlobalIndicesInt())
693 return TLowCommunicationMakeColMapAndReindex<int,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind,colind,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
694 else
695 throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
696 return -1;
697}
698#endif
699
700// =========================================================================
701// =========================================================================
702#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
703int LowCommunicationMakeColMapAndReindex(int N, const int * rowptr, int * colind_LID, long long * colind_GID, const Epetra_Map& domainMap, const int * owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector<int>& RemotePIDs, Epetra_BlockMap & NewColMap) {
704
705
706 if(domainMap.GlobalIndicesLongLong())
707 return TLowCommunicationMakeColMapAndReindex<long long,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind_LID,colind_GID,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
708 else
709 throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
710 return -1;
711}
712#endif
713
714
715}// end namespace Epetra_Import_Util
#define EPETRA_CHK_ERR(a)
T * Epetra_Util_data_ptr(std::vector< T > &vec)
Function that returns either a pointer to the first entry in the vector or, if the vector is empty,...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
long long IndexBase64() const
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
long long GID64(int LID) const
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool ConstantElementSize() const
Returns true if map has constant element size.
int NumMyElements() const
Number of elements on the calling processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
int MyGlobalElementsPtr(int *&MyGlobalElementList) 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...
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Returns internal data pointers associated with Crs matrix format.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor'...
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
int MaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
long long GCID64(int LCID_in) const
void Add(const long long key, const value_type value)
value_type Get(const long long key)
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
static void EPETRA_LIB_DLL_EXPORT 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)
Epetra_Import_Util: The Epetra ImportUtil Wrapper Namespace.
int TUnpackWithOwningPIDsCount(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *, int NumPermuteIDs, const int *, const int *PermuteFromLIDs, int, char *Imports)
int TUnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int, 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 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 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 TPackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, std::vector< int > &pids)
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.
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 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)