EpetraExt  Development
EpetraExt_MMHelpers.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
42 #include <EpetraExt_ConfigDefs.h>
43 #include <EpetraExt_MMHelpers.h>
44 #include <Epetra_Comm.h>
45 #include <Epetra_CrsMatrix.h>
46 #include <Epetra_Import.h>
47 #include <Epetra_Distributor.h>
48 #include <Epetra_HashTable.h>
49 #include <Epetra_Util.h>
50 #include <Epetra_Import_Util.h>
51 #include <Epetra_GIDTypeSerialDenseVector.h>
52 
53 #include <Teuchos_TimeMonitor.hpp>
54 #include <limits>
55 
56 #ifdef HAVE_MPI
57 #include "Epetra_MpiComm.h"
58 #include "Epetra_MpiDistributor.h"
59 #endif
60 #define MIN(x,y) ((x)<(y)?(x):(y))
61 #define MIN3(x,y,z) ((x)<(y)?(MIN(x,z)):(MIN(y,z)))
62 
63 namespace EpetraExt {
64 
65 //------------------------------------
66 // DEBUGGING ROUTINES
67 void debug_print_distor(const char * label, const Epetra_Distributor * Distor, const Epetra_Comm & Comm) {
68 #ifdef HAVE_MPI
69  const Epetra_MpiDistributor * MDistor = dynamic_cast<const Epetra_MpiDistributor*>(Distor);
70  printf("[%d] %s\n",Comm.MyPID(),label);
71  printf("[%d] NumSends = %d NumRecvs = %d\n",Comm.MyPID(),MDistor->NumSends(),MDistor->NumReceives());
72  printf("[%d] ProcsTo = ",Comm.MyPID());
73  for(int ii=0; ii<MDistor->NumSends(); ii++)
74  printf("%d ",MDistor->ProcsTo()[ii]);
75  printf("\n[%d] ProcsFrom = ",Comm.MyPID());
76  for(int ii=0; ii<MDistor->NumReceives(); ii++)
77  printf("%d ",MDistor->ProcsFrom()[ii]);
78  printf("\n");
79  fflush(stdout);
80 #endif
81 }
82 
83 //------------------------------------
84 // DEBUGGING ROUTINES
85 void debug_compare_import(const Epetra_Import * Import1,const Epetra_Import * Import2) {
86  if(!Import1 && !Import2) return;
87  const Epetra_Comm & Comm = (Import1)? Import1->SourceMap().Comm() : Import2->SourceMap().Comm();
88  bool flag=true;
89  int PID=Comm.MyPID();
90  if( (!Import1 && Import2) || (Import2 && !Import1) ) {printf("[%d] DCI: One Import exists, the other does not\n",PID);return;}
91  if(!Import1->SourceMap().SameAs(Import2->SourceMap())) {printf("[%d] DCI: SourceMaps don't match\n",PID);return;}
92  if(!Import1->TargetMap().SameAs(Import2->TargetMap())) {printf("[%d] DCI: TargetMaps don't match\n",PID);return;}
93 
94  if(Import1->NumSameIDs() != Import2->NumSameIDs()) {printf("[%d] DCI NumSameIDs() mismatch %d vs. %d\n",PID,Import1->NumSameIDs(),Import2->NumSameIDs());flag=false;}
95 
96  if(Import1->NumPermuteIDs() != Import2->NumPermuteIDs()) {printf("[%d] DCI NumPermuteIDs() mismatch %d vs. %d\n",PID,Import1->NumPermuteIDs(),Import2->NumPermuteIDs()); flag=false;}
97 
98  if(Import1->NumExportIDs() != Import2->NumExportIDs()) {printf("[%d] DCI NumExportIDs() mismatch %d vs. %d\n",PID,Import1->NumExportIDs(),Import2->NumExportIDs()); flag=false;}
99 
100  if(Import1->NumRemoteIDs() != Import2->NumRemoteIDs()) {printf("[%d] DCI NumRemoteIDs() mismatch %d vs. %d\n",PID,Import1->NumRemoteIDs(),Import2->NumRemoteIDs()); flag=false;}
101 
102  if(Import1->NumSend() != Import2->NumSend()) {printf("[%d] DCI NumSend() mismatch %d vs. %d\n",PID,Import1->NumSend(),Import2->NumSend()); flag=false;}
103 
104  if(Import1->NumRecv() != Import2->NumRecv()) {printf("[%d] DCI NumRecv() mismatch %d vs. %d\n",PID,Import1->NumRecv(),Import2->NumRecv()); flag=false;}
105 
106 
107  if(flag) printf("[%d] DCI Importers compare OK\n",PID);
108  fflush(stdout);
109  Import1->SourceMap().Comm().Barrier();
110  Import1->SourceMap().Comm().Barrier();
111  Import1->SourceMap().Comm().Barrier();
112  if(!flag) exit(1);
113 }
114 
115 
116 
117 //------------------------------------
119  : numRows(0), numEntriesPerRow(NULL), indices(NULL), values(NULL),
120  remote(NULL), numRemote(0), importColMap(NULL), rowMap(NULL), colMap(NULL),
121  domainMap(NULL), importMatrix(NULL), origMatrix(NULL)
122 {
123 }
124 
126 {
127  deleteContents();
128 }
129 
131 {
132  numRows = 0;
133  delete [] numEntriesPerRow; numEntriesPerRow = NULL;
134  delete [] indices; indices = NULL;
135  delete [] values; values = NULL;
136  delete [] remote; remote = NULL;
137  delete importColMap; importColMap = NULL;
138  numRemote = 0;
139  delete importMatrix; importMatrix=0;
140  // origMatrix is not owned by me, so don't delete
141  origMatrix=0;
142  targetMapToOrigRow.resize(0);
143  targetMapToImportRow.resize(0);
144 }
145 
147 {
148  std::cout << "proc " << M.rowMap->Comm().MyPID()<<std::endl;
149  std::cout << "numRows: " << M.numRows<<std::endl;
150  for(int i=0; i<M.numRows; ++i) {
151  for(int j=0; j<M.numEntriesPerRow[i]; ++j) {
152  if (M.remote[i]) {
153  std::cout << " *"<<M.rowMap->GID64(i)<<" "
154  <<M.importColMap->GID64(M.indices[i][j])<<" "<<M.values[i][j]<<std::endl;
155  }
156  else {
157  std::cout << " "<<M.rowMap->GID64(i)<<" "
158  <<M.colMap->GID64(M.indices[i][j])<<" "<<M.values[i][j]<<std::endl;
159  }
160  }
161  }
162  return(0);
163 }
164 
166  : ecrsmat_(epetracrsmatrix)
167 {
168 }
169 
171 {
172 }
173 
174 const Epetra_Map&
176 {
177  return ecrsmat_.RowMap();
178 }
179 
181 {
182  return ecrsmat_.Filled();
183 }
184 
185 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
186 int
187 CrsWrapper_Epetra_CrsMatrix::InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
188 {
189  return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
190 }
191 
192 int
193 CrsWrapper_Epetra_CrsMatrix::SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices)
194 {
195  return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
196 }
197 #endif
198 
199 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
200 int
201 CrsWrapper_Epetra_CrsMatrix::InsertGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices)
202 {
203  return ecrsmat_.InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
204 }
205 
206 int
207 CrsWrapper_Epetra_CrsMatrix::SumIntoGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices)
208 {
209  return ecrsmat_.SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices);
210 }
211 #endif
212 
213 
214 //------------------------------------
215 
216 template<typename int_type>
218  : graph_(),
219  rowmap_(emap),
220  max_row_length_(0)
221 {
222  int num_rows = emap.NumMyElements();
223  int_type* rows = 0;
224  emap.MyGlobalElementsPtr(rows);
225 
226  for(int i=0; i<num_rows; ++i) {
227  graph_[rows[i]] = new std::set<int_type>;
228  }
229 }
230 
231 template<typename int_type>
233 {
234  typename std::map<int_type,std::set<int_type>*>::iterator
235  iter = graph_.begin(), iter_end = graph_.end();
236  for(; iter!=iter_end; ++iter) {
237  delete iter->second;
238  }
239 
240  graph_.clear();
241 }
242 
243 template<typename int_type>
245 {
246  return false;
247 }
248 
249 template<typename int_type>
250 int
251 CrsWrapper_GraphBuilder<int_type>::InsertGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices)
252 {
253  typename std::map<int_type,std::set<int_type>*>::iterator
254  iter = graph_.find(GlobalRow);
255 
256  if (iter == graph_.end()) return(-1);
257 
258  std::set<int_type>& cols = *(iter->second);
259 
260  for(int i=0; i<NumEntries; ++i) {
261  cols.insert(Indices[i]);
262  }
263 
264  int row_length = cols.size();
265  if (row_length > max_row_length_) max_row_length_ = row_length;
266 
267  return(0);
268 }
269 
270 template<typename int_type>
271 int
272 CrsWrapper_GraphBuilder<int_type>::SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices)
273 {
274  return InsertGlobalValues(GlobalRow, NumEntries, Values, Indices);
275 }
276 
277 template<typename int_type>
278 std::map<int_type,std::set<int_type>*>&
280 {
281  return graph_;
282 }
283 
284 template<typename int_type>
286  Epetra_CrsMatrix& C)
287 {
288  int max_row_length = graphbuilder.get_max_row_length();
289  if (max_row_length < 1) return;
290 
291  std::vector<int_type> indices(max_row_length);
292  int_type* indices_ptr = &indices[0];
293  std::vector<double> zeros(max_row_length, 0.0);
294  double* zeros_ptr = &zeros[0];
295 
296  std::map<int_type,std::set<int_type>*>& graph = graphbuilder.get_graph();
297 
298  typename std::map<int_type,std::set<int_type>*>::iterator
299  iter = graph.begin(), iter_end = graph.end();
300 
301  for(; iter!=iter_end; ++iter) {
302  int_type row = iter->first;
303  std::set<int_type>& cols = *(iter->second);
304  int num_entries = cols.size();
305 
306  typename std::set<int_type>::iterator
307  col_iter = cols.begin(), col_end = cols.end();
308  for(int j=0; col_iter!=col_end; ++col_iter, ++j) {
309  indices_ptr[j] = *col_iter;
310  }
311 
312  C.InsertGlobalValues(row, num_entries, zeros_ptr, indices_ptr);
313  }
314 }
315 
316 template<typename int_type>
317 void Tpack_outgoing_rows(const Epetra_CrsMatrix& mtx,
318  const std::vector<int_type>& proc_col_ranges,
319  std::vector<int_type>& send_rows,
320  std::vector<int>& rows_per_send_proc)
321 {
322  const Epetra_Map& rowmap = mtx.RowMap();
323  int numrows = mtx.NumMyRows();
324  const Epetra_CrsGraph& graph = mtx.Graph();
325  int rowlen = 0;
326  int* col_indices = NULL;
327  int_type* Tcol_indices = NULL;
328  int num_col_ranges = proc_col_ranges.size()/2;
329  rows_per_send_proc.resize(num_col_ranges);
330  send_rows.clear();
331  for(int nc=0; nc<num_col_ranges; ++nc) {
332  int_type first_col = proc_col_ranges[nc*2];
333  int_type last_col = proc_col_ranges[nc*2+1];
334  int num_send_rows = 0;
335  for(int i=0; i<numrows; ++i) {
336  int_type grow = (int_type) rowmap.GID64(i);
337  if (mtx.Filled()) {
338  const Epetra_Map& colmap = mtx.ColMap();
339  graph.ExtractMyRowView(i, rowlen, col_indices);
340  for(int j=0; j<rowlen; ++j) {
341  int_type col = (int_type) colmap.GID64(col_indices[j]);
342  if (first_col <= col && last_col >= col) {
343  ++num_send_rows;
344  send_rows.push_back(grow);
345  break;
346  }
347  }
348  }
349  else {
350  graph.ExtractGlobalRowView(grow, rowlen, Tcol_indices);
351  for(int j=0; j<rowlen; ++j) {
352  if (first_col <= Tcol_indices[j] && last_col >= Tcol_indices[j]) {
353  ++num_send_rows;
354  send_rows.push_back(grow);
355  break;
356  }
357  }
358  }
359  }
360  rows_per_send_proc[nc] = num_send_rows;
361  }
362 }
363 
364 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
365 void pack_outgoing_rows(const Epetra_CrsMatrix& mtx,
366  const std::vector<int>& proc_col_ranges,
367  std::vector<int>& send_rows,
368  std::vector<int>& rows_per_send_proc)
369 {
370  if(mtx.RowMatrixRowMap().GlobalIndicesInt()) {
371  Tpack_outgoing_rows<int>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
372  }
373  else {
374  throw "EpetraExt::pack_outgoing_rows: Global indices not int";
375  }
376 }
377 #endif
378 
379 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
380 void pack_outgoing_rows(const Epetra_CrsMatrix& mtx,
381  const std::vector<long long>& proc_col_ranges,
382  std::vector<long long>& send_rows,
383  std::vector<int>& rows_per_send_proc)
384 {
385  if(mtx.RowMatrixRowMap().GlobalIndicesLongLong()) {
386  Tpack_outgoing_rows<long long>(mtx, proc_col_ranges, send_rows, rows_per_send_proc);
387  }
388  else {
389  throw "EpetraExt::pack_outgoing_rows: Global indices not long long";
390  }
391 }
392 #endif
393 
394 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
395 template<>
396 std::pair<int,int> get_col_range<int>(const Epetra_Map& emap)
397  {
398  if(emap.GlobalIndicesInt())
399  return std::make_pair(emap.MinMyGID(),emap.MaxMyGID());
400  else
401  throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_Map";
402 }
403 #endif
404 
405 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
406 template<>
407 std::pair<long long,long long> get_col_range<long long>(const Epetra_Map& emap)
408 {
409  if(emap.GlobalIndicesLongLong())
410  return std::make_pair(emap.MinMyGID64(),emap.MaxMyGID64());
411  else
412  throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_Map";
413 }
414 #endif
415 
416 template<typename int_type>
417 std::pair<int_type,int_type> Tget_col_range(const Epetra_CrsMatrix& mtx)
418 {
419  std::pair<int_type,int_type> col_range;
420  if (mtx.Filled()) {
421  col_range = get_col_range<int_type>(mtx.ColMap());
422  }
423  else {
424  const Epetra_Map& row_map = mtx.RowMap();
425  col_range.first = row_map.MaxMyGID64();
426  col_range.second = row_map.MinMyGID64();
427  int rowlen = 0;
428  int_type* col_indices = NULL;
429  const Epetra_CrsGraph& graph = mtx.Graph();
430  for(int i=0; i<row_map.NumMyElements(); ++i) {
431  graph.ExtractGlobalRowView((int_type) row_map.GID64(i), rowlen, col_indices);
432  for(int j=0; j<rowlen; ++j) {
433  if (col_indices[j] < col_range.first) col_range.first = col_indices[j];
434  if (col_indices[j] > col_range.second) col_range.second = col_indices[j];
435  }
436  }
437  }
438 
439  return col_range;
440 }
441 
442 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
443 template<>
444 std::pair<int,int> get_col_range<int>(const Epetra_CrsMatrix& mtx)
445 {
446  if(mtx.RowMatrixRowMap().GlobalIndicesInt())
447  return Tget_col_range<int>(mtx);
448  else
449  throw "EpetraExt::get_col_range<int>: Unknown Global Indices for Epetra_CrsMatrix";
450 }
451 #endif
452 
453 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
454 template<>
455 std::pair<long long,long long> get_col_range<long long>(const Epetra_CrsMatrix& mtx)
456 {
457  if(mtx.RowMatrixRowMap().GlobalIndicesLongLong())
458  return Tget_col_range<long long>(mtx);
459  else
460  throw "EpetraExt::get_col_range<long long>: Unknown Global Indices for Epetra_CrsMatrix";
461 }
462 #endif
463 
464 
465 /**********************************************************************************************/
466 /**********************************************************************************************/
467 /**********************************************************************************************/
468 #ifdef HAVE_MPI
469 template <typename MyType>
470 void boundary_exchange(const Epetra_MpiComm Comm, MPI_Datatype DataType,
471  int NumSends, const int * SendProcs, const int * SendSizes, MyType* SendBuffer,
472  int NumRecvs, const int * RecvProcs, const int * RecvSizes, MyType* RecvBuffer,int SizeOfPacket,int msg_tag)
473 {
474 
475  MPI_Comm comm = Comm.Comm();
476  std::vector<MPI_Request> requests(NumRecvs);
477  std::vector<MPI_Status> status(NumRecvs);
478 
479  int i,num_waits=0,MyPID=Comm.MyPID();
480  int start, self_recv_len=-1,self_recv_start=-1, self_send_start=-1;
481 
482  // Default send/recv size if the Sizes arrays are NULL.
483  int mysendsize=1, myrecvsize=1;
484 
485  // Post Recvs
486  start=0;
487  for(i=0; i<NumRecvs; i++){
488  if(RecvSizes) myrecvsize=RecvSizes[i]*SizeOfPacket;
489  if(RecvProcs[i] != MyPID) {
490  MPI_Irecv(RecvBuffer + start, myrecvsize, DataType, RecvProcs[i], msg_tag, comm, &requests[num_waits]);
491  num_waits++;
492  }
493  else {
494  self_recv_len = myrecvsize;
495  self_recv_start=start;
496  }
497  start+=myrecvsize;
498  }
499 
500  // Do sends
501  start=0;
502  for(i=0; i<NumSends; i++){
503  if(SendSizes) mysendsize=SendSizes[i]*SizeOfPacket;
504  if(SendProcs[i] != MyPID)
505  MPI_Send(SendBuffer + start, mysendsize,DataType,SendProcs[i],msg_tag,comm);
506  else
507  self_send_start=start;
508  start+=mysendsize;
509  }
510 
511  // Self-copy (if needed)
512  if(self_recv_len != -1)
513  memcpy(RecvBuffer+self_recv_start,SendBuffer+self_send_start,self_recv_len*sizeof(MyType)*SizeOfPacket);
514 
515  // Wait
516  if(NumRecvs > 0)
517  MPI_Waitall(num_waits, &requests[0],&status[0]);
518 }
519 #endif
520 
521 
522 #ifdef HAVE_MPI
523 template <typename MyType>
524 void boundary_exchange_varsize(const Epetra_MpiComm Comm, MPI_Datatype DataType,
525  int NumSends, const int * SendProcs, const int * SendSizes, MyType* SendBuffer,
526  int NumRecvs, const int * RecvProcs, int * RecvSizes, MyType*& RecvBuffer,int SizeOfPacket,int msg_tag)
527 {
528 
529  int i,rbuffersize=0;
530 
531  // Do a first round of boundary exchange with the the SendBuffer sizes
532  boundary_exchange<int>(Comm,MPI_INT,NumSends,SendProcs,(int*)0,const_cast<int*>(SendSizes),NumRecvs,RecvProcs,(int*)0,RecvSizes,1,msg_tag);
533 
534  // Allocate the RecvBuffer
535  for(i=0; i<NumRecvs; i++) rbuffersize+=RecvSizes[i]*SizeOfPacket;
536  RecvBuffer = new MyType[rbuffersize];
537 
538  // Do a second round of boundary exchange to trade the actual values
539  boundary_exchange<MyType>(Comm,DataType,NumSends,SendProcs,SendSizes,SendBuffer,NumRecvs,RecvProcs,RecvSizes,RecvBuffer,SizeOfPacket,msg_tag+100);
540 }
541 #endif
542 
543 
544 //=========================================================================
545 //=========================================================================
546 //=========================================================================
548  Epetra_Data(),
549  IndexBase_(0),
550 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
551  LIDHash_int_(0),
552 #endif
553 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
554  LIDHash_LL_(0),
555 #endif
556  CopyMap_(0)
557 {
558 }
559 //=========================================================================
561 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
562  delete LIDHash_int_;
563 #endif
564 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
565  delete LIDHash_LL_;
566 #endif
567  delete CopyMap_;
568 }
569 
570 //=========================================================================
571 LightweightMap::LightweightMap():Data_(0),IsLongLong(false),IsInt(false){;}
572 
573 //=========================================================================
574 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
575 void LightweightMap::Construct_int(int /* numGlobalElements */, int numMyElements, const int * myGlobalElements, long long /* indexBase */, bool GenerateHash)
576 {
577  Data_=new LightweightMapData();
578  Data_->MyGlobalElements_int_.resize(numMyElements);
579 
580  // Build the hash table
581  if(GenerateHash) Data_->LIDHash_int_ = new Epetra_HashTable<int>(numMyElements + 1 );
582  for(int i=0; i < numMyElements; ++i ) {
583  Data_->MyGlobalElements_int_[i] = myGlobalElements[i];
584  if(GenerateHash) Data_->LIDHash_int_->Add(myGlobalElements[i], i);
585  }
586  IsLongLong = false;
587  IsInt = true;
588 }
589 #endif
590 
591 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
592 void LightweightMap::Construct_LL(long long /* numGlobalElements */, int numMyElements, const long long * myGlobalElements, long long /* indexBase */, bool GenerateHash)
593 {
594  Data_=new LightweightMapData();
595  Data_->MyGlobalElements_LL_.resize(numMyElements);
596 
597  // Build the hash table
598  if(GenerateHash) Data_->LIDHash_LL_ = new Epetra_HashTable<long long>(numMyElements + 1 );
599  for(int i=0; i < numMyElements; ++i ) {
600  Data_->MyGlobalElements_LL_[i] = myGlobalElements[i];
601  if(GenerateHash) Data_->LIDHash_LL_->Add(myGlobalElements[i], i);
602  }
603  IsLongLong = true;
604  IsInt = false;
605 }
606 #endif
607 
608 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
609 LightweightMap::LightweightMap(int numGlobalElements,int numMyElements, const int * myGlobalElements, int indexBase, bool GenerateHash)
610 {
611  Construct_int(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
612 }
613 #endif
614 
615 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
616 LightweightMap::LightweightMap(long long numGlobalElements,int numMyElements, const long long * myGlobalElements, int indexBase, bool GenerateHash)
617 {
618  Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
619 }
620 
621 LightweightMap::LightweightMap(long long numGlobalElements,int numMyElements, const long long * myGlobalElements, long long indexBase, bool GenerateHash)
622 {
623  Construct_LL(numGlobalElements, numMyElements, myGlobalElements, indexBase, GenerateHash);
624 }
625 #endif
626 
627 //=========================================================================
628 LightweightMap::LightweightMap(const Epetra_Map & Map)
629 {
630  Data_=new LightweightMapData();
631  Data_->CopyMap_=new Epetra_Map(Map);
632  IsLongLong = Map.GlobalIndicesLongLong();
633  IsInt = Map.GlobalIndicesInt();
634 }
635 
636 //=========================================================================
638  : Data_(map.Data_)
639 {
640  Data_->IncrementReferenceCount();
641  IsLongLong = map.IsLongLong;
642  IsInt = map.IsInt;
643 }
644 
645 //=========================================================================
647  CleanupData();
648 }
649 
650 //=========================================================================
652 {
653  if((this != &map) && (Data_ != map.Data_)) {
654  CleanupData();
655  Data_ = map.Data_;
656  Data_->IncrementReferenceCount();
657  }
658  IsLongLong = map.IsLongLong;
659  IsInt = map.IsInt;
660  return(*this);
661 }
662 
663 //=========================================================================
664 void LightweightMap::CleanupData(){
665  if(Data_){
666  Data_->DecrementReferenceCount();
667  if(Data_->ReferenceCount() == 0) {
668  delete Data_;
669  }
670  }
671 }
672 
673 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
674 void LightweightMap::MyGlobalElementsPtr(int *& MyGlobalElementList) const {
675  MyGlobalElementList = Data_->MyGlobalElements_int_.size() ? &Data_->MyGlobalElements_int_.front() : 0;
676 }
677 #endif
678 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
679 void LightweightMap::MyGlobalElementsPtr(long long *& MyGlobalElementList) const {
680  MyGlobalElementList = Data_->MyGlobalElements_LL_.size() ? &Data_->MyGlobalElements_LL_.front() : 0;
681 }
682 #endif
683 
684 //=========================================================================
686  if(Data_->CopyMap_) return Data_->CopyMap_->NumMyElements();
687 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
688  if(IsInt)
689  return Data_->MyGlobalElements_int_.size();
690  else
691 #endif
692 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
693  if(IsLongLong)
694  return Data_->MyGlobalElements_LL_.size();
695  else
696 #endif
697  throw "EpetraExt::LightweightMap::NumMyElements: Global indices unknowns";
698 }
699 
700 //=========================================================================
701 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
702 int LightweightMap::LID(int gid) const {
703  if(Data_->CopyMap_) return Data_->CopyMap_->LID(gid);
704  if(IsInt)
705  return Data_->LIDHash_int_->Get(gid);
706  else if(IsLongLong)
707  throw "EpetraExt::LightweightMap::LID: Int version called for long long map";
708  else
709  throw "EpetraExt::LightweightMap::LID: unknown GID type";
710 }
711 #endif
712 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
713 int LightweightMap::LID(long long gid) const {
714 
715  if(Data_->CopyMap_) return Data_->CopyMap_->LID(gid);
716  if(IsLongLong)
717  return Data_->LIDHash_LL_->Get(gid);
718  else if(IsInt)
719  throw "EpetraExt::LightweightMap::LID: Long long version called for int map";
720  else
721  throw "EpetraExt::LightweightMap::LID: unknown GID type";
722 }
723 #endif
724 
725 //=========================================================================
726 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
727 int LightweightMap::GID(int lid) const {
728  if(Data_->CopyMap_) return Data_->CopyMap_->GID(lid);
729  if(lid < 0 || lid > (int)Data_->MyGlobalElements_int_.size()) return -1;
730  return Data_->MyGlobalElements_int_[lid];
731 }
732 #endif
733 long long LightweightMap::GID64(int lid) const {
734  if(Data_->CopyMap_) return Data_->CopyMap_->GID64(lid);
735 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
736  if(IsInt) {
737  if(lid < 0 || lid > (int)Data_->MyGlobalElements_int_.size()) return -1;
738  return Data_->MyGlobalElements_int_[lid];
739  }
740  else
741 #endif
742 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
743  if(IsLongLong) {
744  if(lid < 0 || lid > (int)Data_->MyGlobalElements_LL_.size()) return -1;
745  return Data_->MyGlobalElements_LL_[lid];
746  }
747  else
748 #endif
749  throw "EpetraExt::LightweightMap::GID64: Global indices unknown.";
750 }
751 
752 //=========================================================================
753 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
755  if(Data_->CopyMap_) return Data_->CopyMap_->MyGlobalElements();
756  else if(Data_->MyGlobalElements_int_.size()>0) return const_cast<int*>(&Data_->MyGlobalElements_int_[0]);
757  else return 0;
758 }
759 #endif
760 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
762  if(Data_->CopyMap_) return Data_->CopyMap_->MyGlobalElements64();
763  else if(Data_->MyGlobalElements_LL_.size()>0) return const_cast<long long*>(&Data_->MyGlobalElements_LL_[0]);
764  else return 0;
765 }
766 #endif
767 
768 //=========================================================================
770  if(Data_->CopyMap_) return Data_->CopyMap_->MinLID();
771  else return 0;
772 }
773 
774 //=========================================================================
776  if(Data_->CopyMap_) return Data_->CopyMap_->MaxLID();
777  else {
778 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
779  if(IsInt)
780  return Data_->MyGlobalElements_int_.size() - 1;
781  else
782 #endif
783 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
784  if(IsLongLong)
785  return Data_->MyGlobalElements_LL_.size() - 1;
786  else
787 #endif
788  throw "EpetraExt::LightweightMap::MaxLID: Global indices unknowns";
789  }
790 }
791 
792 
793 //=========================================================================
794 //=========================================================================
795 //=========================================================================
796 RemoteOnlyImport::RemoteOnlyImport(const Epetra_Import & Importer, LightweightMap & RemoteOnlyTargetMap)
797 {
798  int i;
799 
800  // Build an "Importer" that only takes the remote parts of the Importer.
801  SourceMap_=&Importer.SourceMap();
802  TargetMap_=&RemoteOnlyTargetMap;
803 
804  // Pull data from the Importer
805  NumSend_ = Importer.NumSend();
806  NumRemoteIDs_ = Importer.NumRemoteIDs();
807  NumExportIDs_ = Importer.NumExportIDs();
808  Distor_ = &Importer.Distributor();
809  int * OldRemoteLIDs = Importer.RemoteLIDs();
810  int * OldExportLIDs = Importer.ExportLIDs();
811  int * OldExportPIDs = Importer.ExportPIDs();
812 
813  // Sanity Check
814  if(NumRemoteIDs_ != RemoteOnlyTargetMap.NumMyElements())
815  throw std::runtime_error("RemoteOnlyImport: Importer doesn't match RemoteOnlyTargetMap for number of remotes.");
816 
817  // Copy the ExportIDs_, since they don't change
818  ExportLIDs_ = new int[NumExportIDs_];
819  ExportPIDs_ = new int[NumExportIDs_];
820  for(i=0; i<NumExportIDs_; i++) {
821  ExportLIDs_[i] = OldExportLIDs[i];
822  ExportPIDs_[i] = OldExportPIDs[i];
823  }
824 
825  // The RemoteIDs, on the other hand, do change. So let's do this right.
826  // Note: We might be able to bypass the LID call by just indexing off the Same and Permute GIDs, but at the moment this
827  // is fast enough not to worry about it.
828  RemoteLIDs_ = new int[NumRemoteIDs_];
829  if(TargetMap_->GlobalIndicesInt()) {
830  for(i=0; i<NumRemoteIDs_; i++)
831  RemoteLIDs_[i] = TargetMap_->LID( (int) Importer.TargetMap().GID64(OldRemoteLIDs[i]));
832  }
833  else if(TargetMap_->GlobalIndicesLongLong()) {
834  for(i=0; i<NumRemoteIDs_; i++)
835  RemoteLIDs_[i] = TargetMap_->LID(Importer.TargetMap().GID64(OldRemoteLIDs[i]));
836  }
837  else
838  throw std::runtime_error("RemoteOnlyImport: TargetMap_ index type unknown.");
839 
840  // Nowe we make sure these guys are in sorted order. AztecOO, ML and all that jazz.
841  for(i=0; i<NumRemoteIDs_-1; i++)
842  if(RemoteLIDs_[i] > RemoteLIDs_[i+1])
843  throw std::runtime_error("RemoteOnlyImport: Importer and RemoteOnlyTargetMap order don't match.");
844 }
845 
846 //=========================================================================
848 {
849  delete [] ExportLIDs_;
850  delete [] ExportPIDs_;
851  delete [] RemoteLIDs_;
852  // Don't delete the Distributor, SourceMap_ or TargetMap_ - those were shallow copies
853 }
854 
855 //=========================================================================
856 //=========================================================================
857 //=========================================================================
858 
859 template <class GO>
860 void MakeColMapAndReindexSort(int& NumRemoteColGIDs, GO*& RemoteColindices,
861  std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs);
862 
863 template <>
864 void MakeColMapAndReindexSort<int>(int& NumRemoteColGIDs, int*& RemoteColindices,
865  std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
866 {
867  // Sort External column indices so that all columns coming from a given remote processor are contiguous
868  int NLists = 2;
869  int* SortLists[3]; // this array is allocated on the stack, and so we won't need to delete it.
870  if(NumRemoteColGIDs > 0) {
871  SortLists[0] = RemoteColindices;
872  SortLists[1] = &RemotePermuteIDs[0];
873  Epetra_Util::Sort(true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, NLists, SortLists);
874  }
875 
876 
877  bool SortGhostsAssociatedWithEachProcessor_ = false;
878  if (SortGhostsAssociatedWithEachProcessor_) {
879  // Sort external column indices so that columns from a given remote processor are not only contiguous
880  // but also in ascending order. NOTE: I don't know if the number of externals associated
881  // with a given remote processor is known at this point ... so I count them here.
882 
883  NLists=1;
884  int StartCurrent, StartNext;
885  StartCurrent = 0; StartNext = 1;
886  while ( StartNext < NumRemoteColGIDs ) {
887  if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
888  else {
889  SortLists[0] = &RemotePermuteIDs[StartCurrent];
890  Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists);
891  StartCurrent = StartNext; StartNext++;
892  }
893  }
894  SortLists[0] = &RemotePermuteIDs[StartCurrent];
895  Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists);
896  }
897 }
898 
899 template <>
900 void MakeColMapAndReindexSort<long long>(int& NumRemoteColGIDs, long long*& RemoteColindices,
901  std::vector<int>& RemotePermuteIDs, std::vector<int>& RemoteOwningPIDs)
902 {
903  // Sort External column indices so that all columns coming from a given remote processor are contiguous
904  if(NumRemoteColGIDs > 0) {
905  long long* SortLists_LL[1] = {RemoteColindices};
906  int* SortLists_int[1] = {&RemotePermuteIDs[0]};
907  Epetra_Util::Sort(true, NumRemoteColGIDs, &RemoteOwningPIDs[0], 0, 0, 1, SortLists_int, 1, SortLists_LL);
908  }
909 
910  bool SortGhostsAssociatedWithEachProcessor_ = false;
911  if (SortGhostsAssociatedWithEachProcessor_) {
912  // Sort external column indices so that columns from a given remote processor are not only contiguous
913  // but also in ascending order. NOTE: I don't know if the number of externals associated
914  // with a given remote processor is known at this point ... so I count them here.
915 
916  int NLists=1;
917  int StartCurrent, StartNext;
918  StartCurrent = 0; StartNext = 1;
919  while ( StartNext < NumRemoteColGIDs ) {
920  if (RemoteOwningPIDs[StartNext]==RemoteOwningPIDs[StartNext-1]) StartNext++;
921  else {
922  int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
923  Epetra_Util::Sort(true,StartNext-StartCurrent, &(RemoteColindices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
924  StartCurrent = StartNext; StartNext++;
925  }
926  }
927  int* SortLists[1] = {&RemotePermuteIDs[StartCurrent]};
928  Epetra_Util::Sort(true, StartNext-StartCurrent, &(RemoteColindices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
929  }
930 }
931 
932 template <class GO>
933 int LightweightCrsMatrix::MakeColMapAndReindex(std::vector<int> owningPIDs, std::vector<GO> Gcolind)
934 {
935 #ifdef ENABLE_MMM_TIMINGS
936  using Teuchos::TimeMonitor;
937  Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.1")));
938 #endif
939 
940  // Scan all column indices and sort into two groups:
941  // Local: those whose GID matches a GID of the domain map on this processor and
942  // Remote: All others.
943  int numDomainElements = DomainMap_.NumMyElements();
944  std::vector<bool> LocalGIDs(numDomainElements,false);
945 
946  bool DoSizes = !DomainMap_.ConstantElementSize(); // If not constant element size, then error
947  if(DoSizes) EPETRA_CHK_ERR(-1);
948 
949  // In principle it is good to have RemoteGIDs and RemoteGIDList be as long as the number of remote GIDs
950  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
951  // and the number of block rows.
952  int numMyBlockRows;
953  if(use_lw) numMyBlockRows = RowMapLW_->NumMyElements();
954  else numMyBlockRows = RowMapEP_->NumMyElements();
955 
956  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
957  Epetra_HashTable<GO> RemoteGIDs(hashsize);
958  std::vector<GO> RemoteGIDList; RemoteGIDList.reserve(hashsize);
959  std::vector<int> RemoteOwningPIDs; RemoteOwningPIDs.reserve(hashsize);
960 
961  // In order to do the map reindexing inexpensively, we clobber the GIDs during this pass. For *local* GID's we clobber them
962  // with their LID in the domainMap. For *remote* GIDs, we clobber them with (numDomainElements+NumRemoteColGIDs) before the increment of
963  // the remote count. These numberings will be separate because no local LID is greater than numDomainElements.
964  int NumLocalColGIDs = 0;
965  int NumRemoteColGIDs = 0;
966  for(int i = 0; i < numMyBlockRows; i++) {
967  for(int j = rowptr_[i]; j < rowptr_[i+1]; j++) {
968  GO GID = Gcolind[j];
969  // Check if GID matches a row GID
970  int LID = DomainMap_.LID(GID);
971  if(LID != -1) {
972  bool alreadyFound = LocalGIDs[LID];
973  if (!alreadyFound) {
974  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
975  NumLocalColGIDs++;
976  }
977  colind_[j] = LID;
978  }
979  else {
980  GO hash_value=RemoteGIDs.Get(GID);
981  if(hash_value == -1) { // This means its a new remote GID
982  int PID = owningPIDs[j];
983  if(PID==-1) printf("[%d] ERROR: Remote PID should not be -1\n",DomainMap_.Comm().MyPID());
984  colind_[j] = numDomainElements + NumRemoteColGIDs;
985  RemoteGIDs.Add(GID, NumRemoteColGIDs);
986  RemoteGIDList.push_back(GID);
987  RemoteOwningPIDs.push_back(PID);
988  NumRemoteColGIDs++;
989  }
990  else
991  colind_[j] = numDomainElements + hash_value;
992  }
993  }
994  }
995 
996  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
997  if (DomainMap_.Comm().NumProc()==1) {
998  if (NumRemoteColGIDs!=0) {
999  throw "Some column IDs are not in domainMap. If matrix is rectangular, you must pass in domainMap to FillComplete";
1000  // Sanity test: When one processor,there can be no remoteGIDs
1001  }
1002  if (NumLocalColGIDs==numDomainElements) {
1003  ColMap_ = DomainMap_;
1004 
1005  // In this case, we just use the domainMap's indices, which is, not coincidently, what we clobbered colind_ with up above anyway.
1006  // No further reindexing is needed.
1007  return(0);
1008  }
1009  }
1010 
1011  // Now build integer array containing column GIDs
1012  // Build back end, containing remote GIDs, first
1013  int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1014  typename Epetra_GIDTypeSerialDenseVector<GO>::impl Colindices;
1015  if(numMyBlockCols > 0)
1016  Colindices.Size(numMyBlockCols);
1017  GO* RemoteColindices = Colindices.Values() + NumLocalColGIDs; // Points to back end of Colindices
1018 
1019  for(int i = 0; i < NumRemoteColGIDs; i++)
1020  RemoteColindices[i] = RemoteGIDList[i];
1021 
1022  // Build permute array for *remote* reindexing.
1023  std::vector<int> RemotePermuteIDs(NumRemoteColGIDs);
1024  for(int i=0; i<NumRemoteColGIDs; i++) RemotePermuteIDs[i]=i;
1025 
1026  MakeColMapAndReindexSort<GO>(NumRemoteColGIDs, RemoteColindices, RemotePermuteIDs, RemoteOwningPIDs);
1027 
1028  // Reverse the permutation to get the information we actually care about
1029  std::vector<int> ReverseRemotePermuteIDs(NumRemoteColGIDs);
1030  for(int i=0; i<NumRemoteColGIDs; i++) ReverseRemotePermuteIDs[RemotePermuteIDs[i]]=i;
1031 
1032  // Build permute array for *local* reindexing.
1033  bool use_local_permute=false;
1034  std::vector<int> LocalPermuteIDs(numDomainElements);
1035 
1036  // Now fill front end. Two cases:
1037  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
1038  // can simply read the domain GIDs into the front part of Colindices, otherwise
1039  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
1040  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
1041 
1042  if(NumLocalColGIDs == DomainMap_.NumMyElements()) {
1043  DomainMap_.MyGlobalElements(Colindices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
1044  }
1045  else {
1046  GO* MyGlobalElements = 0;
1047  DomainMap_.MyGlobalElementsPtr(MyGlobalElements);
1048  int NumLocalAgain = 0;
1049  use_local_permute = true;
1050  for(int i = 0; i < numDomainElements; i++) {
1051  if(LocalGIDs[i]) {
1052  LocalPermuteIDs[i] = NumLocalAgain;
1053  Colindices[NumLocalAgain++] = MyGlobalElements[i];
1054  }
1055  }
1056  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
1057  }
1058 
1059 
1060  // Copy the remote PID list correctly
1061  ColMapOwningPIDs_.resize(numMyBlockCols);
1062  ColMapOwningPIDs_.assign(numMyBlockCols,DomainMap_.Comm().MyPID());
1063  for(int i=0;i<NumRemoteColGIDs;i++)
1064  ColMapOwningPIDs_[NumLocalColGIDs+i] = RemoteOwningPIDs[i];
1065 
1066 #ifdef ENABLE_MMM_TIMINGS
1067  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.2")));
1068 #endif
1069 
1070  // Make Column map with same element sizes as Domain map
1071  LightweightMap temp((GO) -1, numMyBlockCols, Colindices.Values(), (GO) DomainMap_.IndexBase64());
1072  ColMap_ = temp;
1073 
1074 
1075 #ifdef ENABLE_MMM_TIMINGS
1076  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS-3.3")));
1077 #endif
1078 
1079  // Low-cost reindex of the matrix
1080  for(int i=0; i<numMyBlockRows; i++){
1081  for(int j=rowptr_[i]; j<rowptr_[i+1]; j++){
1082  int ID=colind_[j];
1083  if(ID < numDomainElements){
1084  if(use_local_permute) colind_[j] = LocalPermuteIDs[colind_[j]];
1085  // In the case where use_local_permute==false, we just copy the DomainMap's ordering, which it so happens
1086  // is what we put in colind_ to begin with.
1087  }
1088  else
1089  colind_[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_[j]-numDomainElements];
1090  }
1091  }
1092 
1093  assert((size_t)ColMap_.NumMyElements() == ColMapOwningPIDs_.size());
1094 
1095  return(0);
1096 }
1097 
1098 
1099 // Unused, file-local function doesn't need to exist.
1100 #if 0
1101 //=========================================================================
1102 // Template params are <PID,GID>
1103 static inline bool lessthan12(std::pair<int,int> i, std::pair<int,int> j){
1104  return ((i.first<j.first) || (i.first==j.first && i.second <j.second));
1105 }
1106 #endif // 0
1107 
1108 
1109 template<typename ImportType, typename int_type>
1110 int LightweightCrsMatrix::PackAndPrepareReverseComm(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter,
1111  std::vector<int> &ReverseSendSizes, std::vector<int_type> &ReverseSendBuffer) {
1112 #ifdef HAVE_MPI
1113  // Buffer pairs are in (PID,GID) order
1114  int i,j,k;
1115  const Epetra_Import *MyImporter= SourceMatrix.Importer();
1116  if(MyImporter == 0) return -1;
1117  const Epetra_MpiComm * MpiComm = dynamic_cast<const Epetra_MpiComm*>(&SourceMatrix.Comm());
1118  int MyPID = MpiComm->MyPID();
1119 
1120  // Things related to messages I and sending in forward mode (RowImporter)
1121  int NumExportIDs = RowImporter.NumExportIDs();
1122  int* ExportLIDs = RowImporter.ExportLIDs();
1123  int* ExportPIDs = RowImporter.ExportPIDs();
1124 
1125  // Things related to messages I am sending in reverse mode (MyImporter)
1126  Epetra_Distributor& Distor = MyImporter->Distributor();
1127  const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
1128  int NumRecvs = MDistor->NumReceives();
1129  int* RemoteLIDs = MyImporter->RemoteLIDs();
1130  const int * ProcsFrom = MDistor->ProcsFrom();
1131  const int * LengthsFrom = MDistor->LengthsFrom();
1132 
1133 
1134  // Get the owning pids in a special way, s.t. ProcsFrom[RemotePIDs[i]] is the guy who owns
1135  // RemoteLIDs[j]....
1136  std::vector<int> RemotePIDOrder(SourceMatrix.NumMyCols(),-1);
1137 
1138  // Now, for each remote ID, record which processor (in ProcsFrom ordering) owns it.
1139  for(i=0,j=0;i<NumRecvs;i++){
1140  for(k=0;k<LengthsFrom[i];k++){
1141  int pid=ProcsFrom[i];
1142  if(pid!=MyPID) RemotePIDOrder[RemoteLIDs[j]]=i;
1143  j++;
1144  }
1145  }
1146 
1147  // Step One: Start tacking the (GID,PID) pairs on the std sets
1148  std::vector<std::set<std::pair<int,int_type> > > ReversePGIDs(NumRecvs);
1149  int *rowptr, *colind;
1150  double *vals;
1151  EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1152 
1153 
1154  // Loop over each exported row and add to the temp list
1155  for(i=0; i < NumExportIDs; i++) {
1156  int lid = ExportLIDs[i];
1157  int exp_pid = ExportPIDs[i];
1158  for(j=rowptr[lid]; j<rowptr[lid+1]; j++){
1159  int pid_order = RemotePIDOrder[colind[j]];
1160  if(pid_order!=-1) {
1161  int_type gid = (int_type) SourceMatrix.GCID64(colind[j]);
1162  // This GID is getting shipped off somewhere
1163  ReversePGIDs[pid_order].insert(std::pair<int,int_type>(exp_pid,gid));
1164  }
1165  }
1166  }
1167 
1168  // Step 2: Count sizes (one too large to avoid std::vector errors)
1169  ReverseSendSizes.resize(NumRecvs+1);
1170  int totalsize=0;
1171  for(i=0; i<NumRecvs; i++) {
1172  ReverseSendSizes[i] = 2*ReversePGIDs[i].size();
1173  totalsize += ReverseSendSizes[i];
1174  }
1175 
1176  // Step 3: Alloc and fill the send buffer (one too large to avoid std::vector errors)
1177  ReverseSendBuffer.resize(totalsize+1);
1178  for(i=0, j=0; i<NumRecvs; i++) {
1179  for(typename std::set<std::pair<int,int_type> >::iterator it=ReversePGIDs[i].begin(); it!=ReversePGIDs[i].end(); it++) {
1180  ReverseSendBuffer[j] = it->first;
1181  ReverseSendBuffer[j+1] = it->second;
1182  j+=2;
1183  }
1184  }
1185 #endif
1186 
1187  return 0;
1188 }
1189 
1190 //=========================================================================
1191 
1192 template<typename int_type>
1193 void build_type3_exports_sort(std::vector<int_type>& ExportGID3, std::vector<int> &ExportPID3, int total_length3);
1194 
1195 template<>
1196 void build_type3_exports_sort<int>(std::vector<int>& ExportGID3, std::vector<int> &ExportPID3, int total_length3)
1197 {
1198  int * companion = &ExportGID3[0];
1199  Epetra_Util::Sort(true,total_length3,&ExportPID3[0],0,0,1,&companion,0,0);
1200 }
1201 
1202 template<>
1203 void build_type3_exports_sort<long long>(std::vector<long long>& ExportGID3, std::vector<int> &ExportPID3, int total_length3)
1204 {
1205  long long * companion = &ExportGID3[0];
1206  Epetra_Util::Sort(true,total_length3,&ExportPID3[0],0,0,0,0,1,&companion);
1207 }
1208 
1209 template<typename int_type>
1210 int build_type3_exports(int MyPID,int Nrecv, Epetra_BlockMap &DomainMap, std::vector<int> &ReverseRecvSizes, const int_type *ReverseRecvBuffer, std::vector<int> &ExportLID3, std::vector<int> &ExportPID3){
1211  int i,j;
1212 
1213  // Estimate total length of procs_to for Type 3
1214  int total_length3=0;
1215  for(i=0; i<Nrecv; i++)
1216  total_length3+=ReverseRecvSizes[i]/2;
1217  if(total_length3==0) return 0;
1218 
1219  std::vector<int_type> ExportGID3(total_length3);
1220  ExportLID3.resize(total_length3);
1221  ExportPID3.resize(total_length3);
1222 
1223  // Build a sorted colmap-style list for Type3 (removing any self-sends)
1224  for(i=0,j=0; i<2*total_length3; i+=2) {
1225  if(ReverseRecvBuffer[i] != MyPID){
1226  ExportPID3[j]=ReverseRecvBuffer[i];
1227  ExportGID3[j]=ReverseRecvBuffer[i+1];
1228  j++;
1229  }
1230  }
1231  total_length3=j;
1232 
1233  if(total_length3==0) return 0;
1234 
1235  // Sort (ala Epetra_CrsGraph)
1236  build_type3_exports_sort<int_type>(ExportGID3, ExportPID3, total_length3);
1237  int StartCurrent, StartNext;
1238  StartCurrent = 0; StartNext = 1;
1239  while ( StartNext < total_length3 ) {
1240  if(ExportPID3[StartNext] == ExportPID3[StartNext-1]) StartNext++;
1241  else {
1242  Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1243  StartCurrent = StartNext; StartNext++;
1244  }
1245  }
1246  Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID3[StartCurrent]),0,0,0,0, 0, 0);
1247 
1248 
1249  /* printf("[%d] Type 3 Sorted= ",MyComm.MyPID());
1250  for(i=0; i<total_length3; i++)
1251  printf("(--,%2d,%2d) ",ExportGID3[i],ExportPID3[i]);
1252  printf("\n");
1253  fflush(stdout);*/
1254 
1255 
1256  // Uniq & resize
1257  for(i=1,j=1; i<total_length3; i++){
1258  if(ExportPID3[i]!=ExportPID3[i-1] || ExportGID3[i]!=ExportGID3[i-1]){
1259  ExportPID3[j] = ExportPID3[i];
1260  ExportGID3[j] = ExportGID3[i];
1261  j++;
1262  }
1263  }
1264  ExportPID3.resize(j);
1265  ExportLID3.resize(j);
1266  total_length3=j;
1267 
1268  /* printf("[%d] Type 3 UNIQ = ",MyComm.MyPID());
1269  for(i=0; i<total_length3; i++)
1270  printf("(--,%2d,%2d) ",ExportGID3[i],ExportPID3[i]);
1271  printf("\n");
1272  fflush(stdout);*/
1273 
1274 
1275 
1276  // Now index down to LIDs
1277  for(i=0; i<total_length3; i++) {
1278  ExportLID3[i]=DomainMap.LID(ExportGID3[i]);
1279  if(ExportLID3[i] < 0) throw std::runtime_error("LightweightCrsMatrix:MakeExportLists invalid LID");
1280  }
1281 
1282  /* printf("[%d] Type 3 FINAL = ",MyComm.MyPID());
1283  for(i=0; i<total_length3; i++)
1284  printf("(%2d,%2d,%2d) ",ExportLID3[i],ExportGID3[i],ExportPID3[i]);
1285  printf("\n");
1286  fflush(stdout);*/
1287 
1288  return total_length3;
1289 }
1290 
1291 //=========================================================================
1292 template<typename ImportType, typename int_type>
1293 int build_type2_exports(const Epetra_CrsMatrix & SourceMatrix, ImportType & MyImporter, std::vector<int> &ExportLID2, std::vector<int> &ExportPID2){
1294  int total_length2=0;
1295 #ifdef HAVE_MPI
1296  int i,j;
1297  const Epetra_Import *SourceImporter= SourceMatrix.Importer();
1298 
1299  int *rowptr, *colind;
1300  double *vals;
1301  EPETRA_CHK_ERR(SourceMatrix.ExtractCrsDataPointers(rowptr,colind,vals));
1302 
1303  // Things related to messages I am sending in forward mode (MyImporter)
1304  int NumExportIDs = MyImporter.NumExportIDs();
1305  const int* ExportLIDs = MyImporter.ExportLIDs();
1306  const int* ExportPIDs = MyImporter.ExportPIDs();
1307  if(NumExportIDs==0) return 0;
1308 
1309  Epetra_Distributor& Distor = MyImporter.Distributor();
1310  const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
1311 
1312  // Assume I own all the cols, then flag any cols I don't own
1313  // This allows us to avoid LID calls later...
1314  std::vector<bool> IsOwned(SourceMatrix.NumMyCols(),true);
1315  if(SourceImporter) {
1316  const int * RemoteLIDs = SourceImporter->RemoteLIDs();
1317  // Now flag the cols I don't own
1318  for(i=0; i<SourceImporter->NumRemoteIDs(); i++)
1319  IsOwned[RemoteLIDs[i]]=false;
1320  }
1321 
1322  // Status vector
1323  std::vector<int> SentTo(SourceMatrix.NumMyCols(),-1);
1324 
1325  // Initial allocation: NumUnknowsSent * Max row size (guaranteed to be too large)
1326  for(i=0,total_length2=0; i<MDistor->NumSends(); i++) total_length2+= MDistor->LengthsTo()[i] * SourceMatrix.MaxNumEntries();
1327  std::vector<int_type> ExportGID2(total_length2);
1328 
1329  ExportLID2.resize(total_length2);
1330  ExportPID2.resize(total_length2);
1331 
1332  int current=0, last_start=0, last_pid=ExportPIDs[0];
1333  for(i=0; i<NumExportIDs; i++){
1334  // For each row I have to send via MyImporter...
1335  int row=ExportLIDs[i];
1336  int pid=ExportPIDs[i];
1337 
1338  if(i!=0 && pid>last_pid) {
1339  // We have a new PID, so lets finish up the current one
1340  if(current!=last_start){
1341  int *lids = &ExportLID2[last_start];
1342  Epetra_Util::Sort(true,current-last_start,&ExportGID2[last_start],0,0,1,&lids,0,0);
1343  // Note: we don't need to sort the ExportPIDs since they're the same since last_start
1344  }
1345  // Reset the list
1346  last_pid=pid;
1347  last_start=current;
1348  }
1349  else if(pid < last_pid) {
1350  throw std::runtime_error("build_type2_exports: ExportPIDs are not sorted!");
1351  }
1352 
1353  for(j=rowptr[row]; j<rowptr[row+1]; j++) {
1354  // For each column in that row...
1355  int col=colind[j];
1356  if(IsOwned[col] && SentTo[col]!=pid){
1357  // We haven't added this guy to the list yet.
1358  if(current>= total_length2) throw std::runtime_error("build_type2_exports: More export ids than I thought!");
1359  SentTo[col] = pid;
1360  ExportGID2[current] = (int_type) SourceMatrix.GCID64(col);
1361  ExportLID2[current] = SourceMatrix.DomainMap().LID(ExportGID2[current]);
1362  ExportPID2[current] = pid;
1363  current++;
1364  }
1365  }
1366  }//end main loop
1367 
1368  // Final Sort
1369  int *lids = ExportLID2.size() > (std::size_t) last_start ? &ExportLID2[last_start] : 0;
1370  Epetra_Util::Sort(true,current-last_start,ExportGID2.size() > (std::size_t) last_start ? &ExportGID2[last_start] : 0,0,0,1,&lids,0,0);
1371  // Note: we don't need to sort the ExportPIDs since they're the same since last_start
1372 
1373  total_length2=current;
1374  ExportLID2.resize(total_length2);
1375  ExportPID2.resize(total_length2);
1376 #endif
1377  return total_length2;
1378 }
1379 
1380 //=========================================================================
1381 template<typename int_type>
1382 void build_type1_exports_sort(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int_type>& ExportGID1, int total_length1);
1383 
1384 template<>
1385 void build_type1_exports_sort<int>(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<int>& ExportGID1, int total_length1){
1386  int * companion[2] = {&ExportLID1[0],&ExportGID1[0]};
1387  Epetra_Util::Sort(true,total_length1,&ExportPID1[0],0,0,2,&companion[0],0,0);
1388 }
1389 
1390 template<>
1391 void build_type1_exports_sort<long long>(std::vector<int> &ExportLID1, std::vector<int> &ExportPID1, std::vector<long long>& ExportGID1, int total_length1){
1392  int * companion = &ExportLID1[0];
1393  long long * companion64 = &ExportGID1[0];
1394  Epetra_Util::Sort(true,total_length1,&ExportPID1[0],0,0,1,&companion,1,&companion64);
1395 }
1396 
1397 template<typename int_type>
1398 int build_type1_exports(const Epetra_Import * Importer1, std::vector<int> &ExportLID1, std::vector<int> &ExportPID1){
1399  int i, total_length1=0;
1400  if(!Importer1) return 0;
1401  total_length1 = Importer1->NumSend();
1402  if(total_length1==0) return 0;
1403 
1404  std::vector<int_type> ExportGID1(total_length1);
1405  ExportLID1.resize(total_length1);
1406  ExportPID1.resize(total_length1);
1407  const int * ExportLID1Base = Importer1->ExportLIDs();
1408  const int * ExportPID1Base = Importer1->ExportPIDs();
1409 
1410  for(i=0; i<total_length1; i++){
1411  ExportLID1[i] = ExportLID1Base[i];
1412  ExportPID1[i] = ExportPID1Base[i];
1413  ExportGID1[i] = (int_type) Importer1->SourceMap().GID64(ExportLID1Base[i]);
1414  }
1415 
1416  // Sort (ala Epetra_CrsGraph)
1417  build_type1_exports_sort<int_type>(ExportLID1, ExportPID1, ExportGID1, total_length1);
1418 
1419  int StartCurrent, StartNext;
1420  StartCurrent = 0; StartNext = 1;
1421  while ( StartNext < total_length1 ) {
1422  if(ExportPID1[StartNext] == ExportPID1[StartNext-1]) StartNext++;
1423  else {
1424  int *new_companion = {&ExportLID1[StartCurrent]};
1425  Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1426  StartCurrent = StartNext; StartNext++;
1427  }
1428  }
1429  int *new_companion = {&ExportLID1[StartCurrent]};
1430  Epetra_Util::Sort(true,StartNext-StartCurrent, &(ExportGID1[StartCurrent]),0,0,1,&new_companion, 0, 0);
1431  return total_length1;
1432 }
1433 
1434 //=========================================================================
1435 template<typename ImportType, typename int_type>
1436 int LightweightCrsMatrix::MakeExportLists(const Epetra_CrsMatrix & SourceMatrix, ImportType & Importer2,
1437  std::vector<int> &ReverseRecvSizes, const int_type *ReverseRecvBuffer,
1438  std::vector<int> & ExportPIDs, std::vector<int> & ExportLIDs) {
1439 #ifdef HAVE_MPI
1440  int MyPID = SourceMatrix.Comm().MyPID();
1441 
1442  // This could (legitimately) be zero, in which case we don't have any ReverseRecvs either.
1443  const Epetra_Import *Importer1= SourceMatrix.Importer();
1444 
1445  // So for each entry in the DomainMap, I have to answer the question: Do I need to send this to anybody? And if so, to whom?
1446  //
1447  // This information comes from three sources:
1448  // 1) IDs in my DomainMap that are in someone else's ColMap (e.g. SourceMatrix.Importer()).
1449  // 2) IDs that I own that I sent to another proc in the Forward communication round (e.g. RowImporter).
1450  // 3) IDs that someone else sent on during the Forward round (and they told me about duing the reverse round).
1451  //
1452  // Any of these could legitimately be null.
1453  Epetra_MpiDistributor * Distor1 = (Importer1)?(dynamic_cast<Epetra_MpiDistributor*>(&Importer1->Distributor())):0;
1454 
1455  int Nsend1 = (Distor1)?(Distor1->NumSends()):0; // Also the number of messages we'll need to parse through in build_type3_exports
1456 
1457  std::vector<int> ExportPID3;
1458  std::vector<int> ExportLID3;
1459 
1460  std::vector<int> ExportPID2;
1461  std::vector<int> ExportLID2;
1462 
1463  std::vector<int> ExportPID1;
1464  std::vector<int> ExportLID1;
1465 
1466  // Build (sorted) ExportLID / ExportGID list for Type
1467  int Len1=build_type1_exports<int_type>(Importer1, ExportLID1, ExportPID1);
1468  int Len2=build_type2_exports<ImportType, int_type>(SourceMatrix, Importer2, ExportLID2, ExportPID2);
1469  int Len3=build_type3_exports(MyPID,Nsend1,DomainMap_,ReverseRecvSizes, ReverseRecvBuffer, ExportLID3, ExportPID3);
1470 
1471  // Since everything should now be sorted correctly, we can do a streaming merge of the three Export lists...
1472 #ifdef HAVE_EPETRAEXT_DEBUG
1473  {
1474  int i;
1475  // Now we sanity check the 1 & 2 lists
1476  bool test_passed=true;
1477  for(i=1; i<Len1; i++) {
1478  if(ExportPID1[i] < ExportPID1[i-1] || (ExportPID1[i] == ExportPID1[i-1] && DomainMap_.GID(ExportLID1[i]) < DomainMap_.GID(ExportLID1[i-1])))
1479  test_passed=false;
1480  }
1481  SourceMatrix.Comm().Barrier();
1482  if(!test_passed) {
1483  printf("[%d] Type1 ERRORS = ",SourceMatrix.Comm().MyPID());
1484  for(int i=0; i<Len1; i++)
1485  printf("(%2d,%2d,%2d) ",ExportLID1[i],DomainMap_.GID(ExportLID1[i]),ExportPID1[i]);
1486  printf("\n");
1487  fflush(stdout);
1488  throw std::runtime_error("Importer1 fails the sanity test");
1489  }
1490 
1491  for(i=1; i<Len2; i++) {
1492  if(ExportPID2[i] < ExportPID2[i-1] || (ExportPID2[i] == ExportPID2[i-1] && DomainMap_.GID(ExportLID2[i]) < DomainMap_.GID(ExportLID2[i-1])))
1493  test_passed=false;
1494  }
1495 
1496  SourceMatrix.Comm().Barrier();
1497  if(!test_passed) {
1498  printf("[%d] Type2 ERRORS = ",SourceMatrix.Comm().MyPID());
1499  for(int i=0; i<Len2; i++)
1500  printf("(%2d,%2d,%2d) ",ExportLID2[i],DomainMap_.GID(ExportLID2[i]),ExportPID2[i]);
1501  printf("\n");
1502  fflush(stdout);
1503  throw std::runtime_error("Importer2 fails the sanity test");
1504  }
1505 
1506  for(i=1; i<Len3; i++) {
1507  if(ExportPID3[i] < ExportPID3[i-1] || (ExportPID3[i] == ExportPID3[i-1] && DomainMap_.GID(ExportLID3[i]) < DomainMap_.GID(ExportLID3[i-1])))
1508  test_passed=false;
1509  }
1510 
1511  SourceMatrix.Comm().Barrier();
1512  if(!test_passed) {
1513  printf("[%d] Type3 ERRORS = ",SourceMatrix.Comm().MyPID());
1514  for(int i=0; i<Len3; i++)
1515  printf("(%2d,%2d,%2d) ",ExportLID3[i],DomainMap_.GID(ExportLID3[i]),ExportPID3[i]);
1516  printf("\n");
1517  fflush(stdout);
1518  throw std::runtime_error("Importer3 fails the sanity test");
1519  }
1520 
1521 
1522  }
1523 #endif
1524 
1525  if(Importer1 && !Importer1->SourceMap().SameAs(DomainMap_))
1526  throw std::runtime_error("ERROR: Map Mismatch Importer1");
1527 
1528  if(!Importer2.SourceMap().SameAs(SourceMatrix.RowMap()))
1529  throw std::runtime_error("ERROR: Map Mismatch Importer2");
1530 
1531  int_type InfGID = std::numeric_limits<int_type>::max();
1532  int InfPID = INT_MAX;
1533 
1534  int i1=0, i2=0, i3=0, current=0;
1535 
1536  int MyLen=Len1+Len2+Len3;
1537  ExportLIDs.resize(MyLen);
1538  ExportPIDs.resize(MyLen);
1539 
1540  while(i1 < Len1 || i2 < Len2 || i3 < Len3){
1541  int PID1 = (i1<Len1)?(ExportPID1[i1]):InfPID;
1542  int PID2 = (i2<Len2)?(ExportPID2[i2]):InfPID;
1543  int PID3 = (i3<Len3)?(ExportPID3[i3]):InfPID;
1544 
1545  int_type GID1 = (i1<Len1)?((int_type) DomainMap_.GID64(ExportLID1[i1])):InfGID;
1546  int_type GID2 = (i2<Len2)?((int_type) DomainMap_.GID64(ExportLID2[i2])):InfGID;
1547  int_type GID3 = (i3<Len3)?((int_type) DomainMap_.GID64(ExportLID3[i3])):InfGID;
1548 
1549  int MIN_PID = MIN3(PID1,PID2,PID3);
1550  int_type MIN_GID = MIN3( ((PID1==MIN_PID)?GID1:InfGID), ((PID2==MIN_PID)?GID2:InfGID), ((PID3==MIN_PID)?GID3:InfGID));
1551  bool added_entry=false;
1552 
1553  // Case 1: Add off list 1
1554  if(PID1 == MIN_PID && GID1 == MIN_GID){
1555  ExportLIDs[current] = ExportLID1[i1];
1556  ExportPIDs[current] = ExportPID1[i1];
1557  current++;
1558  i1++;
1559  added_entry=true;
1560  }
1561 
1562  // Case 2: Add off list 2
1563  if(PID2 == MIN_PID && GID2 == MIN_GID){
1564  if(!added_entry) {
1565  ExportLIDs[current] = ExportLID2[i2];
1566  ExportPIDs[current] = ExportPID2[i2];
1567  current++;
1568  added_entry=true;
1569  }
1570  i2++;
1571  }
1572 
1573  // Case 3: Add off list 3
1574  if(PID3 == MIN_PID && GID3 == MIN_GID){
1575  if(!added_entry) {
1576  ExportLIDs[current] = ExportLID3[i3];
1577  ExportPIDs[current] = ExportPID3[i3];
1578  current++;
1579  }
1580  i3++;
1581  }
1582  }// end while
1583  if(current!=MyLen) {
1584  ExportLIDs.resize(current);
1585  ExportPIDs.resize(current);
1586  }
1587 #endif
1588  return 0;
1589 }
1590 
1591 //=========================================================================
1592 template<typename ImportType, typename int_type>
1593 void LightweightCrsMatrix::Construct(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter)
1594 {
1595  // Do we need to use long long for GCIDs?
1596 
1597 #ifdef ENABLE_MMM_TIMINGS
1598  using Teuchos::TimeMonitor;
1599  Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-1")));
1600 #endif
1601 
1602  // Fused constructor, import & FillComplete
1603  int rv=0;
1604  int N;
1605  if(use_lw) N = RowMapLW_->NumMyElements();
1606  else N = RowMapEP_->NumMyElements();
1607 
1608  int MyPID = SourceMatrix.Comm().MyPID();
1609 
1610 #ifdef HAVE_MPI
1611  std::vector<int> ReverseSendSizes;
1612  std::vector<int_type> ReverseSendBuffer;
1613  std::vector<int> ReverseRecvSizes;
1614  int_type * ReverseRecvBuffer=0;
1615 #endif
1616 
1617  bool communication_needed = RowImporter.SourceMap().DistributedGlobal();
1618 
1619  // The basic algorithm here is:
1620  // 1) Call Distor.Do to handle the import.
1621  // 2) Copy all the Imported and Copy/Permuted data into the raw Epetra_CrsMatrix / Epetra_CrsGraphData pointers, still using GIDs.
1622  // 3) Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids) AND
1623  // reindexes to LIDs.
1624 
1625  // Sanity Check
1626  if (!SourceMatrix.RowMap().SameAs(RowImporter.SourceMap()))
1627  throw "LightweightCrsMatrix: Fused copy constructor requires Importer.SourceMap() to match SourceMatrix.RowMap()";
1628 
1629  // Get information from the Importer
1630  int NumSameIDs = RowImporter.NumSameIDs();
1631  int NumPermuteIDs = RowImporter.NumPermuteIDs();
1632  int NumRemoteIDs = RowImporter.NumRemoteIDs();
1633  int NumExportIDs = RowImporter.NumExportIDs();
1634  int* ExportLIDs = RowImporter.ExportLIDs();
1635  int* RemoteLIDs = RowImporter.RemoteLIDs();
1636  int* PermuteToLIDs = RowImporter.PermuteToLIDs();
1637  int* PermuteFromLIDs = RowImporter.PermuteFromLIDs();
1638 
1639 #ifdef HAVE_MPI
1640  Epetra_Distributor& Distor = RowImporter.Distributor();
1641  const Epetra_MpiComm * MpiComm = dynamic_cast<const Epetra_MpiComm*>(&SourceMatrix.Comm());
1642  const Epetra_MpiDistributor * MDistor = dynamic_cast<Epetra_MpiDistributor*>(&Distor);
1643 #endif
1644 
1645  // Allocate memory
1646  rowptr_.resize(N+1);
1647 
1648 
1649  // Owning PIDs
1650  std::vector<int> SourcePids;
1651  std::vector<int> TargetPids;
1652 
1653 
1654  /***************************************************/
1655  /***** 1) From Epetra_DistObject::DoTransfer() *****/
1656  /***************************************************/
1657  // rv=SourceMatrix.CheckSizes(SourceMatrix);
1658 
1659  // NTS: Add CheckSizes stuff here.
1660  if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in CheckSizes()";
1661 
1662  // Buffers & Other Relevant Info
1663  char* Exports_ = 0;
1664  char* Imports_ = 0;
1665  int LenImports_ =0;
1666  int LenExports_ = 0;
1667  int *Sizes_ = 0;
1668 
1669  int SizeOfPacket;
1670  bool VarSizes = false;
1671  if( NumExportIDs > 0) {
1672  Sizes_ = new int[NumExportIDs];
1673  }
1674 
1675 
1676  // Get the owning PIDs
1677  const Epetra_Import *MyImporter= SourceMatrix.Importer();
1678  if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,false);
1679  else {
1680  SourcePids.resize(SourceMatrix.ColMap().NumMyElements());
1681  SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),SourceMatrix.Comm().MyPID());
1682  }
1683 
1684  // Pack & Prepare w/ owning PIDs
1685  rv=Epetra_Import_Util::PackAndPrepareWithOwningPIDs(SourceMatrix,
1686  NumExportIDs,ExportLIDs,
1687  LenExports_,Exports_,SizeOfPacket,
1688  Sizes_,VarSizes,SourcePids);
1689  if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in PackAndPrepare()";
1690 
1691  if (communication_needed) {
1692 #ifdef HAVE_MPI
1693  // Do the exchange of remote data
1694  const int * ExportPIDs = RowImporter.ExportPIDs();
1695 
1696  // Use the fact that the export procs are sorted to avoid building a hash table.
1697  // NOTE: The +1's on the message size lists are to avoid std::vector problems if a proc has no sends or recvs.
1698  std::vector<int> SendSizes(MDistor->NumSends()+1,0);
1699  for(int i=0, curr_pid=0; i<NumExportIDs; i++) {
1700  if(i>0 && ExportPIDs[i] > ExportPIDs[i-1]) curr_pid++;
1701  SendSizes[curr_pid] +=Sizes_[i];
1702 
1703  // sanity check
1704  if(i>0 && ExportPIDs[i] < ExportPIDs[i-1]) throw "ExportPIDs not sorted";
1705  }
1706 
1707  std::vector<int> RecvSizes(MDistor->NumReceives()+1);
1708  int msg_tag=MpiComm->GetMpiTag();
1709  boundary_exchange_varsize<char>(*MpiComm,MPI_CHAR,MDistor->NumSends(),MDistor->ProcsTo(),SendSizes.size() ? &SendSizes[0] : 0,Exports_,
1710  MDistor->NumReceives(),MDistor->ProcsFrom(),RecvSizes.size() ? &RecvSizes[0] : 0,Imports_,SizeOfPacket,msg_tag);
1711 
1712  // If the source matrix doesn't have an importer, then nobody sent data belonging to me in the forward round.
1713  if(SourceMatrix.Importer()) {
1714  Epetra_Import* SourceImporter=const_cast<Epetra_Import*>(SourceMatrix.Importer());
1715  const Epetra_MpiDistributor * MyDistor = dynamic_cast<Epetra_MpiDistributor*>(&SourceImporter->Distributor());
1716 
1717  // Setup the reverse communication
1718  // Note: Buffer pairs are in (PID,GID) order
1719  PackAndPrepareReverseComm<ImportType, int_type>(SourceMatrix,RowImporter,ReverseSendSizes,ReverseSendBuffer);
1720 
1721  // Do the reverse communication
1722  // NOTE: Make the vector one too large to avoid std::vector errors
1723  ReverseRecvSizes.resize(MyDistor->NumSends()+1);
1724  const int msg_tag2 = MpiComm->GetMpiTag ();
1725  MPI_Datatype data_type = sizeof(int_type) == 4 ? MPI_INT : MPI_LONG_LONG;
1726  boundary_exchange_varsize<int_type> (*MpiComm, data_type, MyDistor->NumReceives (),
1727  MyDistor->ProcsFrom (),
1728  ReverseSendSizes.size() ? &ReverseSendSizes[0] : 0,
1729  ReverseSendBuffer.size() ? &ReverseSendBuffer[0] : 0,
1730  MyDistor->NumSends (), MyDistor->ProcsTo (),
1731  ReverseRecvSizes.size() ? &ReverseRecvSizes[0] : 0,
1732  ReverseRecvBuffer, 1, msg_tag2);
1733  }
1734 #endif
1735  }
1736 
1737  if(rv) throw "LightweightCrsMatrix: Fused copy constructor failed in Distor.Do";
1738 
1739  /*********************************************************************/
1740  /**** 2) Copy all of the Same/Permute/Remote data into CSR_arrays ****/
1741  /*********************************************************************/
1742 #ifdef ENABLE_MMM_TIMINGS
1743  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-2")));
1744 #endif
1745 
1746  // Count nnz
1747  int mynnz = Epetra_Import_Util::UnpackWithOwningPIDsCount(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_);
1748  // Presume the rowptr_ is the right size
1749  // Allocate colind_ & vals_ arrays
1750  colind_.resize(mynnz);
1751 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1752  if(sizeof(int_type) == sizeof(long long))
1753  colind_LL_.resize(mynnz);
1754 #endif
1755  vals_.resize(mynnz);
1756 
1757  // Reset the Source PIDs (now with -1 rule)
1758  if(MyImporter) Epetra_Util::GetPids(*MyImporter,SourcePids,true);
1759  else {
1760  SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),-1);
1761  }
1762 
1763  // Unpack into arrays
1764  int * myrowptr = rowptr_.size() ? & rowptr_[0] : 0;
1765  double * myvals = vals_.size() ? & vals_[0] : 0;
1766 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1767  if(sizeof(int_type) == sizeof(int)) {
1768  int * mycolind = colind_.size() ? & colind_[0] : 0;
1769  Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1770  }
1771  else
1772 #endif
1773 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1774  if(sizeof(int_type) == sizeof(long long)) {
1775  long long * mycolind = colind_LL_.size() ? & colind_LL_[0] : 0;
1776  Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,N,mynnz,MyPID,myrowptr,mycolind,myvals,SourcePids,TargetPids);
1777  }
1778  else
1779 #endif
1780  throw "EpetraExt::LightweightCrsMatrix::Construct: sizeof(int_type) error.";
1781 
1782  /**************************************************************/
1783  /**** 3) Call Optimized MakeColMap w/ no Directory Lookups ****/
1784  /**************************************************************/
1785 #ifdef ENABLE_MMM_TIMINGS
1786  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-3")));
1787 #endif
1788 
1789  //Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids).
1790  MakeColMapAndReindex<int_type>(TargetPids,getcolind<int_type>());
1791 
1792  /********************************************/
1793  /**** 4) Make Export Lists for Import ****/
1794  /********************************************/
1795 #ifdef HAVE_MPI
1796  MakeExportLists<ImportType, int_type>(SourceMatrix,RowImporter,ReverseRecvSizes,ReverseRecvBuffer,ExportPIDs_,ExportLIDs_);
1797 #endif
1798 
1799  /********************************************/
1800  /**** 5) Call sort the entries ****/
1801  /********************************************/
1802  // NOTE: If we have no entries the &blah[0] will cause the STL to die in debug mode
1803 #ifdef ENABLE_MMM_TIMINGS
1804  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-4")));
1805 #endif
1806  if(N>0) Epetra_Util::SortCrsEntries(N, &rowptr_[0], colind_.size() ? &colind_[0] : 0, vals_.size() ? &vals_[0] : 0);
1807 
1808  /********************************************/
1809  /**** 6) Cleanup ****/
1810  /********************************************/
1811 #ifdef ENABLE_MMM_TIMINGS
1812  MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS C-5")));
1813 #endif
1814 
1815 #ifdef HAVE_MPI
1816  delete [] ReverseRecvBuffer;
1817 #endif
1818 
1819  delete [] Exports_;
1820  delete [] Imports_;
1821  delete [] Sizes_;
1822 
1823  }// end fused copy constructor
1824 
1825 
1826 
1827 
1828 //=========================================================================
1829 LightweightCrsMatrix::LightweightCrsMatrix(const Epetra_CrsMatrix & SourceMatrix, RemoteOnlyImport & RowImporter):
1830  use_lw(true),
1831  RowMapLW_(0),
1832  RowMapEP_(0),
1833  DomainMap_(SourceMatrix.DomainMap())
1834 {
1835 #ifdef ENABLE_MMM_TIMINGS
1836  using Teuchos::TimeMonitor;
1837  Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: LWCRS Total")));
1838 #endif
1839 
1840  RowMapLW_= new LightweightMap(RowImporter.TargetMap());
1841 
1842 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1843  if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1844  Construct<RemoteOnlyImport, int>(SourceMatrix,RowImporter);
1845  }
1846  else
1847 #endif
1848 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1849  if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1850  Construct<RemoteOnlyImport, long long>(SourceMatrix,RowImporter);
1851  }
1852  else
1853 #endif
1854  throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1855 
1856 }
1857 
1858 
1859 //=========================================================================
1860 LightweightCrsMatrix::LightweightCrsMatrix(const Epetra_CrsMatrix & SourceMatrix, Epetra_Import & RowImporter):
1861  use_lw(false),
1862  RowMapLW_(0),
1863  RowMapEP_(0),
1864  DomainMap_(SourceMatrix.DomainMap())
1865 {
1866  RowMapEP_= new Epetra_BlockMap(RowImporter.TargetMap());
1867 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1868  if(SourceMatrix.RowMatrixRowMap().GlobalIndicesInt()) {
1869  Construct<Epetra_Import, int>(SourceMatrix,RowImporter);
1870  }
1871  else
1872 #endif
1873 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1874  if(SourceMatrix.RowMatrixRowMap().GlobalIndicesLongLong()) {
1875  Construct<Epetra_Import, long long>(SourceMatrix,RowImporter);
1876  }
1877  else
1878 #endif
1879  throw "EpetraExt::LightweightCrsMatrix: ERROR, GlobalIndices type unknown.";
1880 }
1881 
1882 
1884  delete RowMapLW_;
1885  delete RowMapEP_;
1886 }
1887 
1888 
1889 }//namespace EpetraExt
1890 
1891 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1893 #endif
1894 
1895 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1897 #endif
int SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
LightweightCrsMatrix * importMatrix
RemoteOnlyImport(const Epetra_Import &Importer, LightweightMap &RemoteOnlyTargetMap)
Epetra_HashTable< long long > * LIDHash_LL_
std::vector< int > targetMapToOrigRow
long long * MyGlobalElements64() const
void debug_print_distor(const char *label, const Epetra_Distributor *Distor, const Epetra_Comm &Comm)
std::pair< int_type, int_type > Tget_col_range(const Epetra_CrsMatrix &mtx)
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
void pack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int > &proc_col_ranges, std::vector< int > &send_rows, std::vector< int > &rows_per_send_proc)
void MakeColMapAndReindexSort< long long >(int &NumRemoteColGIDs, long long *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
void debug_compare_import(const Epetra_Import *Import1, const Epetra_Import *Import2)
long long GID64(int LID) const
void MakeColMapAndReindexSort< int >(int &NumRemoteColGIDs, int *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
const Epetra_CrsMatrix * origMatrix
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
void build_type1_exports_sort(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< int_type > &ExportGID1, int total_length1)
void build_type1_exports_sort< int >(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< int > &ExportGID1, int total_length1)
int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
std::pair< long long, long long > get_col_range< long long >(const Epetra_Map &emap)
std::vector< int > targetMapToImportRow
LightweightMap & operator=(const LightweightMap &map)
void build_type3_exports_sort< int >(std::vector< int > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
void build_type1_exports_sort< long long >(std::vector< int > &ExportLID1, std::vector< int > &ExportPID1, std::vector< long long > &ExportGID1, int total_length1)
void build_type3_exports_sort(std::vector< int_type > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
int build_type2_exports(const Epetra_CrsMatrix &SourceMatrix, ImportType &MyImporter, std::vector< int > &ExportLID2, std::vector< int > &ExportPID2)
int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
#define MIN3(x, y, z)
LightweightCrsMatrix(const Epetra_CrsMatrix &A, RemoteOnlyImport &RowImporter)
CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix &epetracrsmatrix)
std::vector< int > MyGlobalElements_int_
int build_type3_exports(int MyPID, int Nrecv, Epetra_BlockMap &DomainMap, std::vector< int > &ReverseRecvSizes, const int_type *ReverseRecvBuffer, std::vector< int > &ExportLID3, std::vector< int > &ExportPID3)
const LightweightMap & TargetMap() const
void Tpack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int_type > &proc_col_ranges, std::vector< int_type > &send_rows, std::vector< int > &rows_per_send_proc)
CrsWrapper_GraphBuilder(const Epetra_Map &emap)
Epetra_HashTable< int > * LIDHash_int_
std::vector< long long > MyGlobalElements_LL_
int build_type1_exports(const Epetra_Import *Importer1, std::vector< int > &ExportLID1, std::vector< int > &ExportPID1)
std::map< int_type, std::set< int_type > * > & get_graph()
void build_type3_exports_sort< long long >(std::vector< long long > &ExportGID3, std::vector< int > &ExportPID3, int total_length3)
const Epetra_BlockMap * importColMap
void insert_matrix_locations(CrsWrapper_GraphBuilder< int_type > &graphbuilder, Epetra_CrsMatrix &C)
void MakeColMapAndReindexSort(int &NumRemoteColGIDs, GO *&RemoteColindices, std::vector< int > &RemotePermuteIDs, std::vector< int > &RemoteOwningPIDs)
int dumpCrsMatrixStruct(const CrsMatrixStruct &M)
int InsertGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
std::pair< int, int > get_col_range< int >(const Epetra_Map &emap)