Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
58 namespace Epetra_Import_Util {
59 
60 // =========================================================================
61 // =========================================================================
62 template<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 // =========================================================================
188 template<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 // =========================================================================
261 template<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
689 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) {
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
703 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) {
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
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
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 NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
long long IndexBase64() const
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
value_type Get(const long long key)
bool ConstantElementSize() const
Returns true if map has constant element size.
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...
Definition: Epetra_Util.h:422
#define EPETRA_CHK_ERR(a)
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor&#39;...
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.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumMyElements() const
Number of elements on the calling processor.
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.
long long GID64(int LID) const
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 MaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
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 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)
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int TUnpackWithOwningPIDsCount(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *, int NumPermuteIDs, const int *, const int *PermuteFromLIDs, int, char *Imports)
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 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 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
void Add(const long long key, const value_type value)
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Returns internal data pointers associated with Crs matrix format.