Epetra Package Browser (Single Doxygen Collection)  Development
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 
54 #include <stdexcept>
55 
56 
57 namespace Epetra_Import_Util {
58 
59 // =========================================================================
60 // =========================================================================
61 template<typename int_type>
63  int NumExportIDs,
64  int * ExportLIDs,
65  int & LenExports,
66  char *& Exports,
67  int & SizeOfPacket,
68  int * Sizes,
69  bool & VarSizes,
70  std::vector<int>& pids)
71 {
72  VarSizes = true; //enable variable block size data comm
73 
74  int TotalSendLength = 0;
75  int * IntSizes = 0;
76  if( NumExportIDs>0 ) IntSizes = new int[NumExportIDs];
77 
78  int SizeofIntType = sizeof(int_type);
79 
80  for(int i = 0; i < NumExportIDs; ++i) {
81  int NumEntries;
82  A.NumMyRowEntries( ExportLIDs[i], NumEntries );
83  // Will have NumEntries doubles, 2*NumEntries +2 ints pack them interleaved Sizes[i] = NumEntries;
84  // NTS: We add the owning PID as the SECOND int of the pair for each entry
85  Sizes[i] = NumEntries;
86  // NOTE: Mixing and matching Int Types would be more efficient, BUT what about variable alignment?
87  IntSizes[i] = 1 + (((2*NumEntries+2)*SizeofIntType)/(int)sizeof(double));
88  TotalSendLength += (Sizes[i]+IntSizes[i]);
89  }
90 
91  double * DoubleExports = 0;
92  SizeOfPacket = (int)sizeof(double);
93 
94  //setup buffer locally for memory management by this object
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;
101  }
102 
103  int NumEntries;
104  double * values;
105  double * valptr, * dintptr;
106 
107  // Each segment of Exports will be filled by a packed row of information for each row as follows:
108  // 1st int: GRID of row where GRID is the global row ID for the source matrix
109  // next int: NumEntries, Number of indices in row
110  // next 2*NumEntries: The actual indices and owning [1] PID each for the row in (GID,PID) pairs with the GID first.
111 
112  // [1] Owning is defined in the sense of "Who owns the GID in the DomainMap," aka, who sends the GID in the importer
113 
114  const Epetra_Map & rowMap = A.RowMap();
115  const Epetra_Map & colMap = A.ColMap();
116 
117  if (NumExportIDs > 0) {
118  int_type * Indices;
119  int_type FromRow;
120  int_type * intptr;
121 
122  int maxNumEntries = A.MaxNumEntries();
123  std::vector<int> MyIndices(maxNumEntries);
124  // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
125  // prohibition of unaligned type-punning access.
126  dintptr = (double *) Exports;
127  valptr = dintptr + IntSizes[0];
128  // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
129  // prohibition of unaligned type-punning access.
130  intptr = (int_type *) dintptr;
131  for (int i = 0; i < NumExportIDs; ++i) {
132  FromRow = (int_type) rowMap.GID64(ExportLIDs[i]);
133  intptr[0] = FromRow;
134  values = valptr;
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]); // convert to GIDs
139  Indices[2*j+1] = pids[MyIndices[j]]; // PID owning the entry.
140  }
141  intptr[1] = NumEntries; // Load second slot of segment
142  if (i < (NumExportIDs-1)) {
143  dintptr += (IntSizes[i]+Sizes[i]);
144  valptr = dintptr + IntSizes[i+1];
145  // FIXME (mfh 11 Feb 2015) This probably violates ANSI C++'s
146  // prohibition of unaligned type-punning access.
147  intptr = (int_type *) dintptr;
148  }
149  }
150 
151  for (int i = 0; i < NumExportIDs; ++i) {
152  Sizes[i] += IntSizes[i];
153  }
154  }
155 
156  if( IntSizes ) delete [] IntSizes;
157 
158  return 0;
159 }
160 
161 
162 // =========================================================================
163 // =========================================================================
165  int NumExportIDs,
166  int * ExportLIDs,
167  int & LenExports,
168  char *& Exports,
169  int & SizeOfPacket,
170  int * Sizes,
171  bool & VarSizes,
172  std::vector<int>& SourcePids)
173 {
174  if(SourceMatrix.RowMap().GlobalIndicesInt()) {
175  return TPackAndPrepareWithOwningPIDs<int>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
176  }
177  else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
178  return TPackAndPrepareWithOwningPIDs<long long>(SourceMatrix,NumExportIDs,ExportLIDs,LenExports,Exports,SizeOfPacket,Sizes,VarSizes,SourcePids);
179  }
180  else
181  throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
182 }
183 
184 
185 // =========================================================================
186 // =========================================================================
187 template<typename int_type>
189  int NumSameIDs,
190  int NumRemoteIDs,
191  const int * RemoteLIDs,
192  int NumPermuteIDs,
193  const int *PermuteToLIDs,
194  const int *PermuteFromLIDs,
195  int LenImports,
196  char* Imports)
197 {
198  int i,nnz=0;
199  int SizeofIntType = (int)sizeof(int_type);
200 
201  // SameIDs: Always first, always in the same place
202  for(i=0; i<NumSameIDs; i++)
203  nnz+=SourceMatrix.NumMyEntries(i);
204 
205  // PermuteIDs: Still local, but reordered
206  for(i=0; i<NumPermuteIDs; i++)
207  nnz += SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
208 
209  // RemoteIDs: RemoteLIDs tells us the ID, we need to look up the length the hard way. See UnpackAndCombine for where this code came from
210  if(NumRemoteIDs > 0) {
211  double * dintptr = (double *) Imports;
212  // General version
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++) {
217  nnz += NumEntries;
218 
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));
224  }
225  }
226  }
227 
228  return nnz;
229 }
230 
231 
232 // =========================================================================
233 // =========================================================================
235  int NumSameIDs,
236  int NumRemoteIDs,
237  const int * RemoteLIDs,
238  int NumPermuteIDs,
239  const int *PermuteToLIDs,
240  const int *PermuteFromLIDs,
241  int LenImports,
242  char* Imports)
243 {
244  if(SourceMatrix.RowMap().GlobalIndicesInt()) {
245  return TUnpackWithOwningPIDsCount<int>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
246  PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
247  }
248  else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
249  return TUnpackWithOwningPIDsCount<long long>(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,
250  PermuteToLIDs,PermuteFromLIDs,LenImports,Imports);
251  }
252  else
253  throw std::runtime_error("UnpackWithOwningPIDsCount: Unable to determine source global index type");
254 
255 }
256 
257 
258 // =========================================================================
259 // =========================================================================
260 template<typename int_type>
262  int NumSameIDs,
263  int NumRemoteIDs,
264  const int * RemoteLIDs,
265  int NumPermuteIDs,
266  const int *PermuteToLIDs,
267  const int *PermuteFromLIDs,
268  int LenImports,
269  char* Imports,
270  int TargetNumRows,
271  int TargetNumNonzeros,
272  int MyTargetPID,
273  int * CSR_rowptr,
274  int_type * CSR_colind,
275  double * CSR_vals,
276  const std::vector<int> &SourcePids,
277  std::vector<int> &TargetPids)
278 {
279  // What we really need to know is where in the CSR arrays each row should start (aka the rowptr).
280  // 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
281  // (c) Having each row copied into the right colind / values locations.
282 
283  // From Epetra_CrsMatrix UnpackAndCombine():
284  // Each segment of Exports will be filled by a packed row of information for each row as follows:
285  // 1st int: GRID of row where GRID is the global row ID for the source matrix
286  // next int: NumEntries, Number of indices in row.
287  // next NumEntries: The actual indices for the row.
288 
289  int i,j,rv;
290  int N = TargetNumRows;
291  int mynnz = TargetNumNonzeros;
292  // In the case of reduced communicators, the SourceMatrix won't have the right "MyPID", so thus we have to supply it.
293  int MyPID = MyTargetPID;
294 
295  int SizeofIntType = sizeof(int_type);
296 
297  // Zero the rowptr
298  for(i=0; i<N+1; i++) CSR_rowptr[i]=0;
299 
300  // SameIDs: Always first, always in the same place
301  for(i=0; i<NumSameIDs; i++)
302  CSR_rowptr[i]=SourceMatrix.NumMyEntries(i);
303 
304  // PermuteIDs: Still local, but reordered
305  for(i=0; i<NumPermuteIDs; i++)
306  CSR_rowptr[PermuteToLIDs[i]] = SourceMatrix.NumMyEntries(PermuteFromLIDs[i]);
307 
308  // RemoteIDs: RemoteLIDs tells us the ID, we need to look up the length the hard way. See UnpackAndCombine for where this code came from
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;
316 
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));
322  }
323  }
324  }
325 
326  // If multiple procs contribute to a row;
327  std::vector<int> NewStartRow(N+1);
328 
329  // Turn row length into a real CSR_rowptr
330  int last_len = CSR_rowptr[0];
331  CSR_rowptr[0] = 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];
336  last_len = new_len;
337  }
338 
339  // Preseed TargetPids with -1 for local
340  if(TargetPids.size()!=(size_t)mynnz) TargetPids.resize(mynnz);
341  TargetPids.assign(mynnz,-1);
342 
343  // Grab pointers for SourceMatrix
344  int * Source_rowptr, * Source_colind;
345  double * Source_vals;
346  rv=SourceMatrix.ExtractCrsDataPointers(Source_rowptr,Source_colind,Source_vals);
347  if(rv) throw std::runtime_error("UnpackAndCombineIntoCrsArrays: failed in ExtractCrsDataPointers");
348 
349  // SameIDs: Copy the data over
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];
354 
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;
359  }
360  }
361 
362  // PermuteIDs: Copy the data over
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]];
367 
368  NewStartRow[PermuteToLIDs[i]] += Source_rowptr[FromLID+1]-Source_rowptr[FromLID];
369 
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;
374  }
375  }
376 
377  // RemoteIDs: Loop structure following UnpackAndCombine
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;
384 
385  for (i=0; i<NumRemoteIDs; i++) {
386  int ToLID = RemoteLIDs[i];
387  int StartRow = NewStartRow[ToLID];
388  NewStartRow[ToLID]+=NumEntries;
389 
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;
396  }
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;
403  }
404  }
405  }
406 
407  return 0;
408 }
409 
410 
411 
412 // =========================================================================
413 // =========================================================================
415  int NumSameIDs,
416  int NumRemoteIDs,
417  const int * RemoteLIDs,
418  int NumPermuteIDs,
419  const int *PermuteToLIDs,
420  const int *PermuteFromLIDs,
421  int LenImports,
422  char* Imports,
423  int TargetNumRows,
424  int TargetNumNonzeros,
425  int MyTargetPID,
426  int * CSR_rowptr,
427  int * CSR_colind,
428  double * CSR_values,
429  const std::vector<int> &SourcePids,
430  std::vector<int> &TargetPids)
431 {
432  if(SourceMatrix.RowMap().GlobalIndicesInt()) {
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);
436  }
437  else
438  throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
439 }
440 
441 // =========================================================================
442 // =========================================================================
444  int NumSameIDs,
445  int NumRemoteIDs,
446  const int * RemoteLIDs,
447  int NumPermuteIDs,
448  const int *PermuteToLIDs,
449  const int *PermuteFromLIDs,
450  int LenImports,
451  char* Imports,
452  int TargetNumRows,
453  int TargetNumNonzeros,
454  int MyTargetPID,
455  int * CSR_rowptr,
456  long long * CSR_colind,
457  double * CSR_values,
458  const std::vector<int> &SourcePids,
459  std::vector<int> &TargetPids)
460 {
461  if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
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);
465  }
466  else
467  throw std::runtime_error("UnpackAndCombineIntoCrsArrays: int type not matched.");
468 }
469 
470 
471 // =========================================================================
472 // =========================================================================
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)
475  {
476  int i,j;
477 
478 
479  // Sanity checks
480  bool UseLL;
481  if(domainMap.GlobalIndicesLongLong()) UseLL=true;
482  else if(domainMap.GlobalIndicesInt()) UseLL=false;
483  else throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot detect int type.");
484 
485  // Scan all column indices and sort into two groups:
486  // Local: those whose GID matches a GID of the domain map on this processor and
487  // Remote: All others.
488  int numDomainElements = domainMap.NumMyElements();
489  bool * LocalGIDs = 0;
490  if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
491  for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
492 
493  bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then error
494  if(DoSizes) throw std::runtime_error("LowCommunicationMakeColMapAndReindex: cannot handle non-constant sized domainMap.");
495 
496 
497  // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
498  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
499  // and the number of block rows.
500  const int numMyBlockRows = N;
501  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
502  Epetra_HashTable<int_type> RemoteGIDs(hashsize);
503  std::vector<int_type> RemoteGIDList; RemoteGIDList.reserve(hashsize);
504  std::vector<int> PIDList; PIDList.reserve(hashsize);
505 
506  // Here we start using the *int* colind array. If int_type==int this clobbers the GIDs, if
507  // int_type==long long, then this is the first use of the colind array.
508  // For *local* GID's set colind with with their LID in the domainMap. For *remote* GIDs,
509  // we set colind with (numDomainElements+NumRemoteColGIDs) before the increment of
510  // the remote count. These numberings will be separate because no local LID is greater
511  // than numDomainElements.
512 
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];
518  // Check if GID matches a row GID
519  int LID = domainMap.LID(GID);
520  if(LID != -1) {
521  bool alreadyFound = LocalGIDs[LID];
522  if (!alreadyFound) {
523  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
524  NumLocalColGIDs++;
525  }
526  colind_LID[j] = LID;
527  }
528  else {
529  int_type hash_value=RemoteGIDs.Get(GID);
530  if(hash_value == -1) { // This means its a new remote GID
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);
537  NumRemoteColGIDs++;
538  }
539  else
540  colind_LID[j] = numDomainElements + hash_value;
541  }
542  }
543  }
544 
545  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
546  if (domainMap.Comm().NumProc()==1) {
547 
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");
550  // Sanity test: When one processor,there can be no remoteGIDs
551  }
552  if (NumLocalColGIDs==numDomainElements) {
553  if (LocalGIDs!=0) delete [] LocalGIDs;
554  // In this case, we just use the domainMap's indices, which is, not coincidently, what we clobbered colind with up above anyway.
555  // No further reindexing is needed.
556  NewColMap = domainMap;
557  return 0;
558  }
559  }
560 
561  // Now build integer array containing column GIDs
562  // Build back end, containing remote GIDs, first
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]; // Points to back end of ColIndices
569  else RemoteColIndices=0;
570  }
571 
572  for(i = 0; i < NumRemoteColGIDs; i++)
573  RemoteColIndices[i] = RemoteGIDList[i];
574 
575  // Build permute array for *remote* reindexing.
576  std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
577  for(i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
578 
579  // Sort External column indices so that all columns coming from a given remote processor are contiguous
580  int NumListsInt=0;
581  int NumListsLL =0;
582  int * IntSortLists[2];
583  long long * LLSortLists[2];
584  int * RemotePermuteIDs_ptr = RemotePermuteIDs.size() ? &RemotePermuteIDs[0] : 0;
585  if(!UseLL) {
586  // int version
587  IntSortLists[0] = (int*) RemoteColIndices;
588  IntSortLists[1] = RemotePermuteIDs_ptr;
589  NumListsInt=2;
590  }
591  else {
592  //LL version
593  LLSortLists[0] = (long long*) RemoteColIndices;
594  IntSortLists[0] = RemotePermuteIDs_ptr;
595  NumListsInt = NumListsLL = 1;
596  }
597 
598  int * PIDList_ptr = PIDList.size() ? &PIDList[0] : 0;
599  Epetra_Util::Sort(true, NumRemoteColGIDs, PIDList_ptr, 0, 0, NumListsInt, IntSortLists,NumListsLL,LLSortLists);
600 
601  // Stash the RemotePIDs
602  PIDList.resize(NumRemoteColGIDs);
603  RemotePIDs = PIDList;
604 
605  if (SortGhostsAssociatedWithEachProcessor) {
606  // Sort external column indices so that columns from a given remote processor are not only contiguous
607  // but also in ascending order. NOTE: I don't know if the number of externals associated
608  // with a given remote processor is known at this point ... so I count them here.
609 
610  // NTS: Only sort the RemoteColIndices this time...
611  int StartCurrent, StartNext;
612  StartCurrent = 0; StartNext = 1;
613  while ( StartNext < NumRemoteColGIDs ) {
614  if (PIDList[StartNext]==PIDList[StartNext-1]) StartNext++;
615  else {
616  IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
617  Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,1,IntSortLists,0,0);
618  StartCurrent = StartNext; StartNext++;
619  }
620  }
621  IntSortLists[0] = &RemotePermuteIDs[StartCurrent];
622  Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, 1,IntSortLists,0,0);
623  }
624 
625  // Reverse the permutation to get the information we actually care about
626  std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
627  for(i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
628 
629  // Build permute array for *local* reindexing.
630  bool use_local_permute=false;
631  std::vector<int> LocalPermuteIDs(numDomainElements);
632 
633  // Now fill front end. Two cases:
634  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
635  // can simply read the domain GIDs into the front part of ColIndices, otherwise
636  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
637  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
638 
639  if(NumLocalColGIDs == domainMap.NumMyElements()) {
640  if(NumLocalColGIDs > 0) {
641  domainMap.MyGlobalElements(ColIndices.size() ? &ColIndices[0] : 0); // Load Global Indices into first numMyBlockCols elements column GID list
642  }
643  }
644  else {
645  int_type* MyGlobalElements = 0;
646  domainMap.MyGlobalElementsPtr(MyGlobalElements);
647 
648  int NumLocalAgain = 0;
649  use_local_permute = true;
650  for(i = 0; i < numDomainElements; i++) {
651  if(LocalGIDs[i]) {
652  LocalPermuteIDs[i] = NumLocalAgain;
653  ColIndices[NumLocalAgain++] = MyGlobalElements[i];
654  }
655  }
656  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
657  }
658 
659  // Done with this array
660  if (LocalGIDs!=0) delete [] LocalGIDs;
661 
662  // Make Column map with same element sizes as Domain map
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());
665  NewColMap = temp;
666 
667  // Low-cost reindex of the matrix
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]];
673  // In the case where use_local_permute==false, we just copy the DomainMap's ordering, which it so happens
674  // is what we put in colind to begin with.
675  }
676  else
677  colind_LID[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_LID[j]-numDomainElements];
678  }
679  }
680 
681  return 0;
682 }
683 
684 
685 // =========================================================================
686 // =========================================================================
687 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
688 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) {
689 
690 
691  if(domainMap.GlobalIndicesInt())
692  return TLowCommunicationMakeColMapAndReindex<int,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind,colind,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
693  else
694  throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
695  return -1;
696 }
697 #endif
698 
699 // =========================================================================
700 // =========================================================================
701 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
702 int 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) {
703 
704 
705  if(domainMap.GlobalIndicesLongLong())
706  return TLowCommunicationMakeColMapAndReindex<long long,Epetra_BlockMap,Epetra_Map>(N,rowptr,colind_LID,colind_GID,domainMap,owningPIDs,SortGhostsAssociatedWithEachProcessor,RemotePIDs,NewColMap);
707  else
708  throw std::runtime_error("LowCommunicationMakeColMapAndReindex: GID type mismatch.");
709  return -1;
710 }
711 #endif
712 
713 
714 }// end namespace Epetra_Import_Util
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int MaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on this 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 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.
Definition: Epetra_Map.h:119
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 NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
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.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor&#39;...
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)