Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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 =====================================================================
122template<typename int_type>
123int 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
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
202int 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/*
224int 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//==============================================================================
240
241// private =====================================================================
243 if(CrsGraphData_ != 0) {
245 if(CrsGraphData_->ReferenceCount() == 0) {
246 delete CrsGraphData_;
247 CrsGraphData_ = 0;
248 }
249 }
250}
251
252//==============================================================================
253template<typename int_type>
254int 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
260 int locRow = LRID(Row); // Find local row number for this global row index
261
262 EPETRA_CHK_ERR(InsertIndicesIntoSorted(locRow, NumIndices, indices));
263
265 return(1);
266 else
267 return(0);
268}
269
270//==============================================================================
271#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
272int 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
282int 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//==============================================================================
291int 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
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
315 return(1);
316 else
317 return(0);
318}
320// protected ===================================================================
321template<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
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// =========================================================================
440template<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
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//==============================================================================
538template<typename int_type>
539int 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--;
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
581 return(1);
582 else
583 return(0);
584}
585
586//==============================================================================
587#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
588int 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
598int 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//==============================================================================
607int 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
648 return(1);
649 else
650 return(0);
651}
652
653//==============================================================================
654template<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
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
753 return(1);
754 else
755 return(0);
756}
757
758// protected ===================================================================
759template<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 ===================================================================
839template<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//==============================================================================
969int 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
988 return(1);
989 else
990 return(0);
991}
992
993//==============================================================================
997
998//==============================================================================
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();
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
1106void 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();
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
1181void 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;
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;
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;
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
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
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
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//==============================================================================
2017template<typename int_type>
2018int 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
2049int 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
2059int 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//==============================================================================
2069template<typename int_type>
2070int 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
2091int 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
2101int 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
2123int 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//==============================================================================
2144int 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//==============================================================================
2160int 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{
2202 CrsGraphData_->ColMap_ = newmap;
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//==============================================================================
2225 int rv=0;
2226 if( !NewImporter && ColMap().SameAs(NewDomainMap)) {
2227 CrsGraphData_->DomainMap_ = NewDomainMap;
2228 delete CrsGraphData_->Importer_;
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
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 =====================================================================
2365template<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
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 =====================================================================
2465template<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
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{
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;
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)
2925int *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//==============================================================================
2955void 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();
3065
3066 return(*this);
3067}
3068
3069//=============================================================================
3073
3074//=============================================================================
Epetra_CombineMode
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
void epetra_crsgraph_compress_out_duplicates(int len, int *list, int &newlen)
void epetra_shellsort(int *list, int length)
Epetra_DataAccess
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.
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,...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
long long IndexBase64() const
int RemoteIDList(int NumIDs, const int *GIDList, int *PIDList, int *LIDList) const
Returns the processor IDs and corresponding local index value for a given list of global indices.
int MaxElementSize() const
Maximum element size across all processors.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
bool GlobalIndicesTypeMatch(const Epetra_BlockMap &other) const
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor.
long long GID64(int LID) const
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
long long NumGlobalElements64() const
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
Epetra_BlockMap * ReplaceCommWithSubset(const Epetra_Comm *Comm) const
Replace this BlockMap's communicator with a subset communicator.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool ConstantElementSize() const
Returns true if map has constant element size.
int NumMyElements() const
Number of elements on the calling processor.
bool MyLID(int lid) const
Returns true if the LID passed in belongs to the calling processor in this map, otherwise returns fal...
long long NumGlobalPoints64() const
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition Epetra_Comm.h:73
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
virtual int MyPID() const =0
Return my process ID.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
Epetra_CrsGraphData: The Epetra CrsGraph Data Class.
int ReAllocateAndCast(char *&UserPtr, int &Length, const int IntPacketSizeTimesNumTrans)
called by PackAndPrepare
IndexData< int > * data
IndexData< int_type > & Data()
const Epetra_Import * Importer_
Epetra_IntSerialDenseVector NumIndicesPerRow_
Epetra_IntSerialDenseVector NumAllocatedIndicesPerRow_
int MakeImportExport()
called by FillComplete (and TransformToLocal)
const Epetra_Export * Exporter_
Epetra_IntSerialDenseVector IndexOffset_
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
int MaxNumIndices() const
Returns the maximum number of nonzero entries across all rows on this processor.
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().
int NumMyBlockCols() const
Returns the number of Block matrix columns on this processor.
const Epetra_BlockMap & DomainMap() const
Returns the DomainMap associated with this graph.
int RemoveMyIndices(int LocalRow, int NumIndices, int *Indices)
Remove a list of elements from a specified local row of the graph.
bool StaticProfile() const
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
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.
const Epetra_BlockMap & RowMap() const
Returns the RowMap associated with this graph.
int NumMyDiagonals() const
Returns the number of diagonal entries in the local graph, based on global row/column index compariso...
bool IndicesAreContiguous() const
int NumAllocatedMyIndices(int Row) const
Returns the allocated number of nonzero entries in specified local row on this processor.
bool FindGlobalIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
int MakeColMap(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
int GlobalMaxRowDim() const
Returns the max row dimension of block entries across all processors.
int TRemoveGlobalIndices(long long Row)
Epetra_IntSerialDenseVector & ExpertExtractIndices()
Returns a reference to the Epetra_IntSerialDenseVector used to hold the local All_Indices (CRS colind...
long long NumGlobalEntries64() const
const Epetra_CrsGraphData * DataPtr() const
Returns a pointer to the CrsGraphData instance this CrsGraph uses.
int NumAllocatedGlobalIndices(long long Row) const
Returns the allocated number of nonzero entries in specified global row on this processor.
int MakeIndicesLocal(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
void SetIndicesAreLocal(bool Flag)
int * All_Indices() const
Epetra_IntSerialDenseVector & ExpertExtractIndexOffset()
Returns a reference to the Epetra_IntSerialDenseVector used to hold the local IndexOffsets (CRS rowpt...
void SetIndicesAreContiguous(bool Flag)
int MaxColDim() const
Returns the max column dimension of block entries on the processor.
int SortIndices()
Sort column indices, row-by-row, in ascending order.
int CopyAndPermuteRowMatrix(const Epetra_RowMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
int ** Indices() const
int GlobalMaxNumIndices() const
Returns the maximun number of nonzero entries across all rows across all processors.
bool LowerTriangular() const
If graph is lower triangular in local index space, this query returns true, otherwise it returns fals...
int TAllocate(const int *numIndicesPerRow, int Inc, bool staticProfile)
long long GCID64(int LCID_in) const
Epetra_CrsGraphData * CrsGraphData_
bool HaveColMap() const
Returns true if we have a well-defined ColMap, and returns false otherwise.
long long NumGlobalDiagonals64() const
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 GlobalConstantsComputed() const
bool Sorted() const
If SortIndices() has been called, this query returns true, otherwise it returns false.
int ExtractGlobalRowView(int GlobalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified global row of the graph.
void SetNoRedundancies(bool Flag)
virtual ~Epetra_CrsGraph()
Epetra_CrsGraph Destructor.
bool IndicesAreGlobal() const
If column indices are in global range, this query returns true, otherwise it returns false.
void SetSorted(bool Flag)
long long NumGlobalBlockRows64() const
int RemoveGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Remove a list of elements from a specified global row of the graph.
void SetAllocated(bool Flag)
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 RemoveRedundantIndices()
Removes any redundant column indices in the rows of the graph.
int TransformToLocal()
Use FillComplete() instead.
bool UpperTriangular() const
If graph is upper triangular in local index space, this query returns true, otherwise it returns fals...
long long NumGlobalBlockDiagonals64() const
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this graph.
int InsertIndices(int Row, int NumIndices, int *Indices)
int MakeColMap_int(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not.
int GlobalMaxColDim() const
Returns the max column dimension of block entries across all processors.
bool NoRedundancies() const
If RemoveRedundantIndices() has been called, this query returns true, otherwise it returns false.
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int PackAndPrepareCrsGraph(const Epetra_CrsGraph &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
void SetFilled(bool Flag)
int NumMyEntries() const
Returns the number of entries on this processor.
bool FindMyIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
int NumMyBlockRows() const
Returns the number of block matrix rows on this processor.
Epetra_CrsGraph & operator=(const Epetra_CrsGraph &Source)
Assignment operator.
long long NumGlobalRows64() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified local row of the graph.
virtual void Print(std::ostream &os) const
Print method.
int NumMyRows() const
Returns the number of matrix rows on this processor.
void SetIndicesAreGlobal(bool Flag)
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 ReplaceDomainMapAndImporter(const Epetra_BlockMap &NewDomainMap, const Epetra_Import *NewImporter)
Replaces the current DomainMap & Importer with the user-specified map object.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
int MakeColMap_LL(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
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...
long long NumGlobalNonzeros64() const
int CopyAndPermuteCrsGraph(const Epetra_CrsGraph &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
long long GRID64(int LRID_in) const
int NumMyBlockDiagonals() const
Returns the number of Block diagonal entries in the local graph, based on global row/column index com...
long long NumGlobalBlockCols64() const
int NumGlobalIndices(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
int Allocate(const int *NumIndicesPerRow, int Inc, bool StaticProfile)
int MaxRowDim() const
Returns the max row dimension of block entries on the processor.
void SetGlobalConstantsComputed(bool Flag)
const Epetra_BlockMap & RangeMap() const
Returns the RangeMap associated with this graph.
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().
const Epetra_BlockMap & ColMap() const
Returns the Column Map associated with this graph.
int TInsertIndicesIntoSorted(int Row, int NumIndices, int_type *Indices)
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
int PackAndPrepareRowMatrix(const Epetra_RowMatrix &A, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
int TInsertIndices(int Row, int NumIndices, int_type *Indices)
long long IndexBase64() const
int InsertIndicesIntoSorted(int Row, int NumIndices, int *Indices)
long long NumGlobalCols64() const
int NumMyNonzeros() const
Returns the number of indices in the local graph.
const Epetra_Import * Importer() const
Returns the Importer associated with this graph.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
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...
bool NoDiagonal() const
If graph has no diagonal entries in global index space, this query returns true, otherwise it returns...
int ReplaceRowMap(const Epetra_BlockMap &newmap)
Replaces the current RowMap with the user-specified map object, but only if currentmap->PointSameAs(n...
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
void DecrementReferenceCount()
Decrement reference count.
int ReferenceCount() const
Get reference count.
void IncrementReferenceCount()
Increment reference count.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
const Epetra_Comm * Comm_
Epetra_BlockMap Map_
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements.
void Add(const long long key, const value_type value)
value_type Get(const long long key)
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this importer.
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
int Resize(int Length_in)
Resize a Epetra_IntSerialDenseVector object.
int Size(int Length_in)
Set length of a Epetra_IntSerialDenseVector object; init values to zero.
int Length() const
Returns length of vector.
int * Values()
Returns pointer to the values in vector.
int MakeViewOf(const Epetra_IntSerialDenseVector &Source)
Reset an existing IntSerialDenseVector to point to another Vector.
Epetra_LongLongSerialDenseVector: A class for constructing and using dense vectors.
int Size(int Length_in)
Set length of a Epetra_LongLongSerialDenseVector object; init values to zero.
long long * Values()
Returns pointer to the values in vector.
int Resize(int Length_in)
Resize a Epetra_LongLongSerialDenseVector object.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual const Epetra_Map & RowMatrixColMap() const =0
Returns the Epetra_Map object associated with the columns of this matrix.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
Returns the number of nonzero entries in MyRow.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
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.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
double * Values() const
Returns pointer to the values in vector.
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
virtual const Epetra_BlockMap & Map() const =0
Returns a reference to the Epetra_BlockMap for this object.
Epetra_Util: The Epetra Util Wrapper Class.
Definition Epetra_Util.h:79
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)
Epetra_IntSerialDenseVector All_IndicesPlus1_
Epetra_IntSerialDenseVector All_Indices_