Epetra Package Browser (Single Doxygen Collection)  Development
Epetra_CrsGraph.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_CrsGraph.h"
46 #include "Epetra_Import.h"
47 #include "Epetra_Export.h"
48 #include "Epetra_Distributor.h"
49 #include "Epetra_Util.h"
50 #include "Epetra_Comm.h"
51 #include "Epetra_HashTable.h"
52 #include "Epetra_BlockMap.h"
53 #include "Epetra_Map.h"
54 #include "Epetra_RowMatrix.h"
56 #include "Epetra_SerialComm.h"
57 
58 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
60 #endif
61 
63 #include "Epetra_OffsetIndex.h"
64 
65 //==============================================================================
67  const Epetra_BlockMap& rowMap,
68  const int* numIndicesPerRow, bool staticProfile)
69  : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
70  CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, staticProfile))
71 {
72  Allocate(numIndicesPerRow, 1, staticProfile);
73 }
74 
75 //==============================================================================
77  const Epetra_BlockMap& rowMap,
78  int numIndicesPerRow, bool staticProfile)
79  : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
80  CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, staticProfile))
81 {
82  Allocate(&numIndicesPerRow, 0, staticProfile);
83 }
84 
85 //==============================================================================
87  const Epetra_BlockMap& rowMap,
88  const Epetra_BlockMap& colMap,
89  const int* numIndicesPerRow, bool staticProfile)
90  : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
91  CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, colMap, staticProfile))
92 {
93  if(!rowMap.GlobalIndicesTypeMatch(colMap))
94  throw ReportError("Epetra_CrsGraph::Epetra_CrsGraph: cannot be called with different indices types for rowMap and colMap", -1);
95 
96  Allocate(numIndicesPerRow, 1, staticProfile);
97 }
98 
99 //==============================================================================
101  const Epetra_BlockMap& rowMap,
102  const Epetra_BlockMap& colMap,
103  int numIndicesPerRow, bool staticProfile)
104  : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
105  CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, colMap, staticProfile))
106 {
107  if(!rowMap.GlobalIndicesTypeMatch(colMap))
108  throw ReportError("Epetra_CrsGraph::Epetra_CrsGraph: cannot be called with different indices types for rowMap and colMap", -1);
109 
110  Allocate(&numIndicesPerRow, 0, staticProfile);
111 }
112 
113 //==============================================================================
115  : Epetra_DistObject(Graph),
116  CrsGraphData_(Graph.CrsGraphData_)
117 {
119 }
120 
121 // private =====================================================================
122 template<typename int_type>
123 int Epetra_CrsGraph::TAllocate(const int* numIndicesPerRow, int Inc, bool staticProfile) {
124  int i;
125  const int numMyBlockRows = CrsGraphData_->NumMyBlockRows_;
126 
127  // Portions specific to ISDVs
128  // Note that NumIndicesPerRow_ will be 1 element longer than needed.
129  // This is because, if OptimizeStorage() is called, the storage associated with
130  // NumIndicesPerRow_ will be reused implicitly for the IndexOffset_ array.
131  // Although a bit fragile, for users who care about efficient memory allocation,
132  // the time and memory fragmentation are important: Mike Heroux Feb 2005.
133  int errorcode = CrsGraphData_->NumIndicesPerRow_.Size(numMyBlockRows+1);
134  if(errorcode != 0)
135  throw ReportError("Error with NumIndicesPerRow_ allocation.", -99);
136 
137  errorcode = CrsGraphData_->NumAllocatedIndicesPerRow_.Size(numMyBlockRows);
138  if(errorcode != 0)
139  throw ReportError("Error with NumAllocatedIndicesPerRow_ allocation.", -99);
140 
141  int nnz = 0;
142  if(CrsGraphData_->CV_ == Copy) {
143  if (numIndicesPerRow != 0) {
144  for(i = 0; i < numMyBlockRows; i++) {
145  int nnzr = numIndicesPerRow[i*Inc];
146  nnz += nnzr;
147  }
148  }
149  }
150  CrsGraphData_->NumMyEntries_ = nnz; // Define this now since we have it and might want to use it
151  //***
152 
154 
156 
157  // Allocate and initialize entries if we are copying data
158  if(CrsGraphData_->CV_ == Copy) {
159  if (staticProfile) Data.All_Indices_.Size(nnz);
160  int_type * all_indices = Data.All_Indices_.Values(); // First address of contiguous buffer
161  for(i = 0; i < numMyBlockRows; i++) {
162  const int NumIndices = numIndicesPerRow==0 ? 0 :numIndicesPerRow[i*Inc];
163  const int_type indexBaseMinusOne = (int_type) IndexBase64() - 1;
164 
165  if(NumIndices > 0) {
166  if (staticProfile) {
167  Data.Indices_[i] = all_indices;
168  all_indices += NumIndices;
169  int_type* ColIndices = Data.Indices_[i];
170  for(int j = 0; j < NumIndices; j++)
171  ColIndices[j] = indexBaseMinusOne; // Fill column indices with out-of-range values
172  }
173  else {
174  // reserve memory in the STL vector, and then resize it to zero
175  // again in order to signal the program that no data is in there
176  // yet.
177  Data.SortedEntries_[i].entries_.resize(NumIndices,
178  indexBaseMinusOne);
179  Data.Indices_[i] = Epetra_Util_data_ptr(Data.SortedEntries_[i].entries_);
180  Data.SortedEntries_[i].entries_.resize(0);
181  }
182  }
183  else {
184  Data.Indices_[i] = NULL;
185  }
186 
187  CrsGraphData_->NumAllocatedIndicesPerRow_[i] = NumIndices;
188  }
189  if (staticProfile) assert(Data.All_Indices_.Values()+nnz==all_indices); // Sanity check
190  }
191  else { // CV_ == View
192  for(i = 0; i < numMyBlockRows; i++) {
193  Data.Indices_[i] = 0;
194  }
195  }
196 
197  SetAllocated(true);
198 
199  return(0);
200 }
201 
202 int Epetra_CrsGraph::Allocate(const int* numIndicesPerRow, int Inc, bool staticProfile)
203 {
204  if(RowMap().GlobalIndicesInt()) {
205  return TAllocate<int>(numIndicesPerRow, Inc, staticProfile);
206  }
207 
208  if(RowMap().GlobalIndicesLongLong()) {
209 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
210  TAllocate<int>(numIndicesPerRow, Inc, staticProfile);
211  TAllocate<long long>(numIndicesPerRow, Inc, staticProfile);
212  return 0;
213 #else
214  throw ReportError("Epetra_CrsGraph::Allocate: ERROR, GlobalIndicesLongLong but no API for it.",-1);
215 #endif
216  }
217 
218  throw ReportError("Epetra_CrsGraph::Allocate: Internal error.", -1);
219 }
220 
221 
222 // private =====================================================================
223 /*
224 int Epetra_CrsGraph::ReAllocate() {
225  // Reallocate storage that was deleted in OptimizeStorage
226 
227  // Since NIPR is in Copy mode, NAIPR will become a Copy as well
228  CrsGraphData_->NumAllocatedIndicesPerRow_ = CrsGraphData_->NumIndicesPerRow_;
229 
230  CrsGraphData_->StorageOptimized_ = false;
231 
232  return(0);
233 }
234 */
235 //==============================================================================
237 {
238  CleanupData();
239 }
240 
241 // private =====================================================================
243  if(CrsGraphData_ != 0) {
245  if(CrsGraphData_->ReferenceCount() == 0) {
246  delete CrsGraphData_;
247  CrsGraphData_ = 0;
248  }
249  }
250 }
251 
252 //==============================================================================
253 template<typename int_type>
254 int Epetra_CrsGraph::InsertGlobalIndices(int_type Row, int NumIndices, int_type* indices) {
255  if(IndicesAreLocal())
256  EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
258  EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
259  SetIndicesAreGlobal(true);
260  int locRow = LRID(Row); // Find local row number for this global row index
261 
262  EPETRA_CHK_ERR(InsertIndicesIntoSorted(locRow, NumIndices, indices));
263 
264  if(CrsGraphData_->ReferenceCount() > 1)
265  return(1);
266  else
267  return(0);
268 }
269 
270 //==============================================================================
271 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
272 int Epetra_CrsGraph::InsertGlobalIndices(int Row, int NumIndices, int* indices) {
273 
274  if(RowMap().GlobalIndicesInt())
275  return InsertGlobalIndices<int>(Row, NumIndices, indices);
276  else
277  throw ReportError("Epetra_CrsGraph::InsertGlobalIndices int version called for a graph that is not int.", -1);
278 }
279 #endif
280 //==============================================================================
281 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
282 int Epetra_CrsGraph::InsertGlobalIndices(long long Row, int NumIndices, long long* indices) {
283 
284  if(RowMap().GlobalIndicesLongLong())
285  return InsertGlobalIndices<long long>(Row, NumIndices, indices);
286  else
287  throw ReportError("Epetra_CrsGraph::InsertGlobalIndices long long version called for a graph that is not long long.", -1);
288 }
289 #endif
290 //==============================================================================
291 int Epetra_CrsGraph::InsertMyIndices(int Row, int NumIndices, int* indices) {
292 
293  if(IndicesAreGlobal()) {
294  EPETRA_CHK_ERR(-2); // Cannot insert local indices into a global graph
295  }
297  EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
298 
299  if (CrsGraphData_->HaveColMap_) {
300  SetIndicesAreLocal(true);
301  }
302  else {
303  if (!IndicesAreLocal()) {
304  EPETRA_CHK_ERR(-4);
305  }
306  }
307 
308 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
309  EPETRA_CHK_ERR(InsertIndicesIntoSorted(Row, NumIndices, indices));
310 #else
311  throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
312 #endif
313 
314  if(CrsGraphData_->ReferenceCount() > 1)
315  return(1);
316  else
317  return(0);
318 }
319 
320 // protected ===================================================================
321 template<typename int_type>
323  int NumIndices,
324  int_type* UserIndices)
325 {
326  if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
327 
328  SetSorted(false); // No longer in sorted state.
329  CrsGraphData_->NoRedundancies_ = false; // Redundancies possible.
330  SetGlobalConstantsComputed(false); // No longer have valid global constants.
331 
332  int j;
333  int ierr = 0;
334 
335  if(Row < 0 || Row >= NumMyBlockRows())
336  EPETRA_CHK_ERR(-2); // Not in Row range
337 
338  int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
339  int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
340 
342 
343  if(CrsGraphData_->CV_ == View) {
344  if(Data.Indices_[Row] != 0)
345  ierr = 2; // This row has been defined already. Issue warning.
346  Data.Indices_[Row] = UserIndices;
347  current_numAllocIndices = NumIndices;
348  current_numIndices = NumIndices;
349  }
350  else {
351  // if HaveColMap_ is true, UserIndices is copied into a new array,
352  // and then modified. The UserIndices pointer is updated to point to this
353  // new array. If HaveColMap_ is false, nothing is done. This way,
354  // the same UserIndices pointer can be used later on regardless of whether
355  // changes were made.
356  if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
357  if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
358  delete [] Data.TempColIndices_;
359  Data.TempColIndices_ = new int_type[NumIndices];
360  CrsGraphData_->NumTempColIndices_ = NumIndices;
361  }
362  int_type * tempIndices = Data.TempColIndices_;
363  int loc = 0;
364  if(IndicesAreLocal()) {
365  for(j = 0; j < NumIndices; ++j)
366  if(CrsGraphData_->ColMap_.MyLID(static_cast<int>(UserIndices[j])))
367  tempIndices[loc++] = UserIndices[j];
368  }
369  else {
370  for(j = 0; j < NumIndices; ++j) {
371  const int Index = CrsGraphData_->ColMap_.LID(UserIndices[j]);
372  if (Index > -1)
373  tempIndices[loc++] = UserIndices[j];
374  }
375 
376  }
377  if(loc != NumIndices)
378  ierr = 2; //Some columns excluded
379  NumIndices = loc;
380  UserIndices = tempIndices;
381  }
382 
383  int start = current_numIndices;
384  int stop = start + NumIndices;
386  if(stop > current_numAllocIndices)
387  EPETRA_CHK_ERR(-2); // Cannot reallocate storage if graph created using StaticProfile
388  }
389  else {
390  if (current_numAllocIndices > 0 && stop > current_numAllocIndices)
391  ierr = 3;
392  Data.SortedEntries_[Row].entries_.resize(stop, IndexBase64() - 1);
393  Data.Indices_[Row] = stop>0 ? &Data.SortedEntries_[Row].entries_[0] : NULL;
394 
395  current_numAllocIndices = (int) Data.SortedEntries_[Row].entries_.capacity();
396  }
397 
398  current_numIndices = stop;
399  int_type* RowIndices = Data.Indices_[Row]+start;
400  for(j = 0; j < NumIndices; j++) {
401  RowIndices[j] = UserIndices[j];
402  }
403  }
404 
405  if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
406  CrsGraphData_->MaxNumIndices_ = current_numIndices;
407  }
408  EPETRA_CHK_ERR(ierr);
409 
410 
411  if(CrsGraphData_->ReferenceCount() > 1)
412  return(1); // return 1 if data is shared
413  else
414  return(0);
415 }
416 
417 // =========================================================================
418 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
420  int NumIndices,
421  int* UserIndices)
422 {
423  if(RowMap().GlobalIndicesTypeValid())
424  return InsertIndices<int>(Row, NumIndices, UserIndices);
425  else
426  throw ReportError("Epetra_CrsGraph::InsertIndices global index type unknown.", -1);
427 }
428 #endif
429 // =========================================================================
430 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
432  int NumIndices,
433  long long* UserIndices)
434 {
435  if(RowMap().GlobalIndicesLongLong())
436  return InsertIndices<long long>(Row, NumIndices, UserIndices);
437  else
438  throw ReportError("Epetra_CrsGraph::InsertIndices long long version called for a graph that is not long long.", -1);
439 }
440 #endif
441 
442 // =========================================================================
443 template<typename int_type>
445  int NumIndices,
446  int_type* UserIndices)
447 {
448  // This function is only valid for COPY mode with non-static profile and
449  // sorted entries. Otherwise, go to the other function.
452  return InsertIndices(Row, NumIndices, UserIndices);
453 
454  if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
455 
456  SetGlobalConstantsComputed(false); // No longer have valid global constants.
457 
458  int ierr = 0;
459 
460  if(Row < 0 || Row >= NumMyBlockRows())
461  EPETRA_CHK_ERR(-2); // Not in Row range
462 
463  int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
464  int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
465 
467 
468  // if HaveColMap_ is true, UserIndices filters out excluded indices,
469  // and then modified. The UserIndices pointer is updated to point to this
470  // new array. If HaveColMap_ is false, nothing is done. This way,
471  // the same UserIndices pointer can be used later on regardless of whether
472  // changes were made.
473  if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
474  if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
475  delete [] Data.TempColIndices_;
476  Data.TempColIndices_ = new int_type[NumIndices];
477  CrsGraphData_->NumTempColIndices_ = NumIndices;
478  }
479  int_type * tempIndices = Data.TempColIndices_;
480  int loc = 0;
481  if(IndicesAreLocal()) {
482  for(int j = 0; j < NumIndices; ++j)
483  if(CrsGraphData_->ColMap_.MyLID(static_cast<int>(UserIndices[j])))
484  tempIndices[loc++] = UserIndices[j];
485  }
486  else {
487  for(int j = 0; j < NumIndices; ++j) {
488  if (CrsGraphData_->ColMap_.MyGID(UserIndices[j])) {
489  tempIndices[loc++] = UserIndices[j];
490  }
491  }
492  }
493  if(loc != NumIndices)
494  ierr = 2; //Some columns excluded
495  NumIndices = loc;
496  UserIndices = tempIndices;
497  }
498 
499  // for non-static profile, directly insert into a list that we always
500  // keep sorted.
501  Data.SortedEntries_[Row].AddEntries(NumIndices, UserIndices);
502  current_numIndices = (int) Data.SortedEntries_[Row].entries_.size();
503  current_numAllocIndices = (int) Data.SortedEntries_[Row].entries_.capacity();
504  // reset the pointer to the respective data
505  Data.Indices_[Row] = current_numIndices > 0 ? &Data.SortedEntries_[Row].entries_[0] : NULL;
506 
507  if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
508  CrsGraphData_->MaxNumIndices_ = current_numIndices;
509  }
510  EPETRA_CHK_ERR(ierr);
511 
512  if(CrsGraphData_->ReferenceCount() > 1)
513  return(1); // return 1 if data is shared
514  else
515  return(0);
516 }
517 
518 //==============================================================================
519 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
521  int NumIndices,
522  int* UserIndices)
523 {
524  if(RowMap().GlobalIndicesTypeValid())
525  return InsertIndicesIntoSorted<int>(Row, NumIndices, UserIndices);
526  else
527  throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted global index type unknown.", -1);
528 }
529 #endif
530 //==============================================================================
531 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
533  int NumIndices,
534  long long* UserIndices)
535 {
536  if(RowMap().GlobalIndicesLongLong())
537  return InsertIndicesIntoSorted<long long>(Row, NumIndices, UserIndices);
538  else
539  throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted long long version called for a graph that is not long long.", -1);
540 }
541 #endif
542 //==============================================================================
543 template<typename int_type>
544 int Epetra_CrsGraph::RemoveGlobalIndices(int_type Row, int NumIndices, int_type* indices) {
545  int j;
546  int k;
547  int ierr = 0;
548  int Loc;
549 
551  EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
552 
553  if(IndicesAreLocal())
554  EPETRA_CHK_ERR(-2); // Cannot remove global indices from a filled graph
555 
556  if(CrsGraphData_->CV_ == View)
557  EPETRA_CHK_ERR(-3); // This is a view only. Cannot remove entries.
558 
559  int locRow = LRID(Row); // Normalize row range
560 
561  if(locRow < 0 || locRow >= NumMyBlockRows())
562  EPETRA_CHK_ERR(-1); // Not in Row range
563 
564  int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[locRow];
565 
567 
568  for(j = 0; j < NumIndices; j++) {
569  int_type Index = indices[j];
570  if(FindGlobalIndexLoc(locRow,Index,j,Loc)) {
571  for(k = Loc+1; k < NumCurrentIndices; k++)
572  Data.Indices_[locRow][k-1] = Data.Indices_[locRow][k];
573  NumCurrentIndices--;
574  CrsGraphData_->NumIndicesPerRow_[locRow]--;
576  Data.SortedEntries_[locRow].entries_.pop_back();
577  else
578  Data.Indices_[locRow][NumCurrentIndices-1] = IndexBase64() - 1;
579  }
580  }
581  SetGlobalConstantsComputed(false); // No longer have valid global constants.
582 
583  EPETRA_CHK_ERR(ierr);
584 
585  if(CrsGraphData_->ReferenceCount() > 1)
586  return(1);
587  else
588  return(0);
589 }
590 
591 //==============================================================================
592 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
593 int Epetra_CrsGraph::RemoveGlobalIndices(int Row, int NumIndices, int* indices)
594 {
595  if(RowMap().GlobalIndicesInt())
596  return RemoveGlobalIndices<int>(Row, NumIndices, indices);
597  else
598  throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices int version called for a graph that is not int.", -1);
599 }
600 #endif
601 //==============================================================================
602 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
603 int Epetra_CrsGraph::RemoveGlobalIndices(long long Row, int NumIndices, long long* indices)
604 {
605  if(RowMap().GlobalIndicesLongLong())
606  return RemoveGlobalIndices<long long>(Row, NumIndices, indices);
607  else
608  throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices long long version called for a graph that is not long long.", -1);
609 }
610 #endif
611 //==============================================================================
612 int Epetra_CrsGraph::RemoveMyIndices(int Row, int NumIndices, int* indices) {
613 
615  EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
616 
617  if(IndicesAreGlobal())
618  EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
619 
620  int j;
621  int k;
622  int ierr = 0;
623  int Loc;
624 
625  if(CrsGraphData_->CV_ == View)
626  EPETRA_CHK_ERR(-3); // This is a view only. Cannot remove entries.
627 
628  if(Row < 0 || Row >= NumMyBlockRows())
629  EPETRA_CHK_ERR(-1); // Not in Row range
630 
631  int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[Row];
632 
634 
635  for(j = 0; j < NumIndices; j++) {
636  int Index = indices[j];
637  if(FindMyIndexLoc(Row,Index,j,Loc)) {
638  for(k = Loc + 1; k < NumCurrentIndices; k++)
639  Data.Indices_[Row][k-1] = Data.Indices_[Row][k];
640  NumCurrentIndices--;
643  Data.SortedEntries_[Row].entries_.pop_back();
644  else
645  Data.Indices_[Row][NumCurrentIndices-1] = (int) IndexBase64() - 1;
646  }
647  }
648  SetGlobalConstantsComputed(false); // No longer have valid global constants.
649 
650  EPETRA_CHK_ERR(ierr);
651 
652  if(CrsGraphData_->ReferenceCount() > 1)
653  return(1);
654  else
655  return(0);
656 }
657 
658 //==============================================================================
659 template<typename int_type>
661  int j;
662  int ierr = 0;
663 
665  EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
666 
667  if(IndicesAreLocal())
668  EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
669  if(CrsGraphData_->CV_ == View)
670  EPETRA_CHK_ERR(-3); // This is a view only. Cannot remove entries.
671 
672  // Normalize row range
673 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
674  int locRow = LRID((int) Row);
675 #else
676  int locRow = LRID(Row);
677 #endif
678 
679  if(locRow < 0 || locRow >= NumMyBlockRows())
680  EPETRA_CHK_ERR(-1); // Not in Row range
681 
683 
685  int NumIndices = CrsGraphData_->NumIndicesPerRow_[locRow];
686 
687  const int_type indexBaseMinusOne = (int_type) IndexBase64() - 1;
688  for(j = 0; j < NumIndices; j++)
689  Data.Indices_[locRow][j] = indexBaseMinusOne; // Set to invalid
690  }
691  else
692  Data.SortedEntries_[locRow].entries_.resize(0);
693 
694  CrsGraphData_->NumIndicesPerRow_[locRow] = 0;
695 
696 
697  SetGlobalConstantsComputed(false); // No longer have valid global constants.
698  EPETRA_CHK_ERR(ierr);
699 
700  if(CrsGraphData_->ReferenceCount() > 1)
701  return(1);
702  else
703  return(0);
704 }
705 
707 {
708  if(RowMap().GlobalIndicesLongLong())
709 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
710  return TRemoveGlobalIndices<long long>(Row);
711 #else
712  throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: ERROR, GlobalIndicesLongLong but no API for it.",-1);
713 #endif
714 
715  if(RowMap().GlobalIndicesInt())
716 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
717  return TRemoveGlobalIndices<int>(Row);
718 #else
719  throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: ERROR, GlobalIndicesInt but no API for it.",-1);
720 #endif
721 
722  throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: Internal error.", -1);
723 }
724 
725 //==============================================================================
727 {
728  int ierr = 0;
729 
731  EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
732 
733  if(IndicesAreGlobal())
734  EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
735 
736  if(CrsGraphData_->CV_ == View)
737  EPETRA_CHK_ERR(-3); // This is a view only. Cannot remove entries.
738 
739  if(Row < 0 || Row >= NumMyBlockRows())
740  EPETRA_CHK_ERR(-1); // Not in Row range
741 
743 
745  int NumIndices = CrsGraphData_->NumIndicesPerRow_[Row];
746  for(int j = 0; j < NumIndices; j++)
747  Data.Indices_[Row][j] = -1; // Set to invalid
748  }
749  else
750  Data.SortedEntries_[Row].entries_.resize(0);
751 
753 
754  SetGlobalConstantsComputed(false); // No longer have valid global constants.
755  EPETRA_CHK_ERR(ierr);
756 
757  if(CrsGraphData_->ReferenceCount() > 1)
758  return(1);
759  else
760  return(0);
761 }
762 
763 // protected ===================================================================
764 template<typename int_type>
766  int_type Index,
767  int Start,
768  int& Loc) const
769 {
770  int NumIndices = NumMyIndices(LocalRow);
771 
772  // If we have transformed the column indices, we must map this global Index to local
774  int* locIndices = Indices(LocalRow);
775  int locIndex = LCID(Index);
776  if (CrsGraphData_->Sorted_) {
777  int insertPoint;
778  Loc = Epetra_Util_binary_search(locIndex, locIndices, NumIndices, insertPoint);
779  return( Loc > -1 );
780  }
781  else {
782  int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
783  for(j = 0; j < NumIndices; j++) {
784  if(j0 >= NumIndices)
785  j0 = 0; // wrap around
786  if(locIndices[j0] == locIndex) {
787  Loc = j0;
788  return(true);
789  }
790  j0++;
791  }
792  }
793  }
794  else {
795  int_type* locIndices = TIndices<int_type>(LocalRow);
796  if (CrsGraphData_->Sorted_) {
797  int insertPoint;
798  Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
799  return( Loc > -1 );
800  }
801  else {
802  int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
803  for(j = 0; j < NumIndices; j++) {
804  if(j0 >= NumIndices)
805  j0 = 0; // wrap around
806  if(locIndices[j0] == Index) {
807  Loc = j0;
808  return(true);
809  }
810  j0++;
811  }
812  }
813  }
814 
815  return(false);
816 }
817 
818 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
820  int Index,
821  int Start,
822  int& Loc) const
823 {
824  if(RowMap().GlobalIndicesInt())
825  return FindGlobalIndexLoc<int>(LocalRow, Index, Start, Loc);
826  else
827  throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc int version called for a graph that is not int.", -1);
828 }
829 #endif
830 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
832  long long Index,
833  int Start,
834  int& Loc) const
835 {
836  if(RowMap().GlobalIndicesLongLong())
837  return FindGlobalIndexLoc<long long>(LocalRow, Index, Start, Loc);
838  else
839  throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc long long version called for a graph that is not long long.", -1);
840 }
841 #endif
842 
843 // protected ===================================================================
844 template<typename int_type>
846  const int_type* indices,
847  int_type Index,
848  int Start,
849  int& Loc) const
850 {
851  // If we have transformed the column indices, we must map this global Index to local
853  Index = LCID(Index);
854  }
855 
856  if (CrsGraphData_->Sorted_) {
857  int insertPoint;
858  Loc = Epetra_Util_binary_search(Index, indices, NumIndices, insertPoint);
859  return( Loc > -1 );
860  }
861  else {
862  int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
863  for(j = 0; j < NumIndices; j++) {
864  if(j0 >= NumIndices)
865  j0 = 0; // wrap around
866  if(indices[j0] == Index) {
867  Loc = j0;
868  return(true);
869  }
870  j0++;
871  }
872  }
873  return(false);
874 }
875 
876 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
878  const int* indices,
879  int Index,
880  int Start,
881  int& Loc) const
882 {
883  if(RowMap().GlobalIndicesInt())
884  return FindGlobalIndexLoc<int>(NumIndices, indices, Index, Start, Loc);
885  else
886  throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc int version called for a graph that is not int.", -1);
887 }
888 #endif
889 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
891  const long long* indices,
892  long long Index,
893  int Start,
894  int& Loc) const
895 {
896  if(RowMap().GlobalIndicesLongLong())
897  return FindGlobalIndexLoc<long long>(NumIndices, indices, Index, Start, Loc);
898  else
899  throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc long long version called for a graph that is not long long.", -1);
900 }
901 #endif
902 
903 // protected ===================================================================
905  int Index,
906  int Start,
907  int& Loc) const
908 {
909  int NumIndices = NumMyIndices(LocalRow);
910  int* locIndices = Indices(LocalRow);
911 
913  throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
914  }
915 
916  if (CrsGraphData_->Sorted_) {
917  int insertPoint;
918  Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
919  return( Loc > -1 );
920  }
921  else {
922  int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
923  for(j = 0; j < NumIndices; j++) {
924  if(j0 >= NumIndices)
925  j0 = 0; // wrap around
926  if(locIndices[j0] == Index) {
927  Loc = j0;
928  return(true);
929  }
930  j0++;
931  }
932  }
933  return(false);
934 }
935 
936 // protected ===================================================================
938  const int* indices,
939  int Index,
940  int Start,
941  int& Loc) const
942 {
944  throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
945  }
946 
947  if (CrsGraphData_->Sorted_) {
948  int insertPoint;
949  Loc = Epetra_Util_binary_search(Index, indices, NumIndices, insertPoint);
950  return( Loc > -1 );
951  }
952  else {
953  int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
954  for(j = 0; j < NumIndices; j++) {
955  if(j0 >= NumIndices)
956  j0 = 0; // wrap around
957  if(indices[j0] == Index) {
958  Loc = j0;
959  return(true);
960  }
961  j0++;
962  }
963  }
964  return(false);
965 }
966 
967 //==============================================================================
970  return(0); // uses FillComplete(...)'s shared-ownership test.
971 }
972 
973 //==============================================================================
974 int Epetra_CrsGraph::FillComplete(const Epetra_BlockMap& domainMap, const Epetra_BlockMap& rangeMap) {
975  if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
976  throw ReportError("Epetra_CrsGraph::FillComplete: cannot be called with different indices types for domainMap and rangeMap", -1);
977 
978  if(!RowMap().GlobalIndicesTypeMatch(domainMap))
979  throw ReportError("Epetra_CrsGraph::FillComplete: cannot be called with different indices types for row map and incoming rangeMap", -1);
980 
981  CrsGraphData_->DomainMap_ = domainMap;
982  CrsGraphData_->RangeMap_ = rangeMap;
983 
984  MakeIndicesLocal(domainMap, rangeMap); // Convert global indices to local indices to on each processor
985  SortIndices(); // Sort column entries from smallest to largest
986  RemoveRedundantIndices(); // Get rid of any redundant index values
987  DetermineTriangular(); //determine if lower/upper triangular
988  CrsGraphData_->MakeImportExport(); // Build Import or Export objects
989  ComputeGlobalConstants(); // Compute constants that require communication
990  SetFilled(true);
991 
992  if(CrsGraphData_->ReferenceCount() > 1)
993  return(1);
994  else
995  return(0);
996 }
997 
998 //==============================================================================
1000  return(FillComplete());
1001 }
1002 
1003 //==============================================================================
1005  return(FillComplete(*domainMap, *rangeMap));
1006 }
1007 
1008 // private =====================================================================
1010 {
1012  return(0);
1013 
1014 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1016 #else
1018 #endif
1019  tempvec(8); // Temp space
1020 
1021  const int numMyBlockRows = NumMyBlockRows();
1022 
1023  CrsGraphData_->NumMyEntries_ = 0; // Compute Number of Nonzero entries and max
1025  {for(int i = 0; i < numMyBlockRows; i++) {
1028  }}
1029 
1030  // Case 1: Constant block size (including blocksize = 1)
1031  if(RowMap().ConstantElementSize() && ColMap().ConstantElementSize() && RowMap().ElementSize() == ColMap().ElementSize()) {
1032  // Jim Westfall reported a fix on 22 June 2013 where the second two conditions
1033  // above are necessary. The added conditions check for the case when the row map
1034  // and col map are constant but different as possible with VBR sub matrices used
1035  // in global assemble
1036 
1037  tempvec[0] = CrsGraphData_->NumMyEntries_;
1038  tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
1039 
1040  Comm().SumAll(&tempvec[0], &tempvec[2], 2);
1041  int tmp_MaxNumIndices = CrsGraphData_->MaxNumIndices_;
1042  Comm().MaxAll(&tmp_MaxNumIndices, &CrsGraphData_->GlobalMaxNumIndices_, 1);
1043 
1044  CrsGraphData_->NumGlobalEntries_ = tempvec[2];
1046 
1047  int RowElementSize = RowMap().MaxElementSize();
1048  int ColElementSize = RowElementSize;
1049  CrsGraphData_->NumGlobalDiagonals_ = tempvec[3] * RowElementSize;
1050  CrsGraphData_->NumMyNonzeros_ = CrsGraphData_->NumMyEntries_ * RowElementSize * ColElementSize;
1051  CrsGraphData_->NumGlobalNonzeros_ = CrsGraphData_->NumGlobalEntries_ * RowElementSize * ColElementSize;
1052  CrsGraphData_->MaxNumNonzeros_ = CrsGraphData_->MaxNumIndices_ * RowElementSize * ColElementSize;
1053  CrsGraphData_->GlobalMaxNumNonzeros_ = CrsGraphData_->GlobalMaxNumIndices_ * RowElementSize * ColElementSize;
1054  }
1055 
1056  // Case 2: Variable block size (more work)
1057  else {
1058  CrsGraphData_->NumMyNonzeros_ = 0; // We will count the number of nonzeros here
1059  CrsGraphData_->MaxNumNonzeros_ = 0; // We will determine the max number of nonzeros in any one block row
1060  int* RowElementSizeList = RowMap().ElementSizeList();
1061  int* ColElementSizeList = RowElementSizeList;
1062  if(Importer() != 0) ColElementSizeList = ColMap().ElementSizeList();
1063  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1064  for(int i = 0; i < numMyBlockRows; i++){
1065  int NumEntries = CrsGraphData_->NumIndicesPerRow_[i];
1066  int* indices = intData.Indices_[i];
1067  if(NumEntries > 0) {
1068  int CurNumNonzeros = 0;
1069  int RowDim = RowElementSizeList[i];
1070  for(int j = 0; j < NumEntries; j++) {
1071  int ColDim = ColElementSizeList[indices[j]];
1072  CurNumNonzeros += RowDim*ColDim;
1074  }
1076  CrsGraphData_->NumMyNonzeros_ += CurNumNonzeros;
1077  }
1078  }
1079 
1080  // Sum Up all nonzeros
1081 
1082  tempvec[0] = CrsGraphData_->NumMyEntries_;
1083  tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
1084  tempvec[2] = CrsGraphData_->NumMyDiagonals_;
1085  tempvec[3] = CrsGraphData_->NumMyNonzeros_;
1086 
1087  Comm().SumAll(&tempvec[0], &tempvec[4], 4);
1088 
1089  CrsGraphData_->NumGlobalEntries_ = tempvec[4];
1091  CrsGraphData_->NumGlobalDiagonals_ = tempvec[6];
1092  CrsGraphData_->NumGlobalNonzeros_ = tempvec[7];
1093 
1094  tempvec[0] = CrsGraphData_->MaxNumIndices_;
1095  tempvec[1] = CrsGraphData_->MaxNumNonzeros_;
1096 
1097  Comm().MaxAll(&tempvec[0], &tempvec[2], 2);
1098 
1099  CrsGraphData_->GlobalMaxNumIndices_ = (int) tempvec[2];
1100  CrsGraphData_->GlobalMaxNumNonzeros_ = (int) tempvec[3];
1101  }
1102 
1105 
1107 
1108  return(0);
1109 }
1110 
1111 void epetra_shellsort(int* list, int length)
1112 {
1113  int i, j, j2, temp, istep;
1114  unsigned step;
1115 
1116  step = length/2;
1117  while (step > 0)
1118  {
1119  for (i=step; i < length; i++)
1120  {
1121  istep = step;
1122  j = i;
1123  j2 = j-istep;
1124  temp = list[i];
1125  if (list[j2] > temp) {
1126  while ((j >= istep) && (list[j2] > temp))
1127  {
1128  list[j] = list[j2];
1129  j = j2;
1130  j2 -= istep;
1131  }
1132  list[j] = temp;
1133  }
1134  }
1135 
1136  if (step == 2)
1137  step = 1;
1138  else
1139  step = (int) (step / 2.2);
1140  }
1141 }
1142 
1143 //==============================================================================
1145  if(IndicesAreGlobal())
1146  EPETRA_CHK_ERR(-1);
1147  if(Sorted())
1148  return(0);
1149 
1150  // For each row, sort column entries from smallest to largest.
1151  // Use shell sort, which is fast if indices are already sorted.
1152 
1153  const int numMyBlockRows = NumMyBlockRows();
1154  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1155  for(int i = 0; i < numMyBlockRows; i++){
1156  int n = CrsGraphData_->NumIndicesPerRow_[i];
1157  int* const list = intData.Indices_[i];
1158 
1159  epetra_shellsort(list, n);
1160 // int m = n/2;
1161 
1162 // while(m > 0) {
1163  // int max = n - m;
1164 // for(int j = 0; j < max; j++) {
1165 // int k = j;
1166 // while(k>-1) {
1167 // if(list[k+m] >= list[k])
1168 // break;
1169 // int itemp = list[k+m];
1170 // list[k+m] = list[k];
1171 // list[k] = itemp;
1172 // k-=m;
1173 // }
1174 // }
1175 // m = m/2;
1176 // }
1177  }
1178  SetSorted(true);
1179 
1180  if(CrsGraphData_->ReferenceCount() > 1)
1181  return(1);
1182  else
1183  return(0);
1184 }
1185 
1186 void epetra_crsgraph_compress_out_duplicates(int len, int* list, int& newlen)
1187 {
1188  //
1189  //This function runs the array ('list') checking for
1190  //duplicate entries. Any duplicates that are found are
1191  //removed by sliding subsequent data down in the array,
1192  //over-writing the duplicates. Finally, the new length
1193  //of the array (i.e., the number of unique entries) is
1194  //placed in the output parameter 'newlen'. The array is
1195  //**not** re-allocated.
1196  //
1198  //Important assumption: The array contents are assumed to
1199  //be sorted before this function is called. If the array
1200  //contents are not sorted, then the behavior of this
1201  //function is undefined.
1203  //
1204 
1205  if (len < 2) return;
1206 
1207  int* ptr0 = &list[0];
1208  int* ptr1 = &list[1];
1209 
1210  int* ptr_end = &list[len-1];
1211 
1212  while(*ptr0 != *ptr1 && ptr1 < ptr_end) {
1213  ++ptr0;
1214  ++ptr1;
1215  }
1216 
1217  if (ptr1 < ptr_end) {
1218  //if ptr1 < ptr_end we've found a duplicate...
1219 
1220  ++ptr0;
1221  ++ptr1;
1222 
1223  while(*ptr0 == *ptr1 && ptr1 < ptr_end) ++ptr1;
1224 
1225  while(ptr1 < ptr_end) {
1226 
1227  int val = *ptr1++;
1228 
1229  while(val == *ptr1 && ptr1 < ptr_end) {
1230  ++ptr1;
1231  }
1232 
1233  *ptr0++ = val;
1234  }
1235 
1236  if (*(ptr0-1) != *ptr1) *ptr0++ = *ptr1;
1237 
1238  int num_removed = (int)(ptr_end - ptr0 + 1);
1239  newlen = len - num_removed;
1240  }
1241  else {
1242  if (*ptr0 == *ptr1) newlen = len - 1;
1243  else newlen = len;
1244  }
1245 }
1246 
1247 //==============================================================================
1249 {
1250  if(!Sorted())
1251  EPETRA_CHK_ERR(-1); // Must have sorted index set
1252  if(IndicesAreGlobal())
1253  EPETRA_CHK_ERR(-2); // Indices must be local
1254 
1255  // Note: This function assumes that SortIndices was already called.
1256  // For each row, remove column indices that are repeated.
1257 
1258  const int numMyBlockRows = NumMyBlockRows();
1259  bool found_redundancies = false;
1260 
1261  if(NoRedundancies() == false) {
1262  int* numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
1264  for(int i=0; i<numMyBlockRows; ++i) {
1265  int NumIndices = numIndicesPerRow[i];
1266  int* col_indices = this->Indices(i);
1267 
1268  if(NumIndices > 1) {
1269  epetra_crsgraph_compress_out_duplicates(NumIndices, col_indices,
1270  numIndicesPerRow[i]);
1271  }
1272  if (NumIndices != numIndicesPerRow[i]) {
1273  found_redundancies = true;
1274  }
1275  }
1276  if (found_redundancies && !CrsGraphData_->StaticProfile_)
1277  {
1278  for(int i=0; i<numMyBlockRows; ++i) {
1279  int* col_indices = this->Indices(i);
1280 
1281  // update vector size and address in memory
1282  intData.SortedEntries_[i].entries_.assign(col_indices, col_indices+numIndicesPerRow[i]);
1283  if (numIndicesPerRow[i] > 0) {
1284  intData.Indices_[i] = &(intData.SortedEntries_[i].entries_[0]);
1285  }
1286  else {
1287  intData.Indices_[i] = NULL;
1288  }
1289  }
1290  }
1291  }
1292 
1293  SetNoRedundancies(true);
1294  return 0;
1295 }
1296 
1297 //==============================================================================
1299 {
1300  // determine if graph is upper or lower triangular or has no diagonal
1301 
1302  if(!Sorted())
1303  EPETRA_CHK_ERR(-1); // Must have sorted index set
1304  if(IndicesAreGlobal())
1305  EPETRA_CHK_ERR(-2); // Indices must be local
1306 
1309 
1310  const Epetra_BlockMap& rowMap = RowMap();
1311  const Epetra_BlockMap& colMap = ColMap();
1312 
1313  const int numMyBlockRows = NumMyBlockRows();
1314 
1315  for(int i = 0; i < numMyBlockRows; i++) {
1316  int NumIndices = NumMyIndices(i);
1317  if(NumIndices > 0) {
1318 #if defined(EPETRA_NO_64BIT_GLOBAL_INDICES) && !defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
1319  int ig = rowMap.GID(i);
1320 #else
1321  long long ig = rowMap.GID64(i);
1322 #endif
1323  int* col_indices = this->Indices(i);
1324 
1325  int jl_0 = col_indices[0];
1326  int jl_n = col_indices[NumIndices-1];
1327 
1328  if(jl_n > i) CrsGraphData_->LowerTriangular_ = false;
1329  if(jl_0 < i) CrsGraphData_->UpperTriangular_ = false;
1330 
1331  //jl will be the local-index for the diagonal that we
1332  //want to search for.
1333  int jl = colMap.LID(ig);
1334 
1335  int insertPoint = -1;
1336  if (Epetra_Util_binary_search(jl, col_indices, NumIndices, insertPoint)>-1) {
1339  }
1340  }
1341  }
1342 
1344 
1345  if(CrsGraphData_->ReferenceCount() > 1)
1346  return(1);
1347  else
1348  return(0);
1349 }
1350 
1351 // private =====================================================================
1352 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1354  const Epetra_BlockMap& rangeMap)
1355 {
1356  (void)rangeMap;
1357  int i;
1358  int j;
1359 
1361  return(0); // Already have a Column Map
1362 
1363  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1364  if(IndicesAreLocal())
1365  EPETRA_CHK_ERR(-1); // Return error: Indices must be global
1366 
1367  // Scan all column indices and sort into two groups:
1368  // Local: those whose GID matches a GID of the domain map on this processor and
1369  // Remote: All others.
1370  int numDomainElements = domainMap.NumMyElements();
1371  bool * LocalGIDs = 0;
1372  if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
1373  for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
1374 
1375  // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
1376  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
1377  // and the number of block rows.
1378  const int numMyBlockRows = NumMyBlockRows();
1379  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
1380  //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
1381  Epetra_HashTable<int> RemoteGIDs(hashsize);
1382  Epetra_HashTable<int> RemoteGIDList(hashsize);
1383 
1384  int NumLocalColGIDs = 0;
1385  int NumRemoteColGIDs = 0;
1386  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1387  for(i = 0; i < numMyBlockRows; i++) {
1388  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1389  int* ColIndices = intData.Indices_[i];
1390  for(j = 0; j < NumIndices; j++) {
1391  int GID = ColIndices[j];
1392  // Check if GID matches a row GID
1393  int LID = domainMap.LID(GID);
1394  if(LID != -1) {
1395  bool alreadyFound = LocalGIDs[LID];
1396  if (!alreadyFound) {
1397  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1398  NumLocalColGIDs++;
1399  }
1400  }
1401  else {
1402  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1403  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1404  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1405  }
1406  }
1407  }
1408  }
1409 
1410  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
1411  if (domainMap.Comm().NumProc()==1) {
1412 
1413  if (NumRemoteColGIDs!=0) {
1414  throw ReportError("Some column IDs are not in domainMap. If matrix is rectangular, you must pass in domainMap to FillComplete",-1); // Sanity test: When one processor,there can be no remoteGIDs
1415  }
1416  if (NumLocalColGIDs==numDomainElements) {
1417  CrsGraphData_->ColMap_ = domainMap;
1418  CrsGraphData_->HaveColMap_ = true;
1419  if (LocalGIDs!=0) delete [] LocalGIDs;
1420  return(0);
1421  }
1422  }
1423 
1424  // Now build integer array containing column GIDs
1425  // Build back end, containing remote GIDs, first
1426  int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1427  Epetra_IntSerialDenseVector ColIndices;
1428  if(numMyBlockCols > 0)
1429  ColIndices.Size(numMyBlockCols);
1430 
1431  int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
1432 
1433  for(i = 0; i < NumRemoteColGIDs; i++)
1434  RemoteColIndices[i] = RemoteGIDList.Get(i);
1435 
1436  int NLists = 1;
1438  Epetra_IntSerialDenseVector SizeList;
1439  int* RemoteSizeList = 0;
1440  bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
1441 
1442  if(NumRemoteColGIDs > 0)
1443  PIDList.Size(NumRemoteColGIDs);
1444 
1445  if(DoSizes) {
1446  if(numMyBlockCols > 0)
1447  SizeList.Size(numMyBlockCols);
1448  RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
1449  NLists++;
1450  }
1451  EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
1452 
1453  // Sort External column indices so that all columns coming from a given remote processor are contiguous
1454 
1455  Epetra_Util Util;
1456  int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
1457  SortLists[0] = RemoteColIndices;
1458  SortLists[1] = RemoteSizeList;
1459  Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists, 0, 0);
1461  // Sort external column indices so that columns from a given remote processor are not only contiguous
1462  // but also in ascending order. NOTE: I don't know if the number of externals associated
1463  // with a given remote processor is known at this point ... so I count them here.
1464 
1465  NLists--;
1466  int StartCurrent, StartNext;
1467  StartCurrent = 0; StartNext = 1;
1468  while ( StartNext < NumRemoteColGIDs ) {
1469  if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
1470  else {
1471  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1472  Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
1473  StartCurrent = StartNext; StartNext++;
1474  }
1475  }
1476  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1477  Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
1478  }
1479 
1480  // Now fill front end. Two cases:
1481  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
1482  // can simply read the domain GIDs into the front part of ColIndices, otherwise
1483  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
1484  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
1485 
1486  if(NumLocalColGIDs == domainMap.NumMyElements()) {
1487  domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
1488  if(DoSizes)
1489  domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
1490  }
1491  else {
1492  int NumMyElements = domainMap.NumMyElements();
1493  int* MyGlobalElements = domainMap.MyGlobalElements();
1494  int* ElementSizeList = 0;
1495  if(DoSizes)
1496  ElementSizeList = domainMap.ElementSizeList();
1497  int NumLocalAgain = 0;
1498  for(i = 0; i < NumMyElements; i++) {
1499  if(LocalGIDs[i]) {
1500  if(DoSizes)
1501  SizeList[NumLocalAgain] = ElementSizeList[i];
1502  ColIndices[NumLocalAgain++] = MyGlobalElements[i];
1503  }
1504  }
1505  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
1506  }
1507 
1508  // Done with this array
1509  if (LocalGIDs!=0) delete [] LocalGIDs;
1510 
1511 
1512  // Make Column map with same element sizes as Domain map
1513 
1514  if(domainMap.MaxElementSize() == 1) { // Simple map
1515  Epetra_Map temp((int) -1, numMyBlockCols, ColIndices.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
1516  CrsGraphData_->ColMap_ = temp;
1517  }
1518  else if(domainMap.ConstantElementSize()) { // Constant Block size map
1519  Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(), (int) domainMap.IndexBase64(), domainMap.Comm());
1520  CrsGraphData_->ColMap_ = temp;
1521  }
1522  else { // Most general case where block size is variable.
1523  Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
1524  CrsGraphData_->ColMap_ = temp;
1525  }
1526  CrsGraphData_->HaveColMap_ = true;
1527 
1528  return(0);
1529 }
1530 #endif
1531 //==============================================================================
1532 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1534  const Epetra_BlockMap& rangeMap)
1535 {
1536  (void)rangeMap;
1537  int i;
1538  int j;
1539 
1541  return(0); // Already have a Column Map
1542 
1543  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1544  if(IndicesAreLocal())
1545  EPETRA_CHK_ERR(-1); // Return error: Indices must be global
1546 
1547  // Scan all column indices and sort into two groups:
1548  // Local: those whose GID matches a GID of the domain map on this processor and
1549  // Remote: All others.
1550  int numDomainElements = domainMap.NumMyElements();
1551  bool * LocalGIDs = 0;
1552  if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
1553  for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
1554 
1555  // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
1556  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
1557  // and the number of block rows.
1558  const int numMyBlockRows = NumMyBlockRows();
1559  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
1560  //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
1561  Epetra_HashTable<int> RemoteGIDs(hashsize);
1562  Epetra_HashTable<long long> RemoteGIDList(hashsize);
1563 
1564  int NumLocalColGIDs = 0;
1565  int NumRemoteColGIDs = 0;
1566 
1567  if(IndicesAreLocal())
1568  {
1570 
1571  for(i = 0; i < numMyBlockRows; i++) {
1572  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1573  int* ColIndices = intData.Indices_[i];
1574  for(j = 0; j < NumIndices; j++) {
1575  int GID = ColIndices[j];
1576  // Check if GID matches a row GID
1577  int LID = domainMap.LID(GID);
1578  if(LID != -1) {
1579  bool alreadyFound = LocalGIDs[LID];
1580  if (!alreadyFound) {
1581  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1582  NumLocalColGIDs++;
1583  }
1584  }
1585  else {
1586  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1587  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1588  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1589  }
1590  }
1591  }
1592  }
1593  }
1594  else if(IndicesAreGlobal())
1595  {
1597 
1598  for(i = 0; i < numMyBlockRows; i++) {
1599  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1600  long long* ColIndices = LLData.Indices_[i];
1601  for(j = 0; j < NumIndices; j++) {
1602  long long GID = ColIndices[j];
1603  // Check if GID matches a row GID
1604  int LID = domainMap.LID(GID);
1605  if(LID != -1) {
1606  bool alreadyFound = LocalGIDs[LID];
1607  if (!alreadyFound) {
1608  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1609  NumLocalColGIDs++;
1610  }
1611  }
1612  else {
1613  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1614  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1615  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1616  }
1617  }
1618  }
1619  }
1620  }
1621 
1622  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
1623  if (domainMap.Comm().NumProc()==1) {
1624 
1625  if (NumRemoteColGIDs!=0) {
1626  throw ReportError("Some column IDs are not in domainMap. If matrix is rectangular, you must pass in domainMap to FillComplete",-1); // Sanity test: When one processor,there can be no remoteGIDs
1627  }
1628  if (NumLocalColGIDs==numDomainElements) {
1629  CrsGraphData_->ColMap_ = domainMap;
1630  CrsGraphData_->HaveColMap_ = true;
1631  if (LocalGIDs!=0) delete [] LocalGIDs;
1632  return(0);
1633  }
1634  }
1635 
1636  // Now build integer array containing column GIDs
1637  // Build back end, containing remote GIDs, first
1638  int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1640  if(numMyBlockCols > 0)
1641  ColIndices.Size(numMyBlockCols);
1642 
1643  long long* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
1644 
1645  for(i = 0; i < NumRemoteColGIDs; i++)
1646  RemoteColIndices[i] = RemoteGIDList.Get(i);
1647 
1648  int NLists = 0;
1650  Epetra_IntSerialDenseVector SizeList;
1651  int* RemoteSizeList = 0;
1652  bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
1653 
1654  if(NumRemoteColGIDs > 0)
1655  PIDList.Size(NumRemoteColGIDs);
1656 
1657  if(DoSizes) {
1658  if(numMyBlockCols > 0)
1659  SizeList.Size(numMyBlockCols);
1660  RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
1661  NLists++;
1662  }
1663  EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
1664 
1665  // Sort External column indices so that all columns coming from a given remote processor are contiguous
1666 
1667  Epetra_Util Util;
1668  //int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
1669  //SortLists[0] = RemoteColIndices;
1670  //SortLists[1] = RemoteSizeList;
1671  Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, &RemoteSizeList, 1, &RemoteColIndices);
1673  // Sort external column indices so that columns from a given remote processor are not only contiguous
1674  // but also in ascending order. NOTE: I don't know if the number of externals associated
1675  // with a given remote processor is known at this point ... so I count them here.
1676 
1677  int* SortLists[1] = {0};
1678 
1679  int StartCurrent, StartNext;
1680  StartCurrent = 0; StartNext = 1;
1681  while ( StartNext < NumRemoteColGIDs ) {
1682  if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
1683  else {
1684  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1685  Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
1686  StartCurrent = StartNext; StartNext++;
1687  }
1688  }
1689  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1690  Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
1691  }
1692 
1693  // Now fill front end. Two cases:
1694  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
1695  // can simply read the domain GIDs into the front part of ColIndices, otherwise
1696  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
1697  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
1698 
1699  if(NumLocalColGIDs == domainMap.NumMyElements()) {
1700  domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
1701  if(DoSizes)
1702  domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
1703  }
1704  else {
1705  int NumMyElements = domainMap.NumMyElements();
1706  long long* MyGlobalElements = domainMap.MyGlobalElements64();
1707  int* ElementSizeList = 0;
1708  if(DoSizes)
1709  ElementSizeList = domainMap.ElementSizeList();
1710  int NumLocalAgain = 0;
1711  for(i = 0; i < NumMyElements; i++) {
1712  if(LocalGIDs[i]) {
1713  if(DoSizes)
1714  SizeList[NumLocalAgain] = ElementSizeList[i];
1715  ColIndices[NumLocalAgain++] = MyGlobalElements[i];
1716  }
1717  }
1718  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
1719  }
1720 
1721  // Done with this array
1722  if (LocalGIDs!=0) delete [] LocalGIDs;
1723 
1724 
1725  // Make Column map with same element sizes as Domain map
1726 
1727  if(domainMap.MaxElementSize() == 1) { // Simple map
1728  Epetra_Map temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.IndexBase64(), domainMap.Comm());
1729  CrsGraphData_->ColMap_ = temp;
1730  }
1731  else if(domainMap.ConstantElementSize()) { // Constant Block size map
1732  Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(),domainMap.IndexBase64(), domainMap.Comm());
1733  CrsGraphData_->ColMap_ = temp;
1734  }
1735  else { // Most general case where block size is variable.
1736  Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), domainMap.IndexBase64(), domainMap.Comm());
1737  CrsGraphData_->ColMap_ = temp;
1738  }
1739  CrsGraphData_->HaveColMap_ = true;
1740 
1741  return(0);
1742 }
1743 #endif
1744 
1746  const Epetra_BlockMap& rangeMap)
1747 {
1748  if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
1749  throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for domainMap and rangeMap", -1);
1750 
1751  if(!RowMap().GlobalIndicesTypeMatch(domainMap))
1752  throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for row map and incoming rangeMap", -1);
1753 
1754  if(RowMap().GlobalIndicesInt())
1755 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1756  return MakeColMap_int(domainMap, rangeMap);
1757 #else
1758  throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesInt but no API for it.",-1);
1759 #endif
1760 
1761  if(RowMap().GlobalIndicesLongLong())
1762 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1763  return MakeColMap_LL(domainMap, rangeMap);
1764 #else
1765  throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesLongLong but no API for it.",-1);
1766 #endif
1767 
1768  throw ReportError("Epetra_CrsGraph::MakeColMap: Internal error, unable to determine global index type of maps", -1);
1769 }
1770 
1771 // protected ===================================================================
1773  if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
1774  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for domainMap and rangeMap", -1);
1775 
1776  if(!RowMap().GlobalIndicesTypeMatch(domainMap))
1777  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for row map and incoming rangeMap", -1);
1778 
1779  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1781  EPETRA_CHK_ERR(-1); // Return error: Indices must not be both local and global
1782 
1783  MakeColMap(domainMap, rangeMap); // If user has not prescribed column map, create one from indices
1784  const Epetra_BlockMap& colmap = ColMap();
1785 
1786  // Store number of local columns
1789 
1790  // Transform indices to local index space
1791  const int numMyBlockRows = NumMyBlockRows();
1792 
1793  if(IndicesAreGlobal()) {
1794  // Check if ColMap is monotone. If not, the list will get unsorted.
1795  bool mapMonotone = true;
1796  {
1797  long long oldGID = colmap.GID64(0);
1798  for (int i=1; i<colmap.NumMyElements(); ++i) {
1799  if (oldGID > colmap.GID64(i)) {
1800  mapMonotone = false;
1801  break;
1802  }
1803  oldGID = colmap.GID64(i);
1804  }
1805  }
1806  if (Sorted())
1807  SetSorted(mapMonotone);
1808 
1809  // We don't call Data<int>() here because that would not work (it will throw)
1810  // if GlobalIndicesLongLong() and IndicesAreGlobal(). This is because
1811  // the local flag is not set yet. We are in the middle of the transaction here.
1812  // In all other cases, one should call the function Data<int> or Data<long long>
1813  // instead of obtaining the pointer directly.
1815 
1816  if(RowMap().GlobalIndicesInt())
1817  {
1818  // now comes the actual transformation
1819  for(int i = 0; i < numMyBlockRows; i++) {
1820  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1821  int* ColIndices = intData.Indices_[i];
1822  for(int j = 0; j < NumIndices; j++) {
1823  int GID = ColIndices[j];
1824  int LID = colmap.LID(GID);
1825  if(LID != -1)
1826  ColIndices[j] = LID;
1827  else
1828  throw ReportError("Internal error in FillComplete ",-1);
1829  }
1830  }
1831  }
1832  else if(RowMap().GlobalIndicesLongLong())
1833  {
1834 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1836 
1837  if (!CrsGraphData_->StaticProfile_) {
1838  // Use the resize trick used in TAllocate.
1839  const long long indexBaseMinusOne = IndexBase64() - 1;
1840  for(int i = 0; i < numMyBlockRows; i++) {
1841  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1842  intData.SortedEntries_[i].entries_.resize(NumIndices, indexBaseMinusOne);
1843  intData.Indices_[i] = NumIndices > 0 ? &intData.SortedEntries_[i].entries_[0]: NULL;
1844  intData.SortedEntries_[i].entries_.resize(0);
1845  }
1846  }
1847 
1848  // now comes the actual transformation
1849  for(int i = 0; i < numMyBlockRows; i++) {
1850  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1851  long long* ColIndices = LL_Data.Indices_[i];
1852  int* intColIndices = intData.Indices_[i];
1853  for(int j = 0; j < NumIndices; j++) {
1854  long long GID = ColIndices[j];
1855  int LID = colmap.LID(GID);
1856  if(LID != -1)
1857  intColIndices[j] = LID;
1858  else
1859  throw ReportError("Internal error in FillComplete ",-1);
1860  }
1861  }
1862 
1863  LL_Data.Deallocate(); // deallocate long long data since indices are local now.
1864 #else
1865  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: GlobalIndicesLongLong but no long long API", -1);
1866 #endif
1867  }
1868 
1869  }
1870 
1871  SetIndicesAreLocal(true);
1872  SetIndicesAreGlobal(false);
1873 
1874  if(CrsGraphData_->ReferenceCount() > 1)
1875  return(1); // return 1 if data is shared
1876  else
1877  return(0);
1878 }
1879 
1880 //==============================================================================
1882  int NumIndices;
1883  const int numMyBlockRows = NumMyBlockRows();
1884 
1886 
1887  if(StorageOptimized())
1888  return(0); // Have we been here before?
1889  if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
1890 
1891  bool Contiguous = true; // Assume contiguous is true
1892  for(int i = 1; i < numMyBlockRows; i++) {
1893  NumIndices = CrsGraphData_->NumIndicesPerRow_[i-1];
1894  int NumAllocateIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[i-1];
1895 
1896  // Check if NumIndices is same as NumAllocatedIndices and
1897  // check if end of beginning of current row starts immediately after end of previous row.
1898  if((NumIndices != NumAllocateIndices) ||
1899  (Data.Indices_[i] != Data.Indices_[i-1] + NumIndices)) {
1900  Contiguous = false;
1901  break;
1902  }
1903  }
1904 
1905  // NOTE: At the end of the above loop set, there is a possibility that NumIndices and NumAllocatedIndices
1906  // for the last row could be different, but I don't think it matters.
1907 
1908 
1909  if((CrsGraphData_->CV_ == View) && !Contiguous)
1910  return(3); // This is user data, it's not contiguous and we can't make it so.
1911 
1912  // This next step constructs the scan sum of the number of indices per row. Note that the
1913  // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
1914  // careful with this sum operation
1915 
1918 
1919  int * numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
1920  int curNumIndices = numIndicesPerRow[0];
1921  numIndicesPerRow[0] = 0;
1922  for (int i=0; i<numMyBlockRows; ++i) {
1923  int nextNumIndices = numIndicesPerRow[i+1];
1924  numIndicesPerRow[i+1] = numIndicesPerRow[i]+curNumIndices;
1925  curNumIndices = nextNumIndices;
1926  }
1927 
1928 // *******************************
1929 // Data NOT contiguous, make it so
1930 // *******************************
1931  if(!Contiguous) { // Must pack indices if not already contiguous
1932 
1933  // Allocate one big integer array for all index values
1934  if (!(StaticProfile())) { // If static profile, All_Indices_ is already allocated, only need to pack data
1935  int errorcode = Data.All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
1936  if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
1937  }
1938  // Pack indices into All_Indices_
1939 
1940  int* all_indices = Data.All_Indices_.Values();
1941  int * indexOffset = CrsGraphData_->IndexOffset_.Values();
1942  int ** indices = Data.Indices_;
1943 
1944  if (!(StaticProfile())) {
1945 #ifdef EPETRA_HAVE_OMP
1946 #pragma omp parallel for default(none) shared(indexOffset,all_indices,indices)
1947 #endif
1948  for(int i = 0; i < numMyBlockRows; i++) {
1949  int numColIndices = indexOffset[i+1] - indexOffset[i];
1950  int* ColIndices = indices[i];
1951  int *newColIndices = all_indices+indexOffset[i];
1952  for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
1953  }
1954  for(int i = 0; i < numMyBlockRows; i++) {
1955  if (indices[i]!=0) {
1956  Data.SortedEntries_[i].entries_.clear();
1957  indices[i] = 0;
1958  }
1959  }
1960  } // End of non-contiguous non-static section
1961  else {
1962 
1963  for(int i = 0; i < numMyBlockRows; i++) {
1964  int numColIndices = indexOffset[i+1] - indexOffset[i];
1965  int* ColIndices = indices[i];
1966  int *newColIndices = all_indices+indexOffset[i];
1967  if (ColIndices!=newColIndices) // No need to copy if pointing to same space
1968  for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
1969  indices[i] = 0;
1970  }
1971  } // End of non-Contiguous static section
1972  } // End of non-Contiguous section
1973  else { // Start of Contiguous section
1974  // if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
1975  // Execute the assignment block in parallel using the same pattern as SpMV
1976  // in order to improve page placement
1977  if (numMyBlockRows > 0 && !(StaticProfile())) {
1978  const int numMyNonzeros = NumMyNonzeros();
1979  int errorcode = Data.All_Indices_.Size(numMyNonzeros);
1980  if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
1981  int* new_all_indices = Data.All_Indices_.Values();
1982  int* old_all_indices = Data.Indices_[0];
1983  int * indexOffset = CrsGraphData_->IndexOffset_.Values();
1984 
1985 #ifdef EPETRA_HAVE_OMP
1986 #pragma omp parallel for default(none) shared(indexOffset,old_all_indices,new_all_indices)
1987 #endif
1988  for(int i = 0; i < numMyBlockRows; i++) {
1989  int numColIndices = indexOffset[i+1] - indexOffset[i];
1990  int *oldColIndices = old_all_indices+indexOffset[i];
1991  int *newColIndices = new_all_indices+indexOffset[i];
1992  for(int j = 0; j < numColIndices; j++) newColIndices[j] = oldColIndices[j];
1993  }
1994 
1995 //#ifdef EPETRA_HAVE_OMP
1996 //#pragma omp parallel for default(none) shared(all_indices_values,indices_values)
1997 //#endif
1998 // for(int ii=0; ii<numMyNonzeros; ++ii) {
1999 // all_indices_values[ii] = indices_values[ii];
2000 // }
2001  }
2002  }
2003 
2004  // Delete unneeded storage
2006  delete [] Data.Indices_; Data.Indices_=0;
2007  Data.SortedEntries_.clear();
2008 
2009  SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
2011 
2012 /*
2013 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2014  All_IndicesPlus1(); // see if preemptively calling this improves Multiply timing.
2015 #endif
2016 */
2017 
2018  return(0);
2019 }
2020 
2021 //==============================================================================
2022 template<typename int_type>
2023 int Epetra_CrsGraph::ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
2024 {
2025  int j;
2026 
2027  int locRow = LRID(Row); // Normalize row range
2028 
2029  if(locRow < 0 || locRow >= NumMyBlockRows())
2030  EPETRA_CHK_ERR(-1); // Not in Row range
2031 
2032  NumIndices = NumMyIndices(locRow);
2033  if(LenOfIndices < NumIndices)
2034  EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
2035 
2036  if(IndicesAreLocal())
2037  {
2038  int * srcIndices = TIndices<int>(locRow);
2039  // static_cast is ok because global indices were created from int values and hence must fit ints
2040  for(j = 0; j < NumIndices; j++)
2041  targIndices[j] = static_cast<int_type>(GCID64(srcIndices[j]));
2042  }
2043  else
2044  {
2045  int_type * srcIndices = TIndices<int_type>(locRow);
2046  for(j = 0; j < NumIndices; j++)
2047  targIndices[j] = srcIndices[j];
2048  }
2049 
2050  return(0);
2051 }
2052 
2053 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2054 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
2055 {
2056  if(RowMap().GlobalIndicesInt())
2057  return ExtractGlobalRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
2058  else
2059  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy int version called for a graph that is not int.", -1);
2060 }
2061 #endif
2062 
2063 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2064 int Epetra_CrsGraph::ExtractGlobalRowCopy(long long Row, int LenOfIndices, int& NumIndices, long long* targIndices) const
2065 {
2066  if(RowMap().GlobalIndicesLongLong())
2067  return ExtractGlobalRowCopy<long long>(Row, LenOfIndices, NumIndices, targIndices);
2068  else
2069  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy long long version called for a graph that is not long long.", -1);
2070 }
2071 #endif
2072 
2073 //==============================================================================
2074 template<typename int_type>
2075 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
2076 {
2077  int j;
2078 
2079  if(Row < 0 || Row >= NumMyBlockRows())
2080  EPETRA_CHK_ERR(-1); // Not in Row range
2081 
2082  NumIndices = NumMyIndices(Row);
2083  if(LenOfIndices < NumIndices)
2084  EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
2085 
2086  if(IndicesAreGlobal())
2087  EPETRA_CHK_ERR(-3); // There are no local indices yet
2088 
2089  int * srcIndices = TIndices<int>(Row);
2090  for(j = 0; j < NumIndices; j++)
2091  targIndices[j] = srcIndices[j];
2092 
2093  return(0);
2094 }
2095 
2096 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
2097 {
2098  if(RowMap().GlobalIndicesTypeValid())
2099  return ExtractMyRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
2100  else
2101  throw ReportError("Epetra_CrsGraph::ExtractMyRowCopy graph global index type unknown.", -1);
2102 }
2103 
2104 //==============================================================================
2105 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2106 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const
2107 {
2108  if(!RowMap().GlobalIndicesInt())
2109  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView int version called for a graph that is not int.", -1);
2110 
2111  int locRow = LRID(Row); // Normalize row range
2112 
2113  if(locRow < 0 || locRow >= NumMyBlockRows())
2114  EPETRA_CHK_ERR(-1); // Not in Row range
2115 
2116  if(IndicesAreLocal())
2117  EPETRA_CHK_ERR(-2); // There are no global indices
2118 
2119  NumIndices = NumMyIndices(locRow);
2120 
2121  targIndices = TIndices<int>(locRow);
2122 
2123  return(0);
2124 }
2125 #endif
2126 //==============================================================================
2127 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2128 int Epetra_CrsGraph::ExtractGlobalRowView(long long Row, int& NumIndices, long long*& targIndices) const
2129 {
2130  if(!RowMap().GlobalIndicesLongLong())
2131  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView long long version called for a graph that is not long long.", -1);
2132 
2133  int locRow = LRID(Row); // Normalize row range
2134 
2135  if(locRow < 0 || locRow >= NumMyBlockRows())
2136  EPETRA_CHK_ERR(-1); // Not in Row range
2137 
2138  if(IndicesAreLocal())
2139  EPETRA_CHK_ERR(-2); // There are no global indices
2140 
2141  NumIndices = NumMyIndices(locRow);
2142 
2143  targIndices = TIndices<long long>(locRow);
2144 
2145  return(0);
2146 }
2147 #endif
2148 //==============================================================================
2149 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const
2150 {
2151  if(Row < 0 || Row >= NumMyBlockRows())
2152  EPETRA_CHK_ERR(-1); // Not in Row range
2153 
2154  if(IndicesAreGlobal())
2155  EPETRA_CHK_ERR(-2); // There are no local indices
2156 
2157  NumIndices = NumMyIndices(Row);
2158 
2159  targIndices = TIndices<int>(Row);
2160 
2161  return(0);
2162 }
2163 
2164 //==============================================================================
2165 int Epetra_CrsGraph::NumGlobalIndices(long long Row) const {
2166 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
2167  int locRow = LRID((int) Row);
2168 #else
2169  int locRow = LRID(Row);
2170 #endif
2171  if(locRow != -1)
2172  return(NumMyIndices(locRow));
2173  else
2174  return(0); // No indices for this row on this processor
2175 }
2176 
2177 //==============================================================================
2179 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
2180  int locRow = LRID((int) Row);
2181 #else
2182  int locRow = LRID(Row);
2183 #endif
2184  if(locRow != -1)
2185  return(NumAllocatedMyIndices(locRow));
2186  else
2187  return(0); // No indices allocated for this row on this processor
2188 }
2189 
2190 //==============================================================================
2192 {
2193  if (RowMap().PointSameAs(newmap)) {
2194  Epetra_DistObject::Map_ = newmap;
2195  CrsGraphData_->RowMap_ = newmap;
2197  return(0);
2198  }
2199 
2200  return(-1);
2201 }
2202 
2203 //==============================================================================
2205 {
2206  if (!HaveColMap() && !IndicesAreLocal() && !IndicesAreGlobal() && newmap.GlobalIndicesTypeMatch(RowMap())) {
2207  CrsGraphData_->ColMap_ = newmap;
2213  CrsGraphData_->NumMyCols_ = newmap.NumMyPoints();
2214  CrsGraphData_->HaveColMap_ = true;
2215  return(0);
2216  }
2217 
2218  if(ColMap().PointSameAs(newmap)) {
2219  CrsGraphData_->ColMap_ = newmap;
2221  return(0);
2222  }
2223 
2224  return(-1);
2225 }
2226 
2227 
2228 //==============================================================================
2229 int Epetra_CrsGraph::ReplaceDomainMapAndImporter(const Epetra_BlockMap& NewDomainMap, const Epetra_Import * NewImporter) {
2230  int rv=0;
2231  if( !NewImporter && ColMap().SameAs(NewDomainMap)) {
2232  CrsGraphData_->DomainMap_ = NewDomainMap;
2233  delete CrsGraphData_->Importer_;
2234  CrsGraphData_->Importer_ = 0;
2235  }
2236  else if(NewImporter && ColMap().SameAs(NewImporter->TargetMap()) && NewDomainMap.SameAs(NewImporter->SourceMap())) {
2237  CrsGraphData_->DomainMap_ = NewDomainMap;
2238  delete CrsGraphData_->Importer_;
2239  CrsGraphData_->Importer_ = new Epetra_Import(*NewImporter);
2240  }
2241  else
2242  rv=-1;
2243  return rv;
2244 }
2245 
2246 //==============================================================================
2248  const Epetra_BlockMap *newDomainMap=0, *newRangeMap=0, *newColMap=0;
2249  Epetra_Import * newImport=0;
2250  Epetra_Export * newExport=0;
2251 
2252  const Epetra_Comm *newComm = newMap ? &newMap->Comm() : 0;
2253 
2254  if(DomainMap().SameAs(RowMap())) {
2255  // Common case: original domain and row Maps are identical.
2256  // In that case, we need only replace the original domain Map
2257  // with the new Map. This ensures that the new domain and row
2258  // Maps _stay_ identical.
2259  newDomainMap = newMap;
2260  }
2261  else
2262  newDomainMap = DomainMap().ReplaceCommWithSubset(newComm);
2263 
2264  if(RangeMap().SameAs(RowMap())){
2265  // Common case: original range and row Maps are identical. In
2266  // that case, we need only replace the original range Map with
2267  // the new Map. This ensures that the new range and row Maps
2268  // _stay_ identical.
2269  newRangeMap = newMap;
2270  }
2271  else
2272  newRangeMap = RangeMap().ReplaceCommWithSubset(newComm);
2273 
2274  newColMap=ColMap().ReplaceCommWithSubset(newComm);
2275 
2276  if(newComm) {
2277  // (Re)create the Export and / or Import if necessary.
2278  //
2279  // The operations below are collective on the new communicator.
2280  //
2281  if(RangeMap().DataPtr() != RowMap().DataPtr())
2282  newExport = new Epetra_Export(*newMap,*newRangeMap);
2283 
2284  if(DomainMap().DataPtr() != ColMap().DataPtr())
2285  newImport = new Epetra_Import(*newColMap,*newDomainMap);
2286  }
2287 
2288  // CrsGraphData things
2289  if(CrsGraphData_->ReferenceCount() !=1)
2290  throw ReportError("Epetra_CrsGraph::RemoveEmptyProcessesInPlace does not work for shared CrsGraphData_",-2);
2291 
2292  // Dummy map for the non-active procs
2293 #if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
2294  Epetra_Map dummy;
2295 #else
2296  Epetra_SerialComm SComm;
2297  Epetra_Map dummy(0,0,SComm);
2298 #endif
2299 
2300  delete CrsGraphData_->Importer_;
2301  delete CrsGraphData_->Exporter_;
2302 
2303 
2304  CrsGraphData_->RowMap_ = newMap ? *newMap : dummy;
2305  CrsGraphData_->ColMap_ = newColMap ? *newColMap : dummy;
2306  CrsGraphData_->DomainMap_ = newDomainMap ? *newDomainMap : dummy;
2307  CrsGraphData_->RangeMap_ = newRangeMap ? *newRangeMap : dummy;
2308  CrsGraphData_->Importer_ = newImport;
2309  CrsGraphData_->Exporter_ = newExport;
2310 
2311  // Epetra_DistObject things
2312  if(newMap) {
2313  Map_ = *newMap;
2314  Comm_ = &newMap->Comm();
2315  }
2316 
2317  // Cleanup (newRowMap is always newMap, so don't delete that)
2318  if(newColMap != newMap) delete newColMap;
2319  if(newDomainMap != newMap) delete newDomainMap;
2320  if(newRangeMap != newMap) delete newRangeMap;
2321 
2322  return(0);
2323 }
2324 
2325 // private =====================================================================
2327  try {
2328  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
2329  if(!A.GlobalConstantsComputed())
2330  EPETRA_CHK_ERR(-1); // Must have global constants to proceed
2331  }
2332  catch(...) {
2333  return(0); // No error at this point, object could be a RowMatrix
2334  }
2335  return(0);
2336 }
2337 
2338 // private =====================================================================
2340  int NumSameIDs,
2341  int NumPermuteIDs,
2342  int* PermuteToLIDs,
2343  int* PermuteFromLIDs,
2344  const Epetra_OffsetIndex * Indexor,
2345  Epetra_CombineMode CombineMode)
2346 {
2347  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2348  throw ReportError("Epetra_CrsGraph::CopyAndPermute: Incoming global index type does not match the one for *this",-1);
2349 
2350  try {
2351  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2352  EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
2353  PermuteFromLIDs,Indexor,CombineMode));
2354  }
2355  catch(...) {
2356  try {
2357  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2358  EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
2359  PermuteFromLIDs,Indexor,CombineMode));
2360  }
2361  catch(...) {
2362  EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
2363  }
2364  }
2365 
2366  return(0);
2367 }
2368 
2369 // private =====================================================================
2370 template<typename int_type>
2372  int NumSameIDs,
2373  int NumPermuteIDs,
2374  int* PermuteToLIDs,
2375  int* PermuteFromLIDs,
2376  const Epetra_OffsetIndex * Indexor,
2377  Epetra_CombineMode CombineMode)
2378 {
2379  (void)Indexor;
2380  (void)CombineMode;
2381  int i;
2382  int j;
2383  int NumIndices;
2384  int FromRow;
2385  int_type ToRow;
2386  int maxNumIndices = A.MaxNumEntries();
2387  Epetra_IntSerialDenseVector local_indices_vec;
2388 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2389  Epetra_LongLongSerialDenseVector global_indices_vec;
2390 #endif
2391  Epetra_SerialDenseVector Values;
2392 
2393  int* local_indices = 0;
2394  int_type* global_indices = 0;
2395 
2396  if(maxNumIndices > 0) {
2397  local_indices_vec.Size(maxNumIndices);
2398  local_indices = local_indices_vec.Values();
2399 
2400  if(A.RowMatrixRowMap().GlobalIndicesLongLong())
2401  {
2402 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2403  global_indices_vec.Size(maxNumIndices);
2404  global_indices = reinterpret_cast<int_type*>(global_indices_vec.Values());
2405 #else
2406  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: GlobalIndicesLongLong but no API for long long",-1);
2407 #endif
2408  }
2409  else
2410  {
2411  global_indices = reinterpret_cast<int_type*>(local_indices);
2412  }
2413 
2414  Values.Size(maxNumIndices); // Must extract values even though we discard them
2415  }
2416 
2417  const Epetra_Map& rowMap = A.RowMatrixRowMap();
2418  const Epetra_Map& colMap = A.RowMatrixColMap();
2419 
2420  // Do copy first
2421  for(i = 0; i < NumSameIDs; i++) {
2422  ToRow = (int) rowMap.GID64(i);
2423  EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumIndices, NumIndices, Values.Values(), local_indices));
2424  for(j = 0; j < NumIndices; j++)
2425  global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
2426  // Place into target graph.
2427  int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices);
2428  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2429  }
2430 
2431  // Do local permutation next
2432  for(i = 0; i < NumPermuteIDs; i++) {
2433  FromRow = PermuteFromLIDs[i];
2434  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2435  EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumIndices, NumIndices, Values.Values(), local_indices));
2436  for(j = 0; j < NumIndices; j++)
2437  global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
2438  int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices); // Place into target graph.
2439  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2440  }
2441 
2442  return(0);
2443 }
2444 
2446  int NumSameIDs,
2447  int NumPermuteIDs,
2448  int* PermuteToLIDs,
2449  int* PermuteFromLIDs,
2450  const Epetra_OffsetIndex * Indexor,
2451  Epetra_CombineMode CombineMode)
2452 {
2453  if(!A.RowMatrixRowMap().GlobalIndicesTypeMatch(RowMap()))
2454  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
2455 
2456 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2457  if(A.RowMatrixRowMap().GlobalIndicesInt())
2458  return CopyAndPermuteRowMatrix<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor,CombineMode);
2459  else
2460 #endif
2461 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2462  if(A.RowMatrixRowMap().GlobalIndicesLongLong())
2463  return CopyAndPermuteRowMatrix<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor,CombineMode);
2464  else
2465 #endif
2466  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Unable to determine global index type of map", -1);
2467 }
2468 
2469 // private =====================================================================
2470 template<typename int_type>
2472  int NumSameIDs,
2473  int NumPermuteIDs,
2474  int* PermuteToLIDs,
2475  int* PermuteFromLIDs,
2476  const Epetra_OffsetIndex * Indexor,
2477  Epetra_CombineMode CombineMode)
2478 {
2479  (void)Indexor;
2480  (void)CombineMode;
2481  int i;
2482  int_type Row;
2483  int NumIndices;
2484  int_type* indices = 0;
2485  int_type FromRow, ToRow;
2486  int maxNumIndices = A.MaxNumIndices();
2487  Epetra_IntSerialDenseVector int_IndicesVector;
2488 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2489  Epetra_LongLongSerialDenseVector LL_IndicesVector;
2490 #endif
2491 
2492  if(maxNumIndices > 0 && A.IndicesAreLocal()) {
2493  if(A.RowMap().GlobalIndicesInt())
2494  {
2495  int_IndicesVector.Size(maxNumIndices);
2496  indices = reinterpret_cast<int_type*>(int_IndicesVector.Values());
2497  }
2498  else if(A.RowMap().GlobalIndicesLongLong())
2499  {
2500 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2501  LL_IndicesVector.Size(maxNumIndices);
2502  indices = reinterpret_cast<int_type*>(LL_IndicesVector.Values());
2503 #else
2504  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2505 #endif
2506  }
2507  }
2508 
2509  // Do copy first
2510  if(NumSameIDs > 0) {
2511  if(A.IndicesAreLocal()) {
2512  for(i = 0; i < NumSameIDs; i++) {
2513  Row = (int_type) GRID64(i);
2514  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumIndices, NumIndices, indices));
2515  // Place into target graph.
2516  int ierr = InsertGlobalIndices(Row, NumIndices, indices);
2517  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2518  }
2519  }
2520  else { // A.IndiceAreGlobal()
2521  for(i = 0; i < NumSameIDs; i++) {
2522  Row = (int_type) GRID64(i);
2523  EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, indices));
2524  // Place into target graph.
2525  int ierr = InsertGlobalIndices(Row, NumIndices, indices);
2526  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2527  }
2528  }
2529  }
2530 
2531  // Do local permutation next
2532  if(NumPermuteIDs > 0) {
2533  if(A.IndicesAreLocal()) {
2534  for(i = 0; i < NumPermuteIDs; i++) {
2535  FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
2536  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2537  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2538  // Place into target graph.
2539  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2540  if (ierr < 0) EPETRA_CHK_ERR(ierr);
2541  }
2542  }
2543  else { // A.IndiceAreGlobal()
2544  for(i = 0; i < NumPermuteIDs; i++) {
2545  FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
2546  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2547  EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, indices));
2548  // Place into target graph.
2549  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2550  if (ierr < 0) EPETRA_CHK_ERR(ierr);
2551  }
2552  }
2553  }
2554 
2555  return(0);
2556 }
2557 
2559  int NumSameIDs,
2560  int NumPermuteIDs,
2561  int* PermuteToLIDs,
2562  int* PermuteFromLIDs,
2563  const Epetra_OffsetIndex * Indexor,
2564  Epetra_CombineMode CombineMode)
2565 {
2566  if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
2567  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Incoming global index type does not match the one for *this",-1);
2568 
2569  if(A.RowMap().GlobalIndicesInt())
2570 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2571  return CopyAndPermuteCrsGraph<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2572 #else
2573  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesInt but no API for it.",-1);
2574 #endif
2575 
2576  if(A.RowMap().GlobalIndicesLongLong())
2577 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2578  return CopyAndPermuteCrsGraph<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2579 #else
2580  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2581 #endif
2582 
2583  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Unable to determine global index type of map", -1);
2584 }
2585 
2586 // private =====================================================================
2588  int NumExportIDs,
2589  int* ExportLIDs,
2590  int& LenExports,
2591  char*& Exports,
2592  int& SizeOfPacket,
2593  int * Sizes,
2594  bool& VarSizes,
2595  Epetra_Distributor& Distor)
2596 {
2597  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2598  throw ReportError("Epetra_CrsGraph::PackAndPrepare: Incoming global index type does not match the one for *this",-1);
2599 
2600  int globalMaxNumIndices = 0;
2601  int TotalSendSize = 0;
2602 
2603  VarSizes = true;
2604 
2605  if(Source.Map().GlobalIndicesInt())
2606  SizeOfPacket = (int)sizeof(int);
2607  else if(Source.Map().GlobalIndicesLongLong())
2608  SizeOfPacket = (int)sizeof(long long);
2609  else
2610  throw ReportError("Epetra_CrsGraph::PackAndPrepare: Unable to determine source global index type",-1);
2611 
2612  if(NumExportIDs <= 0) return(0);
2613 
2614  try {
2615  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2616  globalMaxNumIndices = A.GlobalMaxNumIndices();
2617  for( int i = 0; i < NumExportIDs; ++i )
2618  {
2619  Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
2620  TotalSendSize += Sizes[i];
2621  }
2622  }
2623  catch(...) {
2624  try {
2625  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2626  int maxNumIndices = A.MaxNumEntries();
2627  A.Comm().MaxAll(&maxNumIndices, &globalMaxNumIndices, 1);
2628  for( int i = 0; i < NumExportIDs; ++i )
2629  {
2630  int NumEntries;
2631  A.NumMyRowEntries( ExportLIDs[i], NumEntries );
2632  Sizes[i] = (NumEntries + 2);
2633  TotalSendSize += Sizes[i];
2634  }
2635  }
2636  catch(...) {
2637  EPETRA_CHK_ERR(-1); // Bad cast
2638  }
2639  }
2640 
2641  CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
2642 
2643  try {
2644  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2645  EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
2646  SizeOfPacket, Sizes, VarSizes, Distor));
2647  }
2648  catch(...) {
2649  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2650  EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
2651  SizeOfPacket, Sizes, VarSizes, Distor));
2652  }
2653  return(0);
2654 }
2655 
2656 // private =====================================================================
2658  int NumExportIDs,
2659  int* ExportLIDs,
2660  int& LenExports,
2661  char*& Exports,
2662  int& SizeOfPacket,
2663  int* Sizes,
2664  bool& VarSizes,
2665  Epetra_Distributor& Distor)
2666 {
2667  if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
2668  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Incoming global index type does not match the one for *this",-1);
2669 
2670  (void)LenExports;
2671  (void)SizeOfPacket;
2672  (void)Sizes;
2673  (void)VarSizes;
2674  (void)Distor;
2675  int i;
2676  int NumIndices;
2677 
2678  // Each segment of Exports will be filled by a packed row of information for each row as follows:
2679  // 1st int: GRID of row where GRID is the global row ID for the source graph
2680  // next int: NumIndices, Number of indices in row.
2681  // next NumIndices: The actual indices for the row.
2682  // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
2683  // sized segments for current communication routines.
2684  int maxNumIndices = A.MaxNumIndices();
2685  //if( maxNumIndices ) indices = new int[maxNumIndices];
2686 
2687  if(A.RowMap().GlobalIndicesInt()) {
2688  int* indices = 0;
2689  int* intptr = (int*) Exports;
2690  int FromRow;
2691  for(i = 0; i < NumExportIDs; i++) {
2692  FromRow = (int) A.GRID64(ExportLIDs[i]);
2693  *intptr = FromRow;
2694  indices = intptr + 2;
2695  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2696  intptr[1] = NumIndices; // Load second slot of segment
2697  intptr += (NumIndices+2); // Point to next segment
2698  }
2699  }
2700  else if(A.RowMap().GlobalIndicesLongLong()) {
2701 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2702  long long* indices = 0;
2703  long long* LLptr = (long long*) Exports;
2704  long long FromRow;
2705  for(i = 0; i < NumExportIDs; i++) {
2706  FromRow = A.GRID64(ExportLIDs[i]);
2707  *LLptr = FromRow;
2708  indices = LLptr + 2;
2709  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2710  LLptr[1] = NumIndices; // Load second slot of segment
2711  LLptr += (NumIndices+2); // Point to next segment
2712  }
2713 #else
2714  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2715 #endif
2716  }
2717  else
2718  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Unable to determine source global index type",-1);
2719 
2720  //if( indices ) delete [] indices;
2721 
2722  return(0);
2723 }
2724 
2725 // private =====================================================================
2727  int NumExportIDs,
2728  int* ExportLIDs,
2729  int& LenExports,
2730  char*& Exports,
2731  int& SizeOfPacket,
2732  int* Sizes,
2733  bool& VarSizes,
2734  Epetra_Distributor& Distor)
2735 {
2736  if(!A.Map().GlobalIndicesTypeMatch(RowMap()))
2737  throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Incoming global index type does not match the one for *this",-1);
2738 
2739  (void)LenExports;
2740  (void)SizeOfPacket;
2741  (void)Sizes;
2742  (void)VarSizes;
2743  (void)Distor;
2744  int i;
2745  int j;
2746  int NumIndices;
2747  Epetra_SerialDenseVector Values;
2748 
2749  // Each segment of Exports will be filled by a packed row of information for each row as follows:
2750  // 1st int: GRID of row where GRID is the global row ID for the source graph
2751  // next int: NumIndices, Number of indices in row.
2752  // next NumIndices: The actual indices for the row.
2753  // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
2754  // sized segments for current communication routines.
2755  int maxNumIndices = A.MaxNumEntries();
2756  if(maxNumIndices > 0) {
2757  Values.Size(maxNumIndices);
2758 // indices = new int[maxNumIndices];
2759  }
2760  const Epetra_Map& rowMap = A.RowMatrixRowMap();
2761  const Epetra_Map& colMap = A.RowMatrixColMap();
2762 
2763  if(rowMap.GlobalIndicesInt() && colMap.GlobalIndicesInt()) {
2764  int* indices = 0;
2765  int FromRow;
2766  int* intptr = (int*) Exports;
2767  for(i = 0; i < NumExportIDs; i++) {
2768  FromRow = (int) rowMap.GID64(ExportLIDs[i]);
2769  *intptr = FromRow;
2770  indices = intptr + 2;
2771  EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), indices));
2772  for(j = 0; j < NumIndices; j++) indices[j] = (int) colMap.GID64(indices[j]); // convert to GIDs
2773  intptr[1] = NumIndices; // Load second slot of segment
2774  intptr += (NumIndices+2); // Point to next segment
2775  }
2776  }
2777  else if(rowMap.GlobalIndicesLongLong() && colMap.GlobalIndicesLongLong()) {
2778  // Bytes of Exports:
2779  // 12345678.12345678....12345678.12345678 ("." means no spaces)
2780  // FromRow NumIndices id1 id2 id3 id4 <-- before converting to GIDs
2781  // FromRow NumIndices | gid1 | | gid2 | <-- after converting to GIDs
2782 
2783  long long* LL_indices = 0;
2784  long long FromRow;
2785  long long* LLptr = (long long*) Exports;
2786  for(i = 0; i < NumExportIDs; i++) {
2787  FromRow = rowMap.GID64(ExportLIDs[i]);
2788  *LLptr = FromRow;
2789  LL_indices = LLptr + 2;
2790  int* int_indices = reinterpret_cast<int*>(LL_indices);
2791  EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), int_indices));
2792 
2793  // convert to GIDs, start from right.
2794  for(j = NumIndices; j > 0;) {
2795  --j;
2796  LL_indices[j] = colMap.GID64(int_indices[j]);
2797  }
2798 
2799  LLptr[1] = NumIndices; // Load second slot of segment
2800  LLptr += (NumIndices+2); // Point to next segment
2801  }
2802  }
2803  else
2804  throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Unable to determine source global index type",-1);
2805 
2806 
2807 // if( indices ) delete [] indices;
2808 
2809  return(0);
2810 }
2811 
2812 // private =====================================================================
2814  int NumImportIDs,
2815  int* ImportLIDs,
2816  int LenImports,
2817  char* Imports,
2818  int& SizeOfPacket,
2819  Epetra_Distributor& Distor,
2820  Epetra_CombineMode CombineMode,
2821  const Epetra_OffsetIndex * Indexor)
2822 {
2823  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2824  throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Incoming global index type does not match the one for *this",-1);
2825 
2826  (void)Source;
2827  (void)LenImports;
2828  (void)SizeOfPacket;
2829  (void)Distor;
2830  (void)CombineMode;
2831  (void)Indexor;
2832  if(NumImportIDs <= 0)
2833  return(0);
2834 
2835  // Unpack it...
2836 
2837  // Each segment of Sends will be filled by a packed row of information for each row as follows:
2838  // 1st int: GRID of row where GRID is the global row ID for the source graph
2839  // next int: NumIndices, Number of indices in row.
2840  // next NumIndices: The actual indices for the row.
2841 
2842 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2843  if(Source.Map().GlobalIndicesInt()) {
2844  int NumIndices;
2845  int i;
2846  int* indices;
2847  int ToRow;
2848  int* intptr = (int*) Imports;
2849  for(i = 0; i < NumImportIDs; i++) {
2850  ToRow = (int) GRID64(ImportLIDs[i]);
2851  assert((intptr[0])==ToRow); // Sanity check
2852  NumIndices = intptr[1];
2853  indices = intptr + 2;
2854  // Insert indices
2855  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2856  if(ierr < 0)
2857  EPETRA_CHK_ERR(ierr);
2858  intptr += (NumIndices+2); // Point to next segment
2859  }
2860  }
2861  else
2862 #endif
2863 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2864  if(Source.Map().GlobalIndicesLongLong()) {
2865  int NumIndices;
2866  int i;
2867  long long* indices;
2868  long long ToRow;
2869  long long* LLptr = (long long*) Imports;
2870  for(i = 0; i < NumImportIDs; i++) {
2871  ToRow = GRID64(ImportLIDs[i]);
2872  assert((LLptr[0])==ToRow); // Sanity check
2873  NumIndices = (int) LLptr[1];
2874  indices = LLptr + 2;
2875  // Insert indices
2876  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2877  if(ierr < 0)
2878  EPETRA_CHK_ERR(ierr);
2879  LLptr += (NumIndices+2); // Point to next segment
2880  }
2881  }
2882  else
2883 #endif
2884  throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Unable to determine source global index type",-1);
2885 
2886  //destroy buffers since this operation is usually only done once
2887  if( LenExports_ ) {
2888  delete [] Exports_;
2889  Exports_ = 0;
2890  LenExports_ = 0;
2891  }
2892  if( LenImports_ ) {
2893  delete [] Imports_;
2894  Imports_ = 0;
2895  LenImports_ = 0;
2896  }
2897 
2898  return(0);
2899 }
2900 
2901 // protected ===================================================================
2903  int mineComputed = 0;
2904  int allComputed;
2906  mineComputed = 1;
2907  RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
2908  // If allComputed comes back zero then at least one PE need global constants recomputed.
2909  return(allComputed==1);
2910 }
2911 
2912 // private =====================================================================
2914  int myIndicesAreLocal = 0;
2915  int myIndicesAreGlobal = 0;
2917  myIndicesAreLocal = 1;
2919  myIndicesAreGlobal = 1;
2920  int allIndicesAreLocal;
2921  int allIndicesAreGlobal;
2922  RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
2923  RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
2924  CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
2925  CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1); // If indices are global on one PE should be local on all
2926 }
2927 
2928 //==============================================================================
2929 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2930 int *Epetra_CrsGraph::All_IndicesPlus1() const {
2931  // This functionality is needed because MKL "sparse matrix" "dense matrix"
2932  // functions do not work with column-based dense storage and zero-based
2933  // sparse storage. So add "1" to indices and save duplicate data. This means
2934  // we will use one-based indices. This does not affect sparse-matrix and vector
2935  // operations.
2936 
2937  int* ptr = 0;
2938  if (!StorageOptimized()) {
2939  throw ReportError("Epetra_CrsGraph: int *All_IndicesPlus1() cannot be called when StorageOptimized()==false", -1);
2940  }
2941  else {
2943 
2944  if(!vec.Length()) {
2945  int* indices = All_Indices();
2947  ptr = vec.Values();
2948  for(int i = 0; i < CrsGraphData_->NumMyNonzeros_; ++i)
2949  ptr[i] = indices[i] + 1;
2950  }
2951  else {
2952  ptr = vec.Values();
2953  }
2954  }
2955  return ptr;
2956 }
2957 #endif // defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2958 
2959 //==============================================================================
2960 void Epetra_CrsGraph::Print (std::ostream& os) const {
2961  int MyPID = RowMap().Comm().MyPID();
2962  int NumProc = RowMap().Comm().NumProc();
2963 
2964  for(int iproc = 0; iproc < NumProc; iproc++) {
2965  if(MyPID == iproc) {
2966  if(MyPID == 0) {
2967  os << "\nNumber of Global Block Rows = " << NumGlobalBlockRows64() << std::endl;
2968  os << "Number of Global Block Cols = " << NumGlobalBlockCols64() << std::endl;
2969  os << "Number of Global Block Diags = " << NumGlobalBlockDiagonals64() << std::endl;
2970  os << "Number of Global Entries = " << NumGlobalEntries64() << std::endl;
2971  os << "\nNumber of Global Rows = " << NumGlobalRows64() << std::endl;
2972  os << "Number of Global Cols = " << NumGlobalCols64() << std::endl;
2973  os << "Number of Global Diagonals = " << NumGlobalDiagonals64() << std::endl;
2974  os << "Number of Global Nonzeros = " << NumGlobalNonzeros64() << std::endl;
2975  os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim() << std::endl;
2976  os << "Global Maximum Block Col Dim = " << GlobalMaxColDim() << std::endl;
2977  os << "Global Maximum Num Indices = " << GlobalMaxNumIndices() << std::endl;
2978  if(LowerTriangular()) os << " ** Matrix is Lower Triangular **" << std::endl;
2979  if(UpperTriangular()) os << " ** Matrix is Upper Triangular **" << std::endl;
2980  if(NoDiagonal()) os << " ** Matrix has no diagonal **" << std::endl << std::endl;
2981  }
2982  os << "\nNumber of My Block Rows = " << NumMyBlockRows() << std::endl;
2983  os << "Number of My Block Cols = " << NumMyBlockCols() << std::endl;
2984  os << "Number of My Block Diags = " << NumMyBlockDiagonals() << std::endl;
2985  os << "Number of My Entries = " << NumMyEntries() << std::endl;
2986  os << "\nNumber of My Rows = " << NumMyRows() << std::endl;
2987  os << "Number of My Cols = " << NumMyCols() << std::endl;
2988  os << "Number of My Diagonals = " << NumMyDiagonals() << std::endl;
2989  os << "Number of My Nonzeros = " << NumMyNonzeros() << std::endl;
2990  os << "\nMy Maximum Block Row Dim = " << MaxRowDim() << std::endl;
2991  os << "My Maximum Block Col Dim = " << MaxColDim() << std::endl;
2992  os << "My Maximum Num Indices = " << MaxNumIndices() << std::endl << std::endl;
2993 
2994  int NumMyBlockRows1 = NumMyBlockRows();
2995  int MaxNumIndices1 = MaxNumIndices();
2996  Epetra_IntSerialDenseVector Indices1_int(MaxNumIndices1);
2997 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2998  Epetra_LongLongSerialDenseVector Indices1_LL(MaxNumIndices1);
2999 #endif
3000 
3001  if(RowMap().GlobalIndicesInt()) {
3002  Indices1_int.Resize(MaxNumIndices1);
3003  }
3004  else if(RowMap().GlobalIndicesLongLong()) {
3005 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
3006  Indices1_LL.Resize(MaxNumIndices1);
3007 #else
3008  throw ReportError("Epetra_CrsGraph::Print: GlobalIndicesLongLong but no long long API",-1);
3009 #endif
3010  }
3011  else
3012  throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
3013 
3014  int NumIndices1;
3015  int i;
3016  int j;
3017 
3018  os.width(14);
3019  os << " Row Index "; os << " ";
3020  for(j = 0; j < MaxNumIndices(); j++) {
3021  os.width(12);
3022  os << "Col Index"; os << " ";
3023  }
3024  os << std::endl;
3025  for(i = 0; i < NumMyBlockRows1; i++) {
3026  if(RowMap().GlobalIndicesInt()) {
3027  int Row = (int) GRID64(i); // Get global row number
3028  ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_int.Values());
3029  os.width(14);
3030  os << Row ; os << " ";
3031  for(j = 0; j < NumIndices1 ; j++) {
3032  os.width(12);
3033  os << Indices1_int[j]; os << " ";
3034  }
3035  os << std::endl;
3036  }
3037  else if(RowMap().GlobalIndicesLongLong()) {
3038 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
3039  long long Row = GRID64(i); // Get global row number
3040  ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_LL.Values());
3041  os.width(14);
3042  os << Row ; os << " ";
3043  for(j = 0; j < NumIndices1 ; j++) {
3044  os.width(12);
3045  os << Indices1_LL[j]; os << " ";
3046  }
3047  os << std::endl;
3048 #else
3049  throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
3050 #endif
3051  }
3052  }
3053  os << std::flush;
3054  }
3055  // Do a few global ops to give I/O a chance to complete
3056  RowMap().Comm().Barrier();
3057  RowMap().Comm().Barrier();
3058  RowMap().Comm().Barrier();
3059  }
3060 }
3061 
3062 //==============================================================================
3064  if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
3065  return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
3066 
3067  CleanupData();
3068  CrsGraphData_ = Source.CrsGraphData_;
3070 
3071  return(*this);
3072 }
void SetFilled(bool Flag)
int MakeViewOf(const Epetra_IntSerialDenseVector &Source)
Reset an existing IntSerialDenseVector to point to another Vector.
int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
Perform any packing or preparation required for call to DoTransfer().
const Epetra_Export * Exporter_
const Epetra_Import * Importer() const
Returns the Importer associated with this graph.
long long NumGlobalEntries64() const
bool FindMyIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
bool HaveColMap() const
Returns true if we have a well-defined ColMap, and returns false otherwise.
int MakeImportExport()
called by FillComplete (and TransformToLocal)
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
bool MyLID(int lid) const
Returns true if the LID passed in belongs to the calling processor in this map, otherwise returns fal...
virtual ~Epetra_CrsGraph()
Epetra_CrsGraph Destructor.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
Epetra_BlockMap RangeMap_
long long NumGlobalBlockCols64() const
Epetra_IntSerialDenseVector NumIndicesPerRow_
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int SortIndices()
Sort column indices, row-by-row, in ascending order.
int Allocate(const int *NumIndicesPerRow, int Inc, bool StaticProfile)
int Resize(int Length_in)
Resize a Epetra_IntSerialDenseVector object.
int ReferenceCount() const
Get reference count.
Definition: Epetra_Data.cpp:71
bool FindGlobalIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
Epetra_CrsGraph & operator=(const Epetra_CrsGraph &Source)
Assignment operator.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
bool ConstantElementSize() const
Returns true if map has constant element size.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
long long * Values()
Returns pointer to the values in vector.
int TRemoveGlobalIndices(long long Row)
void DecrementReferenceCount()
Decrement reference count.
Definition: Epetra_Data.cpp:66
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int ReAllocateAndCast(char *&UserPtr, int &Length, const int IntPacketSizeTimesNumTrans)
called by PackAndPrepare
int PackAndPrepareCrsGraph(const Epetra_CrsGraph &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
int NumMyDiagonals() const
Returns the number of diagonal entries in the local graph, based on global row/column index compariso...
int NumAllocatedGlobalIndices(long long Row) const
Returns the allocated number of nonzero entries in specified global row on this processor.
value_type Get(const long long key)
int NumMyNonzeros() const
Returns the number of indices in the local graph.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this graph.
int Size(int Length_in)
Set length of a Epetra_IntSerialDenseVector object; init values to zero.
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this importer.
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified local row of the graph.
int NumAllocatedMyIndices(int Row) const
Returns the allocated number of nonzero entries in specified local row on this processor.
const Epetra_CrsGraphData * DataPtr() const
Returns a pointer to the CrsGraphData instance this CrsGraph uses.
double * Values() const
Returns pointer to the values in vector.
long long IndexBase64() const
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:413
Epetra_BlockMap DomainMap_
#define EPETRA_CHK_ERR(a)
int GlobalMaxColDim() const
Returns the max column dimension of block entries across all processors.
int MakeColMap_LL(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements...
Definition: Epetra_Export.h:62
int MakeIndicesLocal(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
int ReplaceRowMap(const Epetra_BlockMap &newmap)
Replaces the current RowMap with the user-specified map object, but only if currentmap->PointSameAs(n...
void SetAllocated(bool Flag)
long long NumGlobalBlockDiagonals64() const
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
bool NoRedundancies() const
If RemoveRedundantIndices() has been called, this query returns true, otherwise it returns false...
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
long long NumGlobalElements64() const
int ExtractGlobalRowView(int GlobalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified global row of the graph.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:63
virtual void Barrier() const =0
Epetra_Comm Barrier function.
const Epetra_BlockMap & RowMap() const
Returns the RowMap associated with this graph.
int ReplaceDomainMapAndImporter(const Epetra_BlockMap &NewDomainMap, const Epetra_Import *NewImporter)
Replaces the current DomainMap & Importer with the user-specified map object.
int Length() const
Returns length of vector.
long long * MyGlobalElements64() const
Epetra_CrsGraphData * CrsGraphData_
virtual int MyPID() const =0
Return my process ID.
int RemoveRedundantIndices()
Removes any redundant column indices in the rows of the graph.
IndexData< int > * data
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int Size(int Length_in)
Set length of a Epetra_LongLongSerialDenseVector object; init values to zero.
Epetra_CrsGraphData: The Epetra CrsGraph Data Class.
int NumMyBlockDiagonals() const
Returns the number of Block diagonal entries in the local graph, based on global row/column index com...
int RemoveGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Remove a list of elements from a specified global row of the graph.
int NumMyRows() const
Returns the number of matrix rows on this processor.
IndexData< int_type > & Data()
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
const Epetra_BlockMap & ColMap() const
Returns the Column Map associated with this graph.
bool LowerTriangular() const
If graph is lower triangular in local index space, this query returns true, otherwise it returns fals...
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
int NumMyBlockRows() const
Returns the number of block matrix rows on this processor.
long long NumGlobalPoints64() const
std::vector< EntriesInOneRow< int > > SortedEntries_
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
int InsertIndices(int Row, int NumIndices, int *Indices)
Epetra_CrsGraph(Epetra_DataAccess CV, const Epetra_BlockMap &RowMap, const int *NumIndicesPerRow, bool StaticProfile=false)
Epetra_CrsGraph constuctor with variable number of indices per row.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
long long GID64(int LID) const
void SetIndicesAreGlobal(bool Flag)
Epetra_IntSerialDenseVector NumAllocatedIndicesPerRow_
int NumMyBlockCols() const
Returns the number of Block matrix columns on this processor.
void SetIndicesAreLocal(bool Flag)
const Epetra_BlockMap & RangeMap() const
Returns the RangeMap associated with this graph.
bool IndicesAreContiguous() const
void SetNoRedundancies(bool Flag)
const Epetra_Comm * Comm_
int MaxRowDim() const
Returns the max row dimension of block entries on the processor.
int CopyAndPermuteCrsGraph(const Epetra_CrsGraph &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
int * All_Indices() const
int MakeColMap(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
void SetIndicesAreContiguous(bool Flag)
int MaxColDim() const
Returns the max column dimension of block entries on the processor.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
virtual const Epetra_BlockMap & Map() const =0
Returns a reference to the Epetra_BlockMap for this object.
int NumMyElements() const
Number of elements on the calling processor.
void epetra_shellsort(int *list, int length)
int PackAndPrepareRowMatrix(const Epetra_RowMatrix &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
int TransformToLocal()
Use FillComplete() instead.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
bool UpperTriangular() const
If graph is upper triangular in local index space, this query returns true, otherwise it returns fals...
Epetra_LongLongSerialDenseVector: A class for constructing and using dense vectors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false...
int GlobalMaxNumIndices() const
Returns the maximun number of nonzero entries across all rows across all processors.
void SetGlobalConstantsComputed(bool Flag)
bool NoDiagonal() const
If graph has no diagonal entries in global index space, this query returns true, otherwise it returns...
int RemoteIDList(int NumIDs, const int *GIDList, int *PIDList, int *LIDList) const
Returns the processor IDs and corresponding local index value for a given list of global indices...
int NumGlobalIndices(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
int Epetra_Util_binary_search(T item, const T *list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
void IncrementReferenceCount()
Increment reference count.
Definition: Epetra_Data.cpp:61
Epetra_IntSerialDenseVector All_Indices_
Epetra_SerialComm: The Epetra Serial Communication Class.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
Epetra_IntSerialDenseVector All_IndicesPlus1_
long long IndexBase64() const
void SetSorted(bool Flag)
int ReplaceColMap(const Epetra_BlockMap &newmap)
Replaces the current ColMap with the user-specified map object, but only if no entries have been inse...
long long NumGlobalNonzeros64() const
virtual int NumProc() const =0
Returns total number of processes.
int RemoveMyIndices(int LocalRow, int NumIndices, int *Indices)
Remove a list of elements from a specified local row of the graph.
const Epetra_BlockMap & DomainMap() const
Returns the DomainMap associated with this graph.
Epetra_IntSerialDenseVector IndexOffset_
bool GlobalIndicesTypeMatch(const Epetra_BlockMap &other) const
Epetra_BlockMap Map_
bool StaticProfile() const
int ExtractMyRowCopy(int LocalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified local row of the graph. Put into storage allocated by calli...
long long GCID64(int LCID_in) const
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int MaxNumIndices() const
Returns the maximum number of nonzero entries across all rows on this processor.
int MaxElementSize() const
Maximum element size across all processors.
long long GRID64(int LRID_in) const
int CopyAndPermuteRowMatrix(const Epetra_RowMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
Epetra_CombineMode
Epetra_DataAccess CV_
int NumMyEntries() const
Returns the number of entries on this processor.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
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 ExtractGlobalRowCopy(int GlobalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified global row of the graph. Put into storage allocated by call...
const Epetra_Import * Importer_
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
long long NumGlobalBlockRows64() const
Epetra_DataAccess
int * Values()
Returns pointer to the values in vector.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
long long NumGlobalCols64() const
int InsertIndicesIntoSorted(int Row, int NumIndices, int *Indices)
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
long long NumGlobalDiagonals64() const
int MakeColMap_int(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
int TAllocate(const int *numIndicesPerRow, int Inc, bool staticProfile)
int n
int Resize(int Length_in)
Resize a Epetra_LongLongSerialDenseVector object.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
virtual void Print(std::ostream &os) const
Print method.
int ** Indices() const
void Add(const long long key, const value_type value)
int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
bool GlobalConstantsComputed() const
void epetra_crsgraph_compress_out_duplicates(int len, int *list, int &newlen)
#define EPETRA_MAX(x, y)
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
Epetra_BlockMap * ReplaceCommWithSubset(const Epetra_Comm *Comm) const
Replace this BlockMap&#39;s communicator with a subset communicator.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
bool IndicesAreGlobal() const
If column indices are in global range, this query returns true, otherwise it returns false...
long long NumGlobalRows64() const
int GlobalMaxRowDim() const
Returns the max row dimension of block entries across all processors.
bool Sorted() const
If SortIndices() has been called, this query returns true, otherwise it returns false.