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