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] = NumIndices > 0 ? &Data.SortedEntries_[i].entries_[0]: NULL;
180  Data.SortedEntries_[i].entries_.resize(0);
181  }
182  }
183  else {
184  Data.Indices_[i] = 0;
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)
1063  ColElementSizeList = ColMap().ElementSizeList();
1064  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1065  for(int i = 0; i < numMyBlockRows; i++){
1066  int NumEntries = CrsGraphData_->NumIndicesPerRow_[i];
1067  int* indices = intData.Indices_[i];
1068  if(NumEntries > 0) {
1069  int CurNumNonzeros = 0;
1070  int RowDim = RowElementSizeList[i];
1071  for(int j = 0; j < NumEntries; j++) {
1072  int ColDim = ColElementSizeList[indices[j]];
1073  CurNumNonzeros += RowDim*ColDim;
1075  }
1077  CrsGraphData_->NumMyNonzeros_ += CurNumNonzeros;
1078  }
1079  }
1080 
1081  // Sum Up all nonzeros
1082 
1083  tempvec[0] = CrsGraphData_->NumMyEntries_;
1084  tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
1085  tempvec[2] = CrsGraphData_->NumMyDiagonals_;
1086  tempvec[3] = CrsGraphData_->NumMyNonzeros_;
1087 
1088  Comm().SumAll(&tempvec[0], &tempvec[4], 4);
1089 
1090  CrsGraphData_->NumGlobalEntries_ = tempvec[4];
1092  CrsGraphData_->NumGlobalDiagonals_ = tempvec[6];
1093  CrsGraphData_->NumGlobalNonzeros_ = tempvec[7];
1094 
1095  tempvec[0] = CrsGraphData_->MaxNumIndices_;
1096  tempvec[1] = CrsGraphData_->MaxNumNonzeros_;
1097 
1098  Comm().MaxAll(&tempvec[0], &tempvec[2], 2);
1099 
1100  CrsGraphData_->GlobalMaxNumIndices_ = (int) tempvec[2];
1101  CrsGraphData_->GlobalMaxNumNonzeros_ = (int) tempvec[3];
1102  }
1103 
1106 
1108 
1109  return(0);
1110 }
1111 
1112 void epetra_shellsort(int* list, int length)
1113 {
1114  int i, j, j2, temp, istep;
1115  unsigned step;
1116 
1117  step = length/2;
1118  while (step > 0)
1119  {
1120  for (i=step; i < length; i++)
1121  {
1122  istep = step;
1123  j = i;
1124  j2 = j-istep;
1125  temp = list[i];
1126  if (list[j2] > temp) {
1127  while ((j >= istep) && (list[j2] > temp))
1128  {
1129  list[j] = list[j2];
1130  j = j2;
1131  j2 -= istep;
1132  }
1133  list[j] = temp;
1134  }
1135  }
1136 
1137  if (step == 2)
1138  step = 1;
1139  else
1140  step = (int) (step / 2.2);
1141  }
1142 }
1143 
1144 //==============================================================================
1146  if(IndicesAreGlobal())
1147  EPETRA_CHK_ERR(-1);
1148  if(Sorted())
1149  return(0);
1150 
1151  // For each row, sort column entries from smallest to largest.
1152  // Use shell sort, which is fast if indices are already sorted.
1153 
1154  const int numMyBlockRows = NumMyBlockRows();
1155  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1156  for(int i = 0; i < numMyBlockRows; i++){
1157  int n = CrsGraphData_->NumIndicesPerRow_[i];
1158  int* const list = intData.Indices_[i];
1159 
1160  epetra_shellsort(list, n);
1161 // int m = n/2;
1162 
1163 // while(m > 0) {
1164  // int max = n - m;
1165 // for(int j = 0; j < max; j++) {
1166 // int k = j;
1167 // while(k>-1) {
1168 // if(list[k+m] >= list[k])
1169 // break;
1170 // int itemp = list[k+m];
1171 // list[k+m] = list[k];
1172 // list[k] = itemp;
1173 // k-=m;
1174 // }
1175 // }
1176 // m = m/2;
1177 // }
1178  }
1179  SetSorted(true);
1180 
1181  if(CrsGraphData_->ReferenceCount() > 1)
1182  return(1);
1183  else
1184  return(0);
1185 }
1186 
1187 void epetra_crsgraph_compress_out_duplicates(int len, int* list, int& newlen)
1188 {
1189  //
1190  //This function runs the array ('list') checking for
1191  //duplicate entries. Any duplicates that are found are
1192  //removed by sliding subsequent data down in the array,
1193  //over-writing the duplicates. Finally, the new length
1194  //of the array (i.e., the number of unique entries) is
1195  //placed in the output parameter 'newlen'. The array is
1196  //**not** re-allocated.
1197  //
1199  //Important assumption: The array contents are assumed to
1200  //be sorted before this function is called. If the array
1201  //contents are not sorted, then the behavior of this
1202  //function is undefined.
1204  //
1205 
1206  if (len < 2) return;
1207 
1208  int* ptr0 = &list[0];
1209  int* ptr1 = &list[1];
1210 
1211  int* ptr_end = &list[len-1];
1212 
1213  while(*ptr0 != *ptr1 && ptr1 < ptr_end) {
1214  ++ptr0;
1215  ++ptr1;
1216  }
1217 
1218  if (ptr1 < ptr_end) {
1219  //if ptr1 < ptr_end we've found a duplicate...
1220 
1221  ++ptr0;
1222  ++ptr1;
1223 
1224  while(*ptr0 == *ptr1 && ptr1 < ptr_end) ++ptr1;
1225 
1226  while(ptr1 < ptr_end) {
1227 
1228  int val = *ptr1++;
1229 
1230  while(val == *ptr1 && ptr1 < ptr_end) {
1231  ++ptr1;
1232  }
1233 
1234  *ptr0++ = val;
1235  }
1236 
1237  if (*(ptr0-1) != *ptr1) *ptr0++ = *ptr1;
1238 
1239  int num_removed = (int)(ptr_end - ptr0 + 1);
1240  newlen = len - num_removed;
1241  }
1242  else {
1243  if (*ptr0 == *ptr1) newlen = len - 1;
1244  else newlen = len;
1245  }
1246 }
1247 
1248 //==============================================================================
1250 {
1251  if(!Sorted())
1252  EPETRA_CHK_ERR(-1); // Must have sorted index set
1253  if(IndicesAreGlobal())
1254  EPETRA_CHK_ERR(-2); // Indices must be local
1255 
1256  // Note: This function assumes that SortIndices was already called.
1257  // For each row, remove column indices that are repeated.
1258 
1259  const int numMyBlockRows = NumMyBlockRows();
1260  bool found_redundancies = false;
1261 
1262  if(NoRedundancies() == false) {
1263  int* numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
1265  for(int i=0; i<numMyBlockRows; ++i) {
1266  int NumIndices = numIndicesPerRow[i];
1267  int* col_indices = this->Indices(i);
1268 
1269  if(NumIndices > 1) {
1270  epetra_crsgraph_compress_out_duplicates(NumIndices, col_indices,
1271  numIndicesPerRow[i]);
1272  }
1273  if (NumIndices != numIndicesPerRow[i]) {
1274  found_redundancies = true;
1275  }
1276  }
1277  if (found_redundancies && !CrsGraphData_->StaticProfile_)
1278  {
1279  for(int i=0; i<numMyBlockRows; ++i) {
1280  int* col_indices = this->Indices(i);
1281 
1282  // update vector size and address in memory
1283  intData.SortedEntries_[i].entries_.assign(col_indices, col_indices+numIndicesPerRow[i]);
1284  if (numIndicesPerRow[i] > 0) {
1285  intData.Indices_[i] = &(intData.SortedEntries_[i].entries_[0]);
1286  }
1287  else {
1288  intData.Indices_[i] = NULL;
1289  }
1290  }
1291  }
1292  }
1293 
1294  SetNoRedundancies(true);
1295  return 0;
1296 }
1297 
1298 //==============================================================================
1300 {
1301  // determine if graph is upper or lower triangular or has no diagonal
1302 
1303  if(!Sorted())
1304  EPETRA_CHK_ERR(-1); // Must have sorted index set
1305  if(IndicesAreGlobal())
1306  EPETRA_CHK_ERR(-2); // Indices must be local
1307 
1310 
1311  const Epetra_BlockMap& rowMap = RowMap();
1312  const Epetra_BlockMap& colMap = ColMap();
1313 
1314  const int numMyBlockRows = NumMyBlockRows();
1315 
1316  for(int i = 0; i < numMyBlockRows; i++) {
1317  int NumIndices = NumMyIndices(i);
1318  if(NumIndices > 0) {
1319 #if defined(EPETRA_NO_64BIT_GLOBAL_INDICES) && !defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
1320  int ig = rowMap.GID(i);
1321 #else
1322  long long ig = rowMap.GID64(i);
1323 #endif
1324  int* col_indices = this->Indices(i);
1325 
1326  int jl_0 = col_indices[0];
1327  int jl_n = col_indices[NumIndices-1];
1328 
1329  if(jl_n > i) CrsGraphData_->LowerTriangular_ = false;
1330  if(jl_0 < i) CrsGraphData_->UpperTriangular_ = false;
1331 
1332  //jl will be the local-index for the diagonal that we
1333  //want to search for.
1334  int jl = colMap.LID(ig);
1335 
1336  int insertPoint = -1;
1337  if (Epetra_Util_binary_search(jl, col_indices, NumIndices, insertPoint)>-1) {
1340  }
1341  }
1342  }
1343 
1345 
1346  if(CrsGraphData_->ReferenceCount() > 1)
1347  return(1);
1348  else
1349  return(0);
1350 }
1351 
1352 // private =====================================================================
1353 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1355  const Epetra_BlockMap& rangeMap)
1356 {
1357  (void)rangeMap;
1358  int i;
1359  int j;
1360 
1362  return(0); // Already have a Column Map
1363 
1364  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1365  if(IndicesAreLocal())
1366  EPETRA_CHK_ERR(-1); // Return error: Indices must be global
1367 
1368  // Scan all column indices and sort into two groups:
1369  // Local: those whose GID matches a GID of the domain map on this processor and
1370  // Remote: All others.
1371  int numDomainElements = domainMap.NumMyElements();
1372  bool * LocalGIDs = 0;
1373  if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
1374  for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
1375 
1376  // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
1377  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
1378  // and the number of block rows.
1379  const int numMyBlockRows = NumMyBlockRows();
1380  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
1381  //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
1382  Epetra_HashTable<int> RemoteGIDs(hashsize);
1383  Epetra_HashTable<int> RemoteGIDList(hashsize);
1384 
1385  int NumLocalColGIDs = 0;
1386  int NumRemoteColGIDs = 0;
1387  const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
1388  for(i = 0; i < numMyBlockRows; i++) {
1389  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1390  int* ColIndices = intData.Indices_[i];
1391  for(j = 0; j < NumIndices; j++) {
1392  int GID = ColIndices[j];
1393  // Check if GID matches a row GID
1394  int LID = domainMap.LID(GID);
1395  if(LID != -1) {
1396  bool alreadyFound = LocalGIDs[LID];
1397  if (!alreadyFound) {
1398  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1399  NumLocalColGIDs++;
1400  }
1401  }
1402  else {
1403  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1404  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1405  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1406  }
1407  }
1408  }
1409  }
1410 
1411  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
1412  if (domainMap.Comm().NumProc()==1) {
1413 
1414  if (NumRemoteColGIDs!=0) {
1415  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
1416  }
1417  if (NumLocalColGIDs==numDomainElements) {
1418  CrsGraphData_->ColMap_ = domainMap;
1419  CrsGraphData_->HaveColMap_ = true;
1420  if (LocalGIDs!=0) delete [] LocalGIDs;
1421  return(0);
1422  }
1423  }
1424 
1425  // Now build integer array containing column GIDs
1426  // Build back end, containing remote GIDs, first
1427  int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1428  Epetra_IntSerialDenseVector ColIndices;
1429  if(numMyBlockCols > 0)
1430  ColIndices.Size(numMyBlockCols);
1431 
1432  int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
1433 
1434  for(i = 0; i < NumRemoteColGIDs; i++)
1435  RemoteColIndices[i] = RemoteGIDList.Get(i);
1436 
1437  int NLists = 1;
1439  Epetra_IntSerialDenseVector SizeList;
1440  int* RemoteSizeList = 0;
1441  bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
1442 
1443  if(NumRemoteColGIDs > 0)
1444  PIDList.Size(NumRemoteColGIDs);
1445 
1446  if(DoSizes) {
1447  if(numMyBlockCols > 0)
1448  SizeList.Size(numMyBlockCols);
1449  RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
1450  NLists++;
1451  }
1452  EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
1453 
1454  // Sort External column indices so that all columns coming from a given remote processor are contiguous
1455 
1456  Epetra_Util Util;
1457  int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
1458  SortLists[0] = RemoteColIndices;
1459  SortLists[1] = RemoteSizeList;
1460  Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists, 0, 0);
1462  // Sort external column indices so that columns from a given remote processor are not only contiguous
1463  // but also in ascending order. NOTE: I don't know if the number of externals associated
1464  // with a given remote processor is known at this point ... so I count them here.
1465 
1466  NLists--;
1467  int StartCurrent, StartNext;
1468  StartCurrent = 0; StartNext = 1;
1469  while ( StartNext < NumRemoteColGIDs ) {
1470  if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
1471  else {
1472  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1473  Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
1474  StartCurrent = StartNext; StartNext++;
1475  }
1476  }
1477  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1478  Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
1479  }
1480 
1481  // Now fill front end. Two cases:
1482  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
1483  // can simply read the domain GIDs into the front part of ColIndices, otherwise
1484  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
1485  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
1486 
1487  if(NumLocalColGIDs == domainMap.NumMyElements()) {
1488  domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
1489  if(DoSizes)
1490  domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
1491  }
1492  else {
1493  int NumMyElements = domainMap.NumMyElements();
1494  int* MyGlobalElements = domainMap.MyGlobalElements();
1495  int* ElementSizeList = 0;
1496  if(DoSizes)
1497  ElementSizeList = domainMap.ElementSizeList();
1498  int NumLocalAgain = 0;
1499  for(i = 0; i < NumMyElements; i++) {
1500  if(LocalGIDs[i]) {
1501  if(DoSizes)
1502  SizeList[NumLocalAgain] = ElementSizeList[i];
1503  ColIndices[NumLocalAgain++] = MyGlobalElements[i];
1504  }
1505  }
1506  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
1507  }
1508 
1509  // Done with this array
1510  if (LocalGIDs!=0) delete [] LocalGIDs;
1511 
1512 
1513  // Make Column map with same element sizes as Domain map
1514 
1515  if(domainMap.MaxElementSize() == 1) { // Simple map
1516  Epetra_Map temp((int) -1, numMyBlockCols, ColIndices.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
1517  CrsGraphData_->ColMap_ = temp;
1518  }
1519  else if(domainMap.ConstantElementSize()) { // Constant Block size map
1520  Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(), (int) domainMap.IndexBase64(), domainMap.Comm());
1521  CrsGraphData_->ColMap_ = temp;
1522  }
1523  else { // Most general case where block size is variable.
1524  Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
1525  CrsGraphData_->ColMap_ = temp;
1526  }
1527  CrsGraphData_->HaveColMap_ = true;
1528 
1529  return(0);
1530 }
1531 #endif
1532 //==============================================================================
1533 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1535  const Epetra_BlockMap& rangeMap)
1536 {
1537  (void)rangeMap;
1538  int i;
1539  int j;
1540 
1542  return(0); // Already have a Column Map
1543 
1544  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1545  if(IndicesAreLocal())
1546  EPETRA_CHK_ERR(-1); // Return error: Indices must be global
1547 
1548  // Scan all column indices and sort into two groups:
1549  // Local: those whose GID matches a GID of the domain map on this processor and
1550  // Remote: All others.
1551  int numDomainElements = domainMap.NumMyElements();
1552  bool * LocalGIDs = 0;
1553  if (numDomainElements>0) LocalGIDs = new bool[numDomainElements];
1554  for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
1555 
1556  // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
1557  // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
1558  // and the number of block rows.
1559  const int numMyBlockRows = NumMyBlockRows();
1560  int hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
1561  //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
1562  Epetra_HashTable<int> RemoteGIDs(hashsize);
1563  Epetra_HashTable<long long> RemoteGIDList(hashsize);
1564 
1565  int NumLocalColGIDs = 0;
1566  int NumRemoteColGIDs = 0;
1567 
1568  if(IndicesAreLocal())
1569  {
1571 
1572  for(i = 0; i < numMyBlockRows; i++) {
1573  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1574  int* ColIndices = intData.Indices_[i];
1575  for(j = 0; j < NumIndices; j++) {
1576  int GID = ColIndices[j];
1577  // Check if GID matches a row GID
1578  int LID = domainMap.LID(GID);
1579  if(LID != -1) {
1580  bool alreadyFound = LocalGIDs[LID];
1581  if (!alreadyFound) {
1582  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1583  NumLocalColGIDs++;
1584  }
1585  }
1586  else {
1587  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1588  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1589  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1590  }
1591  }
1592  }
1593  }
1594  }
1595  else if(IndicesAreGlobal())
1596  {
1598 
1599  for(i = 0; i < numMyBlockRows; i++) {
1600  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1601  long long* ColIndices = LLData.Indices_[i];
1602  for(j = 0; j < NumIndices; j++) {
1603  long long GID = ColIndices[j];
1604  // Check if GID matches a row GID
1605  int LID = domainMap.LID(GID);
1606  if(LID != -1) {
1607  bool alreadyFound = LocalGIDs[LID];
1608  if (!alreadyFound) {
1609  LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
1610  NumLocalColGIDs++;
1611  }
1612  }
1613  else {
1614  if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
1615  RemoteGIDs.Add(GID, NumRemoteColGIDs);
1616  RemoteGIDList.Add(NumRemoteColGIDs++, GID);
1617  }
1618  }
1619  }
1620  }
1621  }
1622 
1623  // Possible short-circuit: If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
1624  if (domainMap.Comm().NumProc()==1) {
1625 
1626  if (NumRemoteColGIDs!=0) {
1627  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
1628  }
1629  if (NumLocalColGIDs==numDomainElements) {
1630  CrsGraphData_->ColMap_ = domainMap;
1631  CrsGraphData_->HaveColMap_ = true;
1632  if (LocalGIDs!=0) delete [] LocalGIDs;
1633  return(0);
1634  }
1635  }
1636 
1637  // Now build integer array containing column GIDs
1638  // Build back end, containing remote GIDs, first
1639  int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
1641  if(numMyBlockCols > 0)
1642  ColIndices.Size(numMyBlockCols);
1643 
1644  long long* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
1645 
1646  for(i = 0; i < NumRemoteColGIDs; i++)
1647  RemoteColIndices[i] = RemoteGIDList.Get(i);
1648 
1649  int NLists = 0;
1651  Epetra_IntSerialDenseVector SizeList;
1652  int* RemoteSizeList = 0;
1653  bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
1654 
1655  if(NumRemoteColGIDs > 0)
1656  PIDList.Size(NumRemoteColGIDs);
1657 
1658  if(DoSizes) {
1659  if(numMyBlockCols > 0)
1660  SizeList.Size(numMyBlockCols);
1661  RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
1662  NLists++;
1663  }
1664  EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
1665 
1666  // Sort External column indices so that all columns coming from a given remote processor are contiguous
1667 
1668  Epetra_Util Util;
1669  //int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
1670  //SortLists[0] = RemoteColIndices;
1671  //SortLists[1] = RemoteSizeList;
1672  Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, &RemoteSizeList, 1, &RemoteColIndices);
1674  // Sort external column indices so that columns from a given remote processor are not only contiguous
1675  // but also in ascending order. NOTE: I don't know if the number of externals associated
1676  // with a given remote processor is known at this point ... so I count them here.
1677 
1678  int* SortLists[1] = {0};
1679 
1680  int StartCurrent, StartNext;
1681  StartCurrent = 0; StartNext = 1;
1682  while ( StartNext < NumRemoteColGIDs ) {
1683  if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
1684  else {
1685  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1686  Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
1687  StartCurrent = StartNext; StartNext++;
1688  }
1689  }
1690  if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
1691  Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
1692  }
1693 
1694  // Now fill front end. Two cases:
1695  // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
1696  // can simply read the domain GIDs into the front part of ColIndices, otherwise
1697  // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
1698  // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
1699 
1700  if(NumLocalColGIDs == domainMap.NumMyElements()) {
1701  domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
1702  if(DoSizes)
1703  domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
1704  }
1705  else {
1706  int NumMyElements = domainMap.NumMyElements();
1707  long long* MyGlobalElements = domainMap.MyGlobalElements64();
1708  int* ElementSizeList = 0;
1709  if(DoSizes)
1710  ElementSizeList = domainMap.ElementSizeList();
1711  int NumLocalAgain = 0;
1712  for(i = 0; i < NumMyElements; i++) {
1713  if(LocalGIDs[i]) {
1714  if(DoSizes)
1715  SizeList[NumLocalAgain] = ElementSizeList[i];
1716  ColIndices[NumLocalAgain++] = MyGlobalElements[i];
1717  }
1718  }
1719  assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
1720  }
1721 
1722  // Done with this array
1723  if (LocalGIDs!=0) delete [] LocalGIDs;
1724 
1725 
1726  // Make Column map with same element sizes as Domain map
1727 
1728  if(domainMap.MaxElementSize() == 1) { // Simple map
1729  Epetra_Map temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.IndexBase64(), domainMap.Comm());
1730  CrsGraphData_->ColMap_ = temp;
1731  }
1732  else if(domainMap.ConstantElementSize()) { // Constant Block size map
1733  Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(),domainMap.IndexBase64(), domainMap.Comm());
1734  CrsGraphData_->ColMap_ = temp;
1735  }
1736  else { // Most general case where block size is variable.
1737  Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), domainMap.IndexBase64(), domainMap.Comm());
1738  CrsGraphData_->ColMap_ = temp;
1739  }
1740  CrsGraphData_->HaveColMap_ = true;
1741 
1742  return(0);
1743 }
1744 #endif
1745 
1747  const Epetra_BlockMap& rangeMap)
1748 {
1749  if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
1750  throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for domainMap and rangeMap", -1);
1751 
1752  if(!RowMap().GlobalIndicesTypeMatch(domainMap))
1753  throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for row map and incoming rangeMap", -1);
1754 
1755  if(RowMap().GlobalIndicesInt())
1756 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1757  return MakeColMap_int(domainMap, rangeMap);
1758 #else
1759  throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesInt but no API for it.",-1);
1760 #endif
1761 
1762  if(RowMap().GlobalIndicesLongLong())
1763 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1764  return MakeColMap_LL(domainMap, rangeMap);
1765 #else
1766  throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesLongLong but no API for it.",-1);
1767 #endif
1768 
1769  throw ReportError("Epetra_CrsGraph::MakeColMap: Internal error, unable to determine global index type of maps", -1);
1770 }
1771 
1772 // protected ===================================================================
1774  if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
1775  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for domainMap and rangeMap", -1);
1776 
1777  if(!RowMap().GlobalIndicesTypeMatch(domainMap))
1778  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for row map and incoming rangeMap", -1);
1779 
1780  ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
1782  EPETRA_CHK_ERR(-1); // Return error: Indices must not be both local and global
1783 
1784  MakeColMap(domainMap, rangeMap); // If user has not prescribed column map, create one from indices
1785  const Epetra_BlockMap& colmap = ColMap();
1786 
1787  // Store number of local columns
1790 
1791  // Transform indices to local index space
1792  const int numMyBlockRows = NumMyBlockRows();
1793 
1794  if(IndicesAreGlobal()) {
1795  // Check if ColMap is monotone. If not, the list will get unsorted.
1796  bool mapMonotone = true;
1797  {
1798  long long oldGID = colmap.GID64(0);
1799  for (int i=1; i<colmap.NumMyElements(); ++i) {
1800  if (oldGID > colmap.GID64(i)) {
1801  mapMonotone = false;
1802  break;
1803  }
1804  oldGID = colmap.GID64(i);
1805  }
1806  }
1807  if (Sorted())
1808  SetSorted(mapMonotone);
1809 
1810  // We don't call Data<int>() here because that would not work (it will throw)
1811  // if GlobalIndicesLongLong() and IndicesAreGlobal(). This is because
1812  // the local flag is not set yet. We are in the middle of the transaction here.
1813  // In all other cases, one should call the function Data<int> or Data<long long>
1814  // instead of obtaining the pointer directly.
1816 
1817  if(RowMap().GlobalIndicesInt())
1818  {
1819  // now comes the actual transformation
1820  for(int i = 0; i < numMyBlockRows; i++) {
1821  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1822  int* ColIndices = intData.Indices_[i];
1823  for(int j = 0; j < NumIndices; j++) {
1824  int GID = ColIndices[j];
1825  int LID = colmap.LID(GID);
1826  if(LID != -1)
1827  ColIndices[j] = LID;
1828  else
1829  throw ReportError("Internal error in FillComplete ",-1);
1830  }
1831  }
1832  }
1833  else if(RowMap().GlobalIndicesLongLong())
1834  {
1835 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1837 
1838  if (!CrsGraphData_->StaticProfile_) {
1839  // Use the resize trick used in TAllocate.
1840  const long long indexBaseMinusOne = IndexBase64() - 1;
1841  for(int i = 0; i < numMyBlockRows; i++) {
1842  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1843  intData.SortedEntries_[i].entries_.resize(NumIndices, indexBaseMinusOne);
1844  intData.Indices_[i] = NumIndices > 0 ? &intData.SortedEntries_[i].entries_[0]: NULL;
1845  intData.SortedEntries_[i].entries_.resize(0);
1846  }
1847  }
1848 
1849  // now comes the actual transformation
1850  for(int i = 0; i < numMyBlockRows; i++) {
1851  const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
1852  long long* ColIndices = LL_Data.Indices_[i];
1853  int* intColIndices = intData.Indices_[i];
1854  for(int j = 0; j < NumIndices; j++) {
1855  long long GID = ColIndices[j];
1856  int LID = colmap.LID(GID);
1857  if(LID != -1)
1858  intColIndices[j] = LID;
1859  else
1860  throw ReportError("Internal error in FillComplete ",-1);
1861  }
1862  }
1863 
1864  LL_Data.Deallocate(); // deallocate long long data since indices are local now.
1865 #else
1866  throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: GlobalIndicesLongLong but no long long API", -1);
1867 #endif
1868  }
1869 
1870  }
1871 
1872  SetIndicesAreLocal(true);
1873  SetIndicesAreGlobal(false);
1874 
1875  if(CrsGraphData_->ReferenceCount() > 1)
1876  return(1); // return 1 if data is shared
1877  else
1878  return(0);
1879 }
1880 
1881 //==============================================================================
1883  int NumIndices;
1884  const int numMyBlockRows = NumMyBlockRows();
1885 
1887 
1888  if(StorageOptimized())
1889  return(0); // Have we been here before?
1890  if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
1891 
1892  bool Contiguous = true; // Assume contiguous is true
1893  for(int i = 1; i < numMyBlockRows; i++) {
1894  NumIndices = CrsGraphData_->NumIndicesPerRow_[i-1];
1895  int NumAllocateIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[i-1];
1896 
1897  // Check if NumIndices is same as NumAllocatedIndices and
1898  // check if end of beginning of current row starts immediately after end of previous row.
1899  if((NumIndices != NumAllocateIndices) ||
1900  (Data.Indices_[i] != Data.Indices_[i-1] + NumIndices)) {
1901  Contiguous = false;
1902  break;
1903  }
1904  }
1905 
1906  // NOTE: At the end of the above loop set, there is a possibility that NumIndices and NumAllocatedIndices
1907  // for the last row could be different, but I don't think it matters.
1908 
1909 
1910  if((CrsGraphData_->CV_ == View) && !Contiguous)
1911  return(3); // This is user data, it's not contiguous and we can't make it so.
1912 
1913  // This next step constructs the scan sum of the number of indices per row. Note that the
1914  // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
1915  // careful with this sum operation
1916 
1919 
1920  int * numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
1921  int curNumIndices = numIndicesPerRow[0];
1922  numIndicesPerRow[0] = 0;
1923  for (int i=0; i<numMyBlockRows; ++i) {
1924  int nextNumIndices = numIndicesPerRow[i+1];
1925  numIndicesPerRow[i+1] = numIndicesPerRow[i]+curNumIndices;
1926  curNumIndices = nextNumIndices;
1927  }
1928 
1929 // *******************************
1930 // Data NOT contiguous, make it so
1931 // *******************************
1932  if(!Contiguous) { // Must pack indices if not already contiguous
1933 
1934  // Allocate one big integer array for all index values
1935  if (!(StaticProfile())) { // If static profile, All_Indices_ is already allocated, only need to pack data
1936  int errorcode = Data.All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
1937  if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
1938  }
1939  // Pack indices into All_Indices_
1940 
1941  int* all_indices = Data.All_Indices_.Values();
1942  int * indexOffset = CrsGraphData_->IndexOffset_.Values();
1943  int ** indices = Data.Indices_;
1944 
1945  if (!(StaticProfile())) {
1946 #ifdef EPETRA_HAVE_OMP
1947 #pragma omp parallel for default(none) shared(indexOffset,all_indices,indices)
1948 #endif
1949  for(int i = 0; i < numMyBlockRows; i++) {
1950  int numColIndices = indexOffset[i+1] - indexOffset[i];
1951  int* ColIndices = indices[i];
1952  int *newColIndices = all_indices+indexOffset[i];
1953  for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
1954  }
1955  for(int i = 0; i < numMyBlockRows; i++) {
1956  if (indices[i]!=0) {
1957  Data.SortedEntries_[i].entries_.clear();
1958  indices[i] = 0;
1959  }
1960  }
1961  } // End of non-contiguous non-static section
1962  else {
1963 
1964  for(int i = 0; i < numMyBlockRows; i++) {
1965  int numColIndices = indexOffset[i+1] - indexOffset[i];
1966  int* ColIndices = indices[i];
1967  int *newColIndices = all_indices+indexOffset[i];
1968  if (ColIndices!=newColIndices) // No need to copy if pointing to same space
1969  for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
1970  indices[i] = 0;
1971  }
1972  } // End of non-Contiguous static section
1973  } // End of non-Contiguous section
1974  else { // Start of Contiguous section
1975  // if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
1976  // Execute the assignment block in parallel using the same pattern as SpMV
1977  // in order to improve page placement
1978  if (numMyBlockRows > 0 && !(StaticProfile())) {
1979  const int numMyNonzeros = NumMyNonzeros();
1980  int errorcode = Data.All_Indices_.Size(numMyNonzeros);
1981  if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
1982  int* new_all_indices = Data.All_Indices_.Values();
1983  int* old_all_indices = Data.Indices_[0];
1984  int * indexOffset = CrsGraphData_->IndexOffset_.Values();
1985 
1986 #ifdef EPETRA_HAVE_OMP
1987 #pragma omp parallel for default(none) shared(indexOffset,old_all_indices,new_all_indices)
1988 #endif
1989  for(int i = 0; i < numMyBlockRows; i++) {
1990  int numColIndices = indexOffset[i+1] - indexOffset[i];
1991  int *oldColIndices = old_all_indices+indexOffset[i];
1992  int *newColIndices = new_all_indices+indexOffset[i];
1993  for(int j = 0; j < numColIndices; j++) newColIndices[j] = oldColIndices[j];
1994  }
1995 
1996 //#ifdef EPETRA_HAVE_OMP
1997 //#pragma omp parallel for default(none) shared(all_indices_values,indices_values)
1998 //#endif
1999 // for(int ii=0; ii<numMyNonzeros; ++ii) {
2000 // all_indices_values[ii] = indices_values[ii];
2001 // }
2002  }
2003  }
2004 
2005  // Delete unneeded storage
2007  delete [] Data.Indices_; Data.Indices_=0;
2008  Data.SortedEntries_.clear();
2009 
2010  SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
2012 
2013 /*
2014 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2015  All_IndicesPlus1(); // see if preemptively calling this improves Multiply timing.
2016 #endif
2017 */
2018 
2019  return(0);
2020 }
2021 
2022 //==============================================================================
2023 template<typename int_type>
2024 int Epetra_CrsGraph::ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
2025 {
2026  int j;
2027 
2028  int locRow = LRID(Row); // Normalize row range
2029 
2030  if(locRow < 0 || locRow >= NumMyBlockRows())
2031  EPETRA_CHK_ERR(-1); // Not in Row range
2032 
2033  NumIndices = NumMyIndices(locRow);
2034  if(LenOfIndices < NumIndices)
2035  EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
2036 
2037  if(IndicesAreLocal())
2038  {
2039  int * srcIndices = TIndices<int>(locRow);
2040  // static_cast is ok because global indices were created from int values and hence must fit ints
2041  for(j = 0; j < NumIndices; j++)
2042  targIndices[j] = static_cast<int_type>(GCID64(srcIndices[j]));
2043  }
2044  else
2045  {
2046  int_type * srcIndices = TIndices<int_type>(locRow);
2047  for(j = 0; j < NumIndices; j++)
2048  targIndices[j] = srcIndices[j];
2049  }
2050 
2051  return(0);
2052 }
2053 
2054 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2055 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
2056 {
2057  if(RowMap().GlobalIndicesInt())
2058  return ExtractGlobalRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
2059  else
2060  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy int version called for a graph that is not int.", -1);
2061 }
2062 #endif
2063 
2064 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2065 int Epetra_CrsGraph::ExtractGlobalRowCopy(long long Row, int LenOfIndices, int& NumIndices, long long* targIndices) const
2066 {
2067  if(RowMap().GlobalIndicesLongLong())
2068  return ExtractGlobalRowCopy<long long>(Row, LenOfIndices, NumIndices, targIndices);
2069  else
2070  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy long long version called for a graph that is not long long.", -1);
2071 }
2072 #endif
2073 
2074 //==============================================================================
2075 template<typename int_type>
2076 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
2077 {
2078  int j;
2079 
2080  if(Row < 0 || Row >= NumMyBlockRows())
2081  EPETRA_CHK_ERR(-1); // Not in Row range
2082 
2083  NumIndices = NumMyIndices(Row);
2084  if(LenOfIndices < NumIndices)
2085  EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
2086 
2087  if(IndicesAreGlobal())
2088  EPETRA_CHK_ERR(-3); // There are no local indices yet
2089 
2090  int * srcIndices = TIndices<int>(Row);
2091  for(j = 0; j < NumIndices; j++)
2092  targIndices[j] = srcIndices[j];
2093 
2094  return(0);
2095 }
2096 
2097 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
2098 {
2099  if(RowMap().GlobalIndicesTypeValid())
2100  return ExtractMyRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
2101  else
2102  throw ReportError("Epetra_CrsGraph::ExtractMyRowCopy graph global index type unknown.", -1);
2103 }
2104 
2105 //==============================================================================
2106 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2107 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const
2108 {
2109  if(!RowMap().GlobalIndicesInt())
2110  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView int version called for a graph that is not int.", -1);
2111 
2112  int locRow = LRID(Row); // Normalize row range
2113 
2114  if(locRow < 0 || locRow >= NumMyBlockRows())
2115  EPETRA_CHK_ERR(-1); // Not in Row range
2116 
2117  if(IndicesAreLocal())
2118  EPETRA_CHK_ERR(-2); // There are no global indices
2119 
2120  NumIndices = NumMyIndices(locRow);
2121 
2122  targIndices = TIndices<int>(locRow);
2123 
2124  return(0);
2125 }
2126 #endif
2127 //==============================================================================
2128 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2129 int Epetra_CrsGraph::ExtractGlobalRowView(long long Row, int& NumIndices, long long*& targIndices) const
2130 {
2131  if(!RowMap().GlobalIndicesLongLong())
2132  throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView long long version called for a graph that is not long long.", -1);
2133 
2134  int locRow = LRID(Row); // Normalize row range
2135 
2136  if(locRow < 0 || locRow >= NumMyBlockRows())
2137  EPETRA_CHK_ERR(-1); // Not in Row range
2138 
2139  if(IndicesAreLocal())
2140  EPETRA_CHK_ERR(-2); // There are no global indices
2141 
2142  NumIndices = NumMyIndices(locRow);
2143 
2144  targIndices = TIndices<long long>(locRow);
2145 
2146  return(0);
2147 }
2148 #endif
2149 //==============================================================================
2150 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const
2151 {
2152  if(Row < 0 || Row >= NumMyBlockRows())
2153  EPETRA_CHK_ERR(-1); // Not in Row range
2154 
2155  if(IndicesAreGlobal())
2156  EPETRA_CHK_ERR(-2); // There are no local indices
2157 
2158  NumIndices = NumMyIndices(Row);
2159 
2160  targIndices = TIndices<int>(Row);
2161 
2162  return(0);
2163 }
2164 
2165 //==============================================================================
2166 int Epetra_CrsGraph::NumGlobalIndices(long long Row) const {
2167 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
2168  int locRow = LRID((int) Row);
2169 #else
2170  int locRow = LRID(Row);
2171 #endif
2172  if(locRow != -1)
2173  return(NumMyIndices(locRow));
2174  else
2175  return(0); // No indices for this row on this processor
2176 }
2177 
2178 //==============================================================================
2180 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
2181  int locRow = LRID((int) Row);
2182 #else
2183  int locRow = LRID(Row);
2184 #endif
2185  if(locRow != -1)
2186  return(NumAllocatedMyIndices(locRow));
2187  else
2188  return(0); // No indices allocated for this row on this processor
2189 }
2190 
2191 //==============================================================================
2193 {
2194  if (RowMap().PointSameAs(newmap)) {
2195  Epetra_DistObject::Map_ = newmap;
2196  CrsGraphData_->RowMap_ = newmap;
2198  return(0);
2199  }
2200 
2201  return(-1);
2202 }
2203 
2204 //==============================================================================
2206 {
2207  if (!HaveColMap() && !IndicesAreLocal() && !IndicesAreGlobal() && newmap.GlobalIndicesTypeMatch(RowMap())) {
2208  CrsGraphData_->ColMap_ = newmap;
2214  CrsGraphData_->NumMyCols_ = newmap.NumMyPoints();
2215  CrsGraphData_->HaveColMap_ = true;
2216  return(0);
2217  }
2218 
2219  if(ColMap().PointSameAs(newmap)) {
2220  CrsGraphData_->ColMap_ = newmap;
2222  return(0);
2223  }
2224 
2225  return(-1);
2226 }
2227 
2228 
2229 //==============================================================================
2230 int Epetra_CrsGraph::ReplaceDomainMapAndImporter(const Epetra_BlockMap& NewDomainMap, const Epetra_Import * NewImporter) {
2231  int rv=0;
2232  if( !NewImporter && ColMap().SameAs(NewDomainMap)) {
2233  CrsGraphData_->DomainMap_ = NewDomainMap;
2234  delete CrsGraphData_->Importer_;
2235  CrsGraphData_->Importer_ = 0;
2236  }
2237  else if(NewImporter && ColMap().SameAs(NewImporter->TargetMap()) && NewDomainMap.SameAs(NewImporter->SourceMap())) {
2238  CrsGraphData_->DomainMap_ = NewDomainMap;
2239  delete CrsGraphData_->Importer_;
2240  CrsGraphData_->Importer_ = new Epetra_Import(*NewImporter);
2241  }
2242  else
2243  rv=-1;
2244  return rv;
2245 }
2246 
2247 //==============================================================================
2249  const Epetra_BlockMap *newDomainMap=0, *newRangeMap=0, *newColMap=0;
2250  Epetra_Import * newImport=0;
2251  Epetra_Export * newExport=0;
2252 
2253  const Epetra_Comm *newComm = newMap ? &newMap->Comm() : 0;
2254 
2255  if(DomainMap().SameAs(RowMap())) {
2256  // Common case: original domain and row Maps are identical.
2257  // In that case, we need only replace the original domain Map
2258  // with the new Map. This ensures that the new domain and row
2259  // Maps _stay_ identical.
2260  newDomainMap = newMap;
2261  }
2262  else
2263  newDomainMap = DomainMap().ReplaceCommWithSubset(newComm);
2264 
2265  if(RangeMap().SameAs(RowMap())){
2266  // Common case: original range and row Maps are identical. In
2267  // that case, we need only replace the original range Map with
2268  // the new Map. This ensures that the new range and row Maps
2269  // _stay_ identical.
2270  newRangeMap = newMap;
2271  }
2272  else
2273  newRangeMap = RangeMap().ReplaceCommWithSubset(newComm);
2274 
2275  newColMap=ColMap().ReplaceCommWithSubset(newComm);
2276 
2277  if(newComm) {
2278  // (Re)create the Export and / or Import if necessary.
2279  //
2280  // The operations below are collective on the new communicator.
2281  //
2282  if(RangeMap().DataPtr() != RowMap().DataPtr())
2283  newExport = new Epetra_Export(*newMap,*newRangeMap);
2284 
2285  if(DomainMap().DataPtr() != ColMap().DataPtr())
2286  newImport = new Epetra_Import(*newColMap,*newDomainMap);
2287  }
2288 
2289  // CrsGraphData things
2290  if(CrsGraphData_->ReferenceCount() !=1)
2291  throw ReportError("Epetra_CrsGraph::RemoveEmptyProcessesInPlace does not work for shared CrsGraphData_",-2);
2292 
2293  // Dummy map for the non-active procs
2294 #if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
2295  Epetra_Map dummy;
2296 #else
2297  Epetra_SerialComm SComm;
2298  Epetra_Map dummy(0,0,SComm);
2299 #endif
2300 
2301  delete CrsGraphData_->Importer_;
2302  delete CrsGraphData_->Exporter_;
2303 
2304 
2305  CrsGraphData_->RowMap_ = newMap ? *newMap : dummy;
2306  CrsGraphData_->ColMap_ = newColMap ? *newColMap : dummy;
2307  CrsGraphData_->DomainMap_ = newDomainMap ? *newDomainMap : dummy;
2308  CrsGraphData_->RangeMap_ = newRangeMap ? *newRangeMap : dummy;
2309  CrsGraphData_->Importer_ = newImport;
2310  CrsGraphData_->Exporter_ = newExport;
2311 
2312  // Epetra_DistObject things
2313  if(newMap) {
2314  Map_ = *newMap;
2315  Comm_ = &newMap->Comm();
2316  }
2317 
2318  // Cleanup (newRowMap is always newMap, so don't delete that)
2319  if(newColMap != newMap) delete newColMap;
2320  if(newDomainMap != newMap) delete newDomainMap;
2321  if(newRangeMap != newMap) delete newRangeMap;
2322 
2323  return(0);
2324 }
2325 
2326 // private =====================================================================
2328  try {
2329  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
2330  if(!A.GlobalConstantsComputed())
2331  EPETRA_CHK_ERR(-1); // Must have global constants to proceed
2332  }
2333  catch(...) {
2334  return(0); // No error at this point, object could be a RowMatrix
2335  }
2336  return(0);
2337 }
2338 
2339 // private =====================================================================
2341  int NumSameIDs,
2342  int NumPermuteIDs,
2343  int* PermuteToLIDs,
2344  int* PermuteFromLIDs,
2345  const Epetra_OffsetIndex * Indexor,
2346  Epetra_CombineMode CombineMode)
2347 {
2348  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2349  throw ReportError("Epetra_CrsGraph::CopyAndPermute: Incoming global index type does not match the one for *this",-1);
2350 
2351  try {
2352  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2353  EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
2354  PermuteFromLIDs,Indexor,CombineMode));
2355  }
2356  catch(...) {
2357  try {
2358  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2359  EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
2360  PermuteFromLIDs,Indexor,CombineMode));
2361  }
2362  catch(...) {
2363  EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
2364  }
2365  }
2366 
2367  return(0);
2368 }
2369 
2370 // private =====================================================================
2371 template<typename int_type>
2373  int NumSameIDs,
2374  int NumPermuteIDs,
2375  int* PermuteToLIDs,
2376  int* PermuteFromLIDs,
2377  const Epetra_OffsetIndex * Indexor,
2378  Epetra_CombineMode CombineMode)
2379 {
2380  (void)Indexor;
2381  (void)CombineMode;
2382  int i;
2383  int j;
2384  int NumIndices;
2385  int FromRow;
2386  int_type ToRow;
2387  int maxNumIndices = A.MaxNumEntries();
2388  Epetra_IntSerialDenseVector local_indices_vec;
2389 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2390  Epetra_LongLongSerialDenseVector global_indices_vec;
2391 #endif
2392  Epetra_SerialDenseVector Values;
2393 
2394  int* local_indices = 0;
2395  int_type* global_indices = 0;
2396 
2397  if(maxNumIndices > 0) {
2398  local_indices_vec.Size(maxNumIndices);
2399  local_indices = local_indices_vec.Values();
2400 
2401  if(A.RowMatrixRowMap().GlobalIndicesLongLong())
2402  {
2403 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2404  global_indices_vec.Size(maxNumIndices);
2405  global_indices = reinterpret_cast<int_type*>(global_indices_vec.Values());
2406 #else
2407  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: GlobalIndicesLongLong but no API for long long",-1);
2408 #endif
2409  }
2410  else
2411  {
2412  global_indices = reinterpret_cast<int_type*>(local_indices);
2413  }
2414 
2415  Values.Size(maxNumIndices); // Must extract values even though we discard them
2416  }
2417 
2418  const Epetra_Map& rowMap = A.RowMatrixRowMap();
2419  const Epetra_Map& colMap = A.RowMatrixColMap();
2420 
2421  // Do copy first
2422  for(i = 0; i < NumSameIDs; i++) {
2423  ToRow = (int) rowMap.GID64(i);
2424  EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumIndices, NumIndices, Values.Values(), local_indices));
2425  for(j = 0; j < NumIndices; j++)
2426  global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
2427  // Place into target graph.
2428  int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices);
2429  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2430  }
2431 
2432  // Do local permutation next
2433  for(i = 0; i < NumPermuteIDs; i++) {
2434  FromRow = PermuteFromLIDs[i];
2435  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2436  EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumIndices, NumIndices, Values.Values(), local_indices));
2437  for(j = 0; j < NumIndices; j++)
2438  global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
2439  int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices); // Place into target graph.
2440  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2441  }
2442 
2443  return(0);
2444 }
2445 
2447  int NumSameIDs,
2448  int NumPermuteIDs,
2449  int* PermuteToLIDs,
2450  int* PermuteFromLIDs,
2451  const Epetra_OffsetIndex * Indexor,
2452  Epetra_CombineMode CombineMode)
2453 {
2454  if(!A.RowMatrixRowMap().GlobalIndicesTypeMatch(RowMap()))
2455  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
2456 
2457 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2458  if(A.RowMatrixRowMap().GlobalIndicesInt())
2459  return CopyAndPermuteRowMatrix<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor,CombineMode);
2460  else
2461 #endif
2462 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2463  if(A.RowMatrixRowMap().GlobalIndicesLongLong())
2464  return CopyAndPermuteRowMatrix<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor,CombineMode);
2465  else
2466 #endif
2467  throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Unable to determine global index type of map", -1);
2468 }
2469 
2470 // private =====================================================================
2471 template<typename int_type>
2473  int NumSameIDs,
2474  int NumPermuteIDs,
2475  int* PermuteToLIDs,
2476  int* PermuteFromLIDs,
2477  const Epetra_OffsetIndex * Indexor,
2478  Epetra_CombineMode CombineMode)
2479 {
2480  (void)Indexor;
2481  (void)CombineMode;
2482  int i;
2483  int_type Row;
2484  int NumIndices;
2485  int_type* indices = 0;
2486  int_type FromRow, ToRow;
2487  int maxNumIndices = A.MaxNumIndices();
2488  Epetra_IntSerialDenseVector int_IndicesVector;
2489 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2490  Epetra_LongLongSerialDenseVector LL_IndicesVector;
2491 #endif
2492 
2493  if(maxNumIndices > 0 && A.IndicesAreLocal()) {
2494  if(A.RowMap().GlobalIndicesInt())
2495  {
2496  int_IndicesVector.Size(maxNumIndices);
2497  indices = reinterpret_cast<int_type*>(int_IndicesVector.Values());
2498  }
2499  else if(A.RowMap().GlobalIndicesLongLong())
2500  {
2501 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2502  LL_IndicesVector.Size(maxNumIndices);
2503  indices = reinterpret_cast<int_type*>(LL_IndicesVector.Values());
2504 #else
2505  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2506 #endif
2507  }
2508  }
2509 
2510  // Do copy first
2511  if(NumSameIDs > 0) {
2512  if(A.IndicesAreLocal()) {
2513  for(i = 0; i < NumSameIDs; i++) {
2514  Row = (int_type) GRID64(i);
2515  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumIndices, NumIndices, indices));
2516  // Place into target graph.
2517  int ierr = InsertGlobalIndices(Row, NumIndices, indices);
2518  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2519  }
2520  }
2521  else { // A.IndiceAreGlobal()
2522  for(i = 0; i < NumSameIDs; i++) {
2523  Row = (int_type) GRID64(i);
2524  EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, indices));
2525  // Place into target graph.
2526  int ierr = InsertGlobalIndices(Row, NumIndices, indices);
2527  if(ierr < 0) EPETRA_CHK_ERR(ierr);
2528  }
2529  }
2530  }
2531 
2532  // Do local permutation next
2533  if(NumPermuteIDs > 0) {
2534  if(A.IndicesAreLocal()) {
2535  for(i = 0; i < NumPermuteIDs; i++) {
2536  FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
2537  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2538  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2539  // Place into target graph.
2540  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2541  if (ierr < 0) EPETRA_CHK_ERR(ierr);
2542  }
2543  }
2544  else { // A.IndiceAreGlobal()
2545  for(i = 0; i < NumPermuteIDs; i++) {
2546  FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
2547  ToRow = (int_type) GRID64(PermuteToLIDs[i]);
2548  EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, indices));
2549  // Place into target graph.
2550  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2551  if (ierr < 0) EPETRA_CHK_ERR(ierr);
2552  }
2553  }
2554  }
2555 
2556  return(0);
2557 }
2558 
2560  int NumSameIDs,
2561  int NumPermuteIDs,
2562  int* PermuteToLIDs,
2563  int* PermuteFromLIDs,
2564  const Epetra_OffsetIndex * Indexor,
2565  Epetra_CombineMode CombineMode)
2566 {
2567  if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
2568  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Incoming global index type does not match the one for *this",-1);
2569 
2570  if(A.RowMap().GlobalIndicesInt())
2571 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2572  return CopyAndPermuteCrsGraph<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2573 #else
2574  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesInt but no API for it.",-1);
2575 #endif
2576 
2577  if(A.RowMap().GlobalIndicesLongLong())
2578 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2579  return CopyAndPermuteCrsGraph<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2580 #else
2581  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2582 #endif
2583 
2584  throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Unable to determine global index type of map", -1);
2585 }
2586 
2587 // private =====================================================================
2589  int NumExportIDs,
2590  int* ExportLIDs,
2591  int& LenExports,
2592  char*& Exports,
2593  int& SizeOfPacket,
2594  int * Sizes,
2595  bool& VarSizes,
2596  Epetra_Distributor& Distor)
2597 {
2598  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2599  throw ReportError("Epetra_CrsGraph::PackAndPrepare: Incoming global index type does not match the one for *this",-1);
2600 
2601  int globalMaxNumIndices = 0;
2602  int TotalSendSize = 0;
2603 
2604  VarSizes = true;
2605 
2606  if(Source.Map().GlobalIndicesInt())
2607  SizeOfPacket = (int)sizeof(int);
2608  else if(Source.Map().GlobalIndicesLongLong())
2609  SizeOfPacket = (int)sizeof(long long);
2610  else
2611  throw ReportError("Epetra_CrsGraph::PackAndPrepare: Unable to determine source global index type",-1);
2612 
2613  if(NumExportIDs <= 0) return(0);
2614 
2615  try {
2616  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2617  globalMaxNumIndices = A.GlobalMaxNumIndices();
2618  for( int i = 0; i < NumExportIDs; ++i )
2619  {
2620  Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
2621  TotalSendSize += Sizes[i];
2622  }
2623  }
2624  catch(...) {
2625  try {
2626  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2627  int maxNumIndices = A.MaxNumEntries();
2628  A.Comm().MaxAll(&maxNumIndices, &globalMaxNumIndices, 1);
2629  for( int i = 0; i < NumExportIDs; ++i )
2630  {
2631  int NumEntries;
2632  A.NumMyRowEntries( ExportLIDs[i], NumEntries );
2633  Sizes[i] = (NumEntries + 2);
2634  TotalSendSize += Sizes[i];
2635  }
2636  }
2637  catch(...) {
2638  EPETRA_CHK_ERR(-1); // Bad cast
2639  }
2640  }
2641 
2642  CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
2643 
2644  try {
2645  const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
2646  EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
2647  SizeOfPacket, Sizes, VarSizes, Distor));
2648  }
2649  catch(...) {
2650  const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
2651  EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
2652  SizeOfPacket, Sizes, VarSizes, Distor));
2653  }
2654  return(0);
2655 }
2656 
2657 // private =====================================================================
2659  int NumExportIDs,
2660  int* ExportLIDs,
2661  int& LenExports,
2662  char*& Exports,
2663  int& SizeOfPacket,
2664  int* Sizes,
2665  bool& VarSizes,
2666  Epetra_Distributor& Distor)
2667 {
2668  if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
2669  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Incoming global index type does not match the one for *this",-1);
2670 
2671  (void)LenExports;
2672  (void)SizeOfPacket;
2673  (void)Sizes;
2674  (void)VarSizes;
2675  (void)Distor;
2676  int i;
2677  int NumIndices;
2678 
2679  // Each segment of Exports will be filled by a packed row of information for each row as follows:
2680  // 1st int: GRID of row where GRID is the global row ID for the source graph
2681  // next int: NumIndices, Number of indices in row.
2682  // next NumIndices: The actual indices for the row.
2683  // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
2684  // sized segments for current communication routines.
2685  int maxNumIndices = A.MaxNumIndices();
2686  //if( maxNumIndices ) indices = new int[maxNumIndices];
2687 
2688  if(A.RowMap().GlobalIndicesInt()) {
2689  int* indices = 0;
2690  int* intptr = (int*) Exports;
2691  int FromRow;
2692  for(i = 0; i < NumExportIDs; i++) {
2693  FromRow = (int) A.GRID64(ExportLIDs[i]);
2694  *intptr = FromRow;
2695  indices = intptr + 2;
2696  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2697  intptr[1] = NumIndices; // Load second slot of segment
2698  intptr += (NumIndices+2); // Point to next segment
2699  }
2700  }
2701  else if(A.RowMap().GlobalIndicesLongLong()) {
2702 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2703  long long* indices = 0;
2704  long long* LLptr = (long long*) Exports;
2705  long long FromRow;
2706  for(i = 0; i < NumExportIDs; i++) {
2707  FromRow = A.GRID64(ExportLIDs[i]);
2708  *LLptr = FromRow;
2709  indices = LLptr + 2;
2710  EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
2711  LLptr[1] = NumIndices; // Load second slot of segment
2712  LLptr += (NumIndices+2); // Point to next segment
2713  }
2714 #else
2715  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2716 #endif
2717  }
2718  else
2719  throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Unable to determine source global index type",-1);
2720 
2721  //if( indices ) delete [] indices;
2722 
2723  return(0);
2724 }
2725 
2726 // private =====================================================================
2728  int NumExportIDs,
2729  int* ExportLIDs,
2730  int& LenExports,
2731  char*& Exports,
2732  int& SizeOfPacket,
2733  int* Sizes,
2734  bool& VarSizes,
2735  Epetra_Distributor& Distor)
2736 {
2737  if(!A.Map().GlobalIndicesTypeMatch(RowMap()))
2738  throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Incoming global index type does not match the one for *this",-1);
2739 
2740  (void)LenExports;
2741  (void)SizeOfPacket;
2742  (void)Sizes;
2743  (void)VarSizes;
2744  (void)Distor;
2745  int i;
2746  int j;
2747  int NumIndices;
2748  Epetra_SerialDenseVector Values;
2749 
2750  // Each segment of Exports will be filled by a packed row of information for each row as follows:
2751  // 1st int: GRID of row where GRID is the global row ID for the source graph
2752  // next int: NumIndices, Number of indices in row.
2753  // next NumIndices: The actual indices for the row.
2754  // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
2755  // sized segments for current communication routines.
2756  int maxNumIndices = A.MaxNumEntries();
2757  if(maxNumIndices > 0) {
2758  Values.Size(maxNumIndices);
2759 // indices = new int[maxNumIndices];
2760  }
2761  const Epetra_Map& rowMap = A.RowMatrixRowMap();
2762  const Epetra_Map& colMap = A.RowMatrixColMap();
2763 
2764  if(rowMap.GlobalIndicesInt() && colMap.GlobalIndicesInt()) {
2765  int* indices = 0;
2766  int FromRow;
2767  int* intptr = (int*) Exports;
2768  for(i = 0; i < NumExportIDs; i++) {
2769  FromRow = (int) rowMap.GID64(ExportLIDs[i]);
2770  *intptr = FromRow;
2771  indices = intptr + 2;
2772  EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), indices));
2773  for(j = 0; j < NumIndices; j++) indices[j] = (int) colMap.GID64(indices[j]); // convert to GIDs
2774  intptr[1] = NumIndices; // Load second slot of segment
2775  intptr += (NumIndices+2); // Point to next segment
2776  }
2777  }
2778  else if(rowMap.GlobalIndicesLongLong() && colMap.GlobalIndicesLongLong()) {
2779  // Bytes of Exports:
2780  // 12345678.12345678....12345678.12345678 ("." means no spaces)
2781  // FromRow NumIndices id1 id2 id3 id4 <-- before converting to GIDs
2782  // FromRow NumIndices | gid1 | | gid2 | <-- after converting to GIDs
2783 
2784  long long* LL_indices = 0;
2785  long long FromRow;
2786  long long* LLptr = (long long*) Exports;
2787  for(i = 0; i < NumExportIDs; i++) {
2788  FromRow = rowMap.GID64(ExportLIDs[i]);
2789  *LLptr = FromRow;
2790  LL_indices = LLptr + 2;
2791  int* int_indices = reinterpret_cast<int*>(LL_indices);
2792  EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), int_indices));
2793 
2794  // convert to GIDs, start from right.
2795  for(j = NumIndices; j > 0;) {
2796  --j;
2797  LL_indices[j] = colMap.GID64(int_indices[j]);
2798  }
2799 
2800  LLptr[1] = NumIndices; // Load second slot of segment
2801  LLptr += (NumIndices+2); // Point to next segment
2802  }
2803  }
2804  else
2805  throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Unable to determine source global index type",-1);
2806 
2807 
2808 // if( indices ) delete [] indices;
2809 
2810  return(0);
2811 }
2812 
2813 // private =====================================================================
2815  int NumImportIDs,
2816  int* ImportLIDs,
2817  int LenImports,
2818  char* Imports,
2819  int& SizeOfPacket,
2820  Epetra_Distributor& Distor,
2821  Epetra_CombineMode CombineMode,
2822  const Epetra_OffsetIndex * Indexor)
2823 {
2824  if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
2825  throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Incoming global index type does not match the one for *this",-1);
2826 
2827  (void)Source;
2828  (void)LenImports;
2829  (void)SizeOfPacket;
2830  (void)Distor;
2831  (void)CombineMode;
2832  (void)Indexor;
2833  if(NumImportIDs <= 0)
2834  return(0);
2835 
2836  // Unpack it...
2837 
2838  // Each segment of Sends will be filled by a packed row of information for each row as follows:
2839  // 1st int: GRID of row where GRID is the global row ID for the source graph
2840  // next int: NumIndices, Number of indices in row.
2841  // next NumIndices: The actual indices for the row.
2842 
2843 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2844  if(Source.Map().GlobalIndicesInt()) {
2845  int NumIndices;
2846  int i;
2847  int* indices;
2848  int ToRow;
2849  int* intptr = (int*) Imports;
2850  for(i = 0; i < NumImportIDs; i++) {
2851  ToRow = (int) GRID64(ImportLIDs[i]);
2852  assert((intptr[0])==ToRow); // Sanity check
2853  NumIndices = intptr[1];
2854  indices = intptr + 2;
2855  // Insert indices
2856  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2857  if(ierr < 0)
2858  EPETRA_CHK_ERR(ierr);
2859  intptr += (NumIndices+2); // Point to next segment
2860  }
2861  }
2862  else
2863 #endif
2864 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2865  if(Source.Map().GlobalIndicesLongLong()) {
2866  int NumIndices;
2867  int i;
2868  long long* indices;
2869  long long ToRow;
2870  long long* LLptr = (long long*) Imports;
2871  for(i = 0; i < NumImportIDs; i++) {
2872  ToRow = GRID64(ImportLIDs[i]);
2873  assert((LLptr[0])==ToRow); // Sanity check
2874  NumIndices = (int) LLptr[1];
2875  indices = LLptr + 2;
2876  // Insert indices
2877  int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
2878  if(ierr < 0)
2879  EPETRA_CHK_ERR(ierr);
2880  LLptr += (NumIndices+2); // Point to next segment
2881  }
2882  }
2883  else
2884 #endif
2885  throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Unable to determine source global index type",-1);
2886 
2887  //destroy buffers since this operation is usually only done once
2888  if( LenExports_ ) {
2889  delete [] Exports_;
2890  Exports_ = 0;
2891  LenExports_ = 0;
2892  }
2893  if( LenImports_ ) {
2894  delete [] Imports_;
2895  Imports_ = 0;
2896  LenImports_ = 0;
2897  }
2898 
2899  return(0);
2900 }
2901 
2902 // protected ===================================================================
2904  int mineComputed = 0;
2905  int allComputed;
2907  mineComputed = 1;
2908  RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
2909  // If allComputed comes back zero then at least one PE need global constants recomputed.
2910  return(allComputed==1);
2911 }
2912 
2913 // private =====================================================================
2915  int myIndicesAreLocal = 0;
2916  int myIndicesAreGlobal = 0;
2918  myIndicesAreLocal = 1;
2920  myIndicesAreGlobal = 1;
2921  int allIndicesAreLocal;
2922  int allIndicesAreGlobal;
2923  RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
2924  RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
2925  CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
2926  CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1); // If indices are global on one PE should be local on all
2927 }
2928 
2929 //==============================================================================
2930 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2931 int *Epetra_CrsGraph::All_IndicesPlus1() const {
2932  // This functionality is needed because MKL "sparse matrix" "dense matrix"
2933  // functions do not work with column-based dense storage and zero-based
2934  // sparse storage. So add "1" to indices and save duplicate data. This means
2935  // we will use one-based indices. This does not affect sparse-matrix and vector
2936  // operations.
2937 
2938  int* ptr = 0;
2939  if (!StorageOptimized()) {
2940  throw ReportError("Epetra_CrsGraph: int *All_IndicesPlus1() cannot be called when StorageOptimized()==false", -1);
2941  }
2942  else {
2944 
2945  if(!vec.Length()) {
2946  int* indices = All_Indices();
2948  ptr = vec.Values();
2949  for(int i = 0; i < CrsGraphData_->NumMyNonzeros_; ++i)
2950  ptr[i] = indices[i] + 1;
2951  }
2952  else {
2953  ptr = vec.Values();
2954  }
2955  }
2956  return ptr;
2957 }
2958 #endif // defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
2959 
2960 //==============================================================================
2961 void Epetra_CrsGraph::Print (std::ostream& os) const {
2962  int MyPID = RowMap().Comm().MyPID();
2963  int NumProc = RowMap().Comm().NumProc();
2964 
2965  for(int iproc = 0; iproc < NumProc; iproc++) {
2966  if(MyPID == iproc) {
2967  if(MyPID == 0) {
2968  os << "\nNumber of Global Block Rows = " << NumGlobalBlockRows64() << std::endl;
2969  os << "Number of Global Block Cols = " << NumGlobalBlockCols64() << std::endl;
2970  os << "Number of Global Block Diags = " << NumGlobalBlockDiagonals64() << std::endl;
2971  os << "Number of Global Entries = " << NumGlobalEntries64() << std::endl;
2972  os << "\nNumber of Global Rows = " << NumGlobalRows64() << std::endl;
2973  os << "Number of Global Cols = " << NumGlobalCols64() << std::endl;
2974  os << "Number of Global Diagonals = " << NumGlobalDiagonals64() << std::endl;
2975  os << "Number of Global Nonzeros = " << NumGlobalNonzeros64() << std::endl;
2976  os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim() << std::endl;
2977  os << "Global Maximum Block Col Dim = " << GlobalMaxColDim() << std::endl;
2978  os << "Global Maximum Num Indices = " << GlobalMaxNumIndices() << std::endl;
2979  if(LowerTriangular()) os << " ** Matrix is Lower Triangular **" << std::endl;
2980  if(UpperTriangular()) os << " ** Matrix is Upper Triangular **" << std::endl;
2981  if(NoDiagonal()) os << " ** Matrix has no diagonal **" << std::endl << std::endl;
2982  }
2983  os << "\nNumber of My Block Rows = " << NumMyBlockRows() << std::endl;
2984  os << "Number of My Block Cols = " << NumMyBlockCols() << std::endl;
2985  os << "Number of My Block Diags = " << NumMyBlockDiagonals() << std::endl;
2986  os << "Number of My Entries = " << NumMyEntries() << std::endl;
2987  os << "\nNumber of My Rows = " << NumMyRows() << std::endl;
2988  os << "Number of My Cols = " << NumMyCols() << std::endl;
2989  os << "Number of My Diagonals = " << NumMyDiagonals() << std::endl;
2990  os << "Number of My Nonzeros = " << NumMyNonzeros() << std::endl;
2991  os << "\nMy Maximum Block Row Dim = " << MaxRowDim() << std::endl;
2992  os << "My Maximum Block Col Dim = " << MaxColDim() << std::endl;
2993  os << "My Maximum Num Indices = " << MaxNumIndices() << std::endl << std::endl;
2994 
2995  int NumMyBlockRows1 = NumMyBlockRows();
2996  int MaxNumIndices1 = MaxNumIndices();
2997  Epetra_IntSerialDenseVector Indices1_int(MaxNumIndices1);
2998 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2999  Epetra_LongLongSerialDenseVector Indices1_LL(MaxNumIndices1);
3000 #endif
3001 
3002  if(RowMap().GlobalIndicesInt()) {
3003  Indices1_int.Resize(MaxNumIndices1);
3004  }
3005  else if(RowMap().GlobalIndicesLongLong()) {
3006 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
3007  Indices1_LL.Resize(MaxNumIndices1);
3008 #else
3009  throw ReportError("Epetra_CrsGraph::Print: GlobalIndicesLongLong but no long long API",-1);
3010 #endif
3011  }
3012  else
3013  throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
3014 
3015  int NumIndices1;
3016  int i;
3017  int j;
3018 
3019  os.width(14);
3020  os << " Row Index "; os << " ";
3021  for(j = 0; j < MaxNumIndices(); j++) {
3022  os.width(12);
3023  os << "Col Index"; os << " ";
3024  }
3025  os << std::endl;
3026  for(i = 0; i < NumMyBlockRows1; i++) {
3027  if(RowMap().GlobalIndicesInt()) {
3028  int Row = (int) GRID64(i); // Get global row number
3029  ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_int.Values());
3030  os.width(14);
3031  os << Row ; os << " ";
3032  for(j = 0; j < NumIndices1 ; j++) {
3033  os.width(12);
3034  os << Indices1_int[j]; os << " ";
3035  }
3036  os << std::endl;
3037  }
3038  else if(RowMap().GlobalIndicesLongLong()) {
3039 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
3040  long long Row = GRID64(i); // Get global row number
3041  ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_LL.Values());
3042  os.width(14);
3043  os << Row ; os << " ";
3044  for(j = 0; j < NumIndices1 ; j++) {
3045  os.width(12);
3046  os << Indices1_LL[j]; os << " ";
3047  }
3048  os << std::endl;
3049 #else
3050  throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
3051 #endif
3052  }
3053  }
3054  os << std::flush;
3055  }
3056  // Do a few global ops to give I/O a chance to complete
3057  RowMap().Comm().Barrier();
3058  RowMap().Comm().Barrier();
3059  RowMap().Comm().Barrier();
3060  }
3061 }
3062 
3063 //==============================================================================
3065  if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
3066  return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
3067 
3068  CleanupData();
3069  CrsGraphData_ = Source.CrsGraphData_;
3071 
3072  return(*this);
3073 }
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
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.