Epetra Package Browser (Single Doxygen Collection) Development
Epetra_CrsGraph.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2001 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_CrsGraph.h"
00045 #include "Epetra_Import.h"
00046 #include "Epetra_Export.h"
00047 #include "Epetra_Distributor.h"
00048 #include "Epetra_Util.h"
00049 #include "Epetra_Comm.h"
00050 #include "Epetra_HashTable.h"
00051 #include "Epetra_BlockMap.h"
00052 #include "Epetra_Map.h"
00053 #include "Epetra_RowMatrix.h"
00054 #include "Epetra_IntSerialDenseVector.h"
00055 #include "Epetra_SerialDenseVector.h"
00056 #include "Epetra_OffsetIndex.h"
00057 
00058 //==============================================================================
00059 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV, 
00060         const Epetra_BlockMap& RowMap, 
00061         const int* NumIndicesPerRow, bool StaticProfile) 
00062   : Epetra_DistObject(RowMap, "Epetra::CrsGraph"),
00063     CrsGraphData_(new Epetra_CrsGraphData(CV, RowMap, StaticProfile))
00064 {
00065   Allocate(NumIndicesPerRow, 1, StaticProfile);
00066 }
00067 
00068 //==============================================================================
00069 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV, 
00070         const Epetra_BlockMap& RowMap, 
00071         int NumIndicesPerRow, bool StaticProfile) 
00072   : Epetra_DistObject(RowMap, "Epetra::CrsGraph"),
00073     CrsGraphData_(new Epetra_CrsGraphData(CV, RowMap, StaticProfile))
00074 {
00075   Allocate(&NumIndicesPerRow, 0, StaticProfile);
00076 }
00077 
00078 //==============================================================================
00079 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV, 
00080          const Epetra_BlockMap& RowMap, 
00081          const Epetra_BlockMap& ColMap, 
00082          const int* NumIndicesPerRow, bool StaticProfile) 
00083   : Epetra_DistObject(RowMap, "Epetra::CrsGraph"),
00084     CrsGraphData_(new Epetra_CrsGraphData(CV, RowMap, ColMap, StaticProfile))
00085 {
00086   Allocate(NumIndicesPerRow, 1, StaticProfile);
00087 }
00088 
00089 //==============================================================================
00090 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV, 
00091          const Epetra_BlockMap& RowMap, 
00092          const Epetra_BlockMap& ColMap, 
00093          int NumIndicesPerRow, bool StaticProfile) 
00094   : Epetra_DistObject(RowMap, "Epetra::CrsGraph"),
00095     CrsGraphData_(new Epetra_CrsGraphData(CV, RowMap, ColMap, StaticProfile))
00096 {
00097   Allocate(&NumIndicesPerRow, 0, StaticProfile);
00098 }
00099 
00100 //==============================================================================
00101 Epetra_CrsGraph::Epetra_CrsGraph(const Epetra_CrsGraph& Graph)
00102   : Epetra_DistObject(Graph),
00103     CrsGraphData_(Graph.CrsGraphData_)
00104 {
00105   CrsGraphData_->IncrementReferenceCount();
00106 }
00107 
00108 // private =====================================================================
00109 int Epetra_CrsGraph::Allocate(const int* NumIndicesPerRow, int Inc, bool StaticProfile) {
00110   int i;
00111   const int numMyBlockRows = CrsGraphData_->NumMyBlockRows_;
00112   
00113   // Portions specific to ISDVs
00114   // Note that NumIndicesPerRow_ will be 1 element longer than needed.
00115   // This is because, if OptimizeStorage() is called, the storage associated with
00116   // NumIndicesPerRow_ will be reused implicitly for the IndexOffset_ array.
00117   // Although a bit fragile, for users who care about efficient memory allocation,
00118   // the time and memory fragmentation are important: Mike Heroux Feb 2005.
00119   int errorcode = CrsGraphData_->NumIndicesPerRow_.Size(numMyBlockRows+1); 
00120   if(errorcode != 0) 
00121     throw ReportError("Error with NumIndicesPerRow_ allocation.", -99);
00122 
00123   errorcode = CrsGraphData_->NumAllocatedIndicesPerRow_.Size(numMyBlockRows);
00124   if(errorcode != 0) 
00125     throw ReportError("Error with NumAllocatedIndicesPerRow_ allocation.", -99);
00126 
00127   int nnz = 0;
00128   if(CrsGraphData_->CV_ == Copy) {
00129     if (NumIndicesPerRow != 0) {
00130       for(i = 0; i < numMyBlockRows; i++) {
00131   int nnzr = NumIndicesPerRow[i*Inc];
00132   nnz += nnzr;
00133       }
00134     }
00135   }
00136   CrsGraphData_->NumMyEntries_ = nnz; // Define this now since we have it and might want to use it
00137   //***
00138 
00139   CrsGraphData_->MaxNumIndices_ = 0;
00140   
00141   // Allocate and initialize entries if we are copying data
00142   if(CrsGraphData_->CV_ == Copy) {
00143     if (StaticProfile) CrsGraphData_->All_Indices_.Size(nnz);
00144     int * All_Indices = CrsGraphData_->All_Indices_.Values(); // First address of contiguous buffer
00145     for(i = 0; i < numMyBlockRows; i++) {
00146       const int NumIndices = NumIndicesPerRow==0 ? 0 :NumIndicesPerRow[i*Inc];
00147       const int indexBaseMinusOne = IndexBase() - 1;
00148 
00149       if(NumIndices > 0) {
00150   if (StaticProfile) {
00151     CrsGraphData_->Indices_[i] = All_Indices;
00152     All_Indices += NumIndices;
00153     int* ColIndices = CrsGraphData_->Indices_[i];
00154     for(int j = 0; j < NumIndices; j++) 
00155       ColIndices[j] = indexBaseMinusOne; // Fill column indices with out-of-range values
00156   }
00157   else {
00158     // reserve memory in the STL vector, and then resize it to zero
00159     // again in order to signal the program that no data is in there
00160     // yet.
00161     CrsGraphData_->SortedEntries_[i].entries_.resize(NumIndices,
00162                  indexBaseMinusOne);
00163     CrsGraphData_->Indices_[i] = NumIndices > 0 ? &CrsGraphData_->SortedEntries_[i].entries_[0]: NULL;
00164     CrsGraphData_->SortedEntries_[i].entries_.resize(0);
00165   }
00166       }
00167       else {
00168   CrsGraphData_->Indices_[i] = 0;
00169       }
00170 
00171       CrsGraphData_->NumAllocatedIndicesPerRow_[i] = NumIndices;
00172     }
00173     if (StaticProfile) assert(CrsGraphData_->All_Indices_.Values()+nnz==All_Indices); // Sanity check
00174   }  
00175   else { // CV_ == View
00176     for(i = 0; i < numMyBlockRows; i++) {
00177       CrsGraphData_->Indices_[i] = 0;
00178     }
00179   }
00180 
00181   SetAllocated(true);
00182 
00183   return(0);
00184 }
00185 
00186 // private =====================================================================
00187 /*
00188 int Epetra_CrsGraph::ReAllocate() {
00189   // Reallocate storage that was deleted in OptimizeStorage
00190 
00191   // Since NIPR is in Copy mode, NAIPR will become a Copy as well
00192   CrsGraphData_->NumAllocatedIndicesPerRow_ = CrsGraphData_->NumIndicesPerRow_;
00193 
00194   CrsGraphData_->StorageOptimized_ = false;
00195 
00196   return(0);
00197 }
00198 */
00199 //==============================================================================
00200 Epetra_CrsGraph::~Epetra_CrsGraph()
00201 {
00202   CleanupData();
00203 }
00204 
00205 // private =====================================================================
00206 void Epetra_CrsGraph::CleanupData() {
00207   if(CrsGraphData_ != 0) {
00208     CrsGraphData_->DecrementReferenceCount();
00209     if(CrsGraphData_->ReferenceCount() == 0) {
00210       delete CrsGraphData_;
00211       CrsGraphData_ = 0;
00212     }
00213   }
00214 }
00215 
00216 //==============================================================================
00217 int Epetra_CrsGraph::InsertGlobalIndices(int Row, int NumIndices, int* Indices) {
00218   if(IndicesAreLocal()) 
00219     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00220   if(IndicesAreContiguous()) 
00221     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00222   SetIndicesAreGlobal(true);
00223   Row = LRID(Row); // Find local row number for this global row index
00224 
00225   EPETRA_CHK_ERR(InsertIndicesIntoSorted(Row, NumIndices, Indices));
00226 
00227   if(CrsGraphData_->ReferenceCount() > 1)
00228     return(1);
00229   else
00230     return(0);
00231 }
00232 
00233 //==============================================================================
00234 int Epetra_CrsGraph::InsertMyIndices(int Row, int NumIndices, int* Indices) {
00235 
00236   if(IndicesAreGlobal()) {
00237     EPETRA_CHK_ERR(-2); // Cannot insert local indices into a global graph
00238   }
00239   if(IndicesAreContiguous()) 
00240     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00241 
00242   if (CrsGraphData_->HaveColMap_) {
00243     SetIndicesAreLocal(true);
00244   }
00245   else {
00246      if (!IndicesAreLocal()) {
00247        EPETRA_CHK_ERR(-4);
00248      }
00249   }
00250 
00251   EPETRA_CHK_ERR(InsertIndicesIntoSorted(Row, NumIndices, Indices));
00252 
00253   if(CrsGraphData_->ReferenceCount() > 1)
00254     return(1);
00255   else
00256     return(0);
00257 }
00258 
00259 // protected ===================================================================
00260 int Epetra_CrsGraph::InsertIndices(int Row,
00261            int NumIndices,
00262            int* UserIndices)
00263 {
00264   if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
00265 
00266   SetSorted(false); // No longer in sorted state.
00267   CrsGraphData_->NoRedundancies_ = false; // Redundancies possible.
00268   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00269 
00270   int j;
00271   int ierr = 0;
00272 
00273   if(Row < 0 || Row >= NumMyBlockRows()) 
00274     EPETRA_CHK_ERR(-2); // Not in Row range
00275     
00276   int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
00277   int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00278 
00279   if(CrsGraphData_->CV_ == View) {
00280     if(CrsGraphData_->Indices_[Row] != 0) 
00281       ierr = 2; // This row has been defined already.  Issue warning.
00282     CrsGraphData_->Indices_[Row] = UserIndices;
00283     current_numAllocIndices = NumIndices;
00284     current_numIndices = NumIndices;
00285   }
00286   else {
00287     // if HaveColMap_ is true, UserIndices is copied into a new array,
00288     // and then modified. The UserIndices pointer is updated to point to this 
00289     // new array. If HaveColMap_ is false, nothing is done. This way,
00290     // the same UserIndices pointer can be used later on regardless of whether
00291     // changes were made.
00292     if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
00293       if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
00294         delete [] CrsGraphData_->TempColIndices_;
00295         CrsGraphData_->TempColIndices_ = new int[NumIndices];
00296         CrsGraphData_->NumTempColIndices_ = NumIndices;
00297       }
00298       int * tempIndices = CrsGraphData_->TempColIndices_;
00299       int loc = 0;
00300       if(IndicesAreLocal()) {
00301         for(j = 0; j < NumIndices; ++j)
00302           if(CrsGraphData_->ColMap_.MyLID(UserIndices[j])) 
00303             tempIndices[loc++] = UserIndices[j];
00304       }
00305       else {
00306         for(j = 0; j < NumIndices; ++j) {
00307           const int Index = CrsGraphData_->ColMap_.LID(UserIndices[j]);
00308           if (Index > -1)
00309             tempIndices[loc++] = UserIndices[j];
00310         }
00311 
00312       }
00313       if(loc != NumIndices) 
00314         ierr = 2; //Some columns excluded
00315       NumIndices = loc;
00316       UserIndices = tempIndices;
00317     }
00318 
00319     int start = current_numIndices;
00320     int stop = start + NumIndices;
00321     if (CrsGraphData_->StaticProfile_) {
00322       if(stop > current_numAllocIndices)
00323         EPETRA_CHK_ERR(-2); // Cannot reallocate storage if graph created using StaticProfile
00324     }
00325     else {
00326       if (current_numAllocIndices > 0 && stop > current_numAllocIndices)
00327         ierr = 3;
00328       CrsGraphData_->SortedEntries_[Row].entries_.resize(stop, IndexBase() - 1);
00329       CrsGraphData_->Indices_[Row] = stop>0 ? &CrsGraphData_->SortedEntries_[Row].entries_[0] : NULL;
00330 
00331       current_numAllocIndices =  CrsGraphData_->SortedEntries_[Row].entries_.capacity();    
00332     }
00333 
00334     current_numIndices = stop;
00335     int* RowIndices = CrsGraphData_->Indices_[Row]+start;
00336     for(j = 0; j < NumIndices; j++) {
00337       RowIndices[j] = UserIndices[j];
00338     }
00339   }
00340 
00341   if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
00342     CrsGraphData_->MaxNumIndices_ = current_numIndices;
00343   }
00344   EPETRA_CHK_ERR(ierr);
00345 
00346 
00347   if(CrsGraphData_->ReferenceCount() > 1)
00348     return(1); // return 1 if data is shared
00349   else
00350     return(0);
00351 }
00352 
00353 // =========================================================================
00354 int Epetra_CrsGraph::InsertIndicesIntoSorted(int Row,
00355               int NumIndices,
00356               int* UserIndices)
00357 {
00358   // This function is only valid for COPY mode with non-static profile and
00359   // sorted entries. Otherwise, go to the other function.
00360   if (!CrsGraphData_->NoRedundancies_ || !CrsGraphData_->Sorted_  || 
00361       CrsGraphData_->StaticProfile_ || CrsGraphData_->CV_ == View ) 
00362     return InsertIndices(Row, NumIndices, UserIndices);
00363 
00364   if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
00365 
00366   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00367 
00368   int ierr = 0;
00369 
00370   if(Row < 0 || Row >= NumMyBlockRows()) 
00371     EPETRA_CHK_ERR(-2); // Not in Row range
00372     
00373   int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
00374   int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00375 
00376   // if HaveColMap_ is true, UserIndices filters out excluded indices,
00377   // and then modified. The UserIndices pointer is updated to point to this 
00378   // new array. If HaveColMap_ is false, nothing is done. This way,
00379   // the same UserIndices pointer can be used later on regardless of whether
00380   // changes were made.
00381   if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
00382     if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
00383       delete [] CrsGraphData_->TempColIndices_;
00384       CrsGraphData_->TempColIndices_ = new int[NumIndices];
00385       CrsGraphData_->NumTempColIndices_ = NumIndices;
00386     }
00387     int * tempIndices = CrsGraphData_->TempColIndices_;
00388     int loc = 0;
00389     if(IndicesAreLocal()) {
00390       for(int j = 0; j < NumIndices; ++j)
00391         if(CrsGraphData_->ColMap_.MyLID(UserIndices[j]))
00392           tempIndices[loc++] = UserIndices[j];
00393     }
00394     else {
00395       for(int j = 0; j < NumIndices; ++j) {
00396         if (CrsGraphData_->ColMap_.MyGID(UserIndices[j])) {
00397           tempIndices[loc++] = UserIndices[j];
00398         }
00399       }
00400     }
00401     if(loc != NumIndices) 
00402       ierr = 2; //Some columns excluded
00403     NumIndices = loc;
00404     UserIndices = tempIndices;
00405   }
00406 
00407   // for non-static profile, directly insert into a list that we always
00408   // keep sorted.
00409   CrsGraphData_->SortedEntries_[Row].AddEntries(NumIndices, UserIndices);
00410   current_numIndices = CrsGraphData_->SortedEntries_[Row].entries_.size();
00411   current_numAllocIndices = CrsGraphData_->SortedEntries_[Row].entries_.capacity();
00412   // reset the pointer to the respective data
00413   CrsGraphData_->Indices_[Row] = current_numIndices > 0 ? &CrsGraphData_->SortedEntries_[Row].entries_[0] : NULL;
00414 
00415   if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
00416     CrsGraphData_->MaxNumIndices_ = current_numIndices;
00417   }
00418   EPETRA_CHK_ERR(ierr);
00419 
00420   if(CrsGraphData_->ReferenceCount() > 1)
00421     return(1); // return 1 if data is shared
00422   else
00423     return(0);
00424 }
00425 
00426 //==============================================================================
00427 int Epetra_CrsGraph::RemoveGlobalIndices(int Row, int NumIndices, int* Indices) {
00428   int j;
00429   int k;
00430   int ierr = 0;
00431   int Loc;
00432 
00433   if(IndicesAreContiguous() || StorageOptimized()) 
00434     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00435 
00436   if(IndicesAreLocal()) 
00437     EPETRA_CHK_ERR(-2); // Cannot remove global indices from a filled graph
00438 
00439   if(CrsGraphData_->CV_ == View) 
00440     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00441 
00442   Row = LRID(Row); // Normalize row range
00443     
00444   if(Row < 0 || Row >= NumMyBlockRows()) 
00445     EPETRA_CHK_ERR(-1); // Not in Row range
00446     
00447   int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00448   
00449   for(j = 0; j < NumIndices; j++) {
00450     int Index = Indices[j];
00451     if(FindGlobalIndexLoc(Row,Index,j,Loc)) {
00452       for(k = Loc+1; k < NumCurrentIndices; k++) 
00453         CrsGraphData_->Indices_[Row][k-1] = CrsGraphData_->Indices_[Row][k];
00454       NumCurrentIndices--;
00455       CrsGraphData_->NumIndicesPerRow_[Row]--;
00456       if (!CrsGraphData_->StaticProfile_)
00457         CrsGraphData_->SortedEntries_[Row].entries_.pop_back();
00458       else
00459         CrsGraphData_->Indices_[Row][NumCurrentIndices-1] = IndexBase() - 1;
00460     }
00461   }
00462   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00463 
00464   EPETRA_CHK_ERR(ierr);
00465 
00466   if(CrsGraphData_->ReferenceCount() > 1)
00467     return(1);
00468   else
00469     return(0);
00470 }
00471 
00472 //==============================================================================
00473 int Epetra_CrsGraph::RemoveMyIndices(int Row, int NumIndices, int* Indices) {
00474 
00475   if(IndicesAreContiguous() || StorageOptimized()) 
00476     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00477 
00478   if(IndicesAreGlobal()) 
00479     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00480 
00481   int j;
00482   int k;
00483   int ierr = 0;
00484   int Loc;
00485 
00486   if(CrsGraphData_->CV_ == View) 
00487     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00488 
00489   if(Row < 0 || Row >= NumMyBlockRows()) 
00490     EPETRA_CHK_ERR(-1); // Not in Row range
00491     
00492   int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00493 
00494   for(j = 0; j < NumIndices; j++) {
00495     int Index = Indices[j];
00496     if(FindMyIndexLoc(Row,Index,j,Loc)) {
00497       for(k = Loc + 1; k < NumCurrentIndices; k++) 
00498         CrsGraphData_->Indices_[Row][k-1] = CrsGraphData_->Indices_[Row][k];
00499       NumCurrentIndices--;
00500       CrsGraphData_->NumIndicesPerRow_[Row]--;
00501       if (!CrsGraphData_->StaticProfile_)
00502         CrsGraphData_->SortedEntries_[Row].entries_.pop_back();
00503       else
00504         CrsGraphData_->Indices_[Row][NumCurrentIndices-1] = IndexBase() - 1;
00505     }
00506   }
00507   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00508 
00509   EPETRA_CHK_ERR(ierr);
00510 
00511   if(CrsGraphData_->ReferenceCount() > 1)
00512     return(1);
00513   else
00514     return(0);
00515 }
00516 
00517 //==============================================================================
00518 int Epetra_CrsGraph::RemoveGlobalIndices(int Row) {
00519   int j;
00520   int ierr = 0;
00521 
00522   if(IndicesAreContiguous() || StorageOptimized()) 
00523     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00524 
00525   if(IndicesAreLocal()) 
00526     EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
00527   if(CrsGraphData_->CV_ == View) 
00528     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00529 
00530   Row = LRID(Row); // Normalize row range
00531     
00532   if(Row < 0 || Row >= NumMyBlockRows()) 
00533     EPETRA_CHK_ERR(-1); // Not in Row range
00534 
00535   if (CrsGraphData_->StaticProfile_) {
00536     int NumIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00537   
00538     const int indexBaseMinusOne = IndexBase() - 1;
00539     for(j = 0; j < NumIndices; j++) 
00540       CrsGraphData_->Indices_[Row][j] = indexBaseMinusOne; // Set to invalid 
00541   }
00542   else
00543     CrsGraphData_->SortedEntries_[Row].entries_.resize(0);
00544  
00545   CrsGraphData_->NumIndicesPerRow_[Row] = 0;
00546 
00547 
00548   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00549   EPETRA_CHK_ERR(ierr);
00550 
00551   if(CrsGraphData_->ReferenceCount() > 1)
00552     return(1);
00553   else
00554     return(0);
00555 }
00556 
00557 //==============================================================================
00558 int Epetra_CrsGraph::RemoveMyIndices(int Row)
00559 {
00560   int ierr = 0;
00561 
00562   if(IndicesAreContiguous() || StorageOptimized()) 
00563     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00564 
00565   if(IndicesAreGlobal()) 
00566     EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
00567 
00568   if(CrsGraphData_->CV_ == View) 
00569     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00570 
00571   if(Row < 0 || Row >= NumMyBlockRows()) 
00572     EPETRA_CHK_ERR(-1); // Not in Row range
00573     
00574   if (CrsGraphData_->StaticProfile_) {
00575     int NumIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00576     for(int j = 0; j < NumIndices; j++) 
00577       CrsGraphData_->Indices_[Row][j] = -1; // Set to invalid 
00578   }
00579   else
00580     CrsGraphData_->SortedEntries_[Row].entries_.resize(0);
00581 
00582   CrsGraphData_->NumIndicesPerRow_[Row] = 0;
00583 
00584   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00585   EPETRA_CHK_ERR(ierr);
00586 
00587   if(CrsGraphData_->ReferenceCount() > 1)
00588     return(1);
00589   else
00590     return(0);
00591 }
00592 
00593 // protected ===================================================================
00594 bool Epetra_CrsGraph::FindGlobalIndexLoc(int LocalRow,
00595            int Index,
00596            int Start,
00597            int& Loc) const
00598 {
00599   int NumIndices = NumMyIndices(LocalRow);
00600   int* locIndices = Indices(LocalRow);
00601 
00602   // If we have transformed the column indices, we must map this global Index to local
00603   if(CrsGraphData_->IndicesAreLocal_) {
00604     Index = LCID(Index);
00605   }
00606 
00607   if (CrsGraphData_->Sorted_) {
00608     int insertPoint;
00609     Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
00610     return( Loc > -1 );
00611   }
00612   else {
00613     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00614     for(j = 0; j < NumIndices; j++) {
00615       if(j0 >= NumIndices) 
00616   j0 = 0; // wrap around
00617       if(locIndices[j0] == Index) {
00618   Loc = j0;
00619   return(true);
00620       }
00621       j0++;
00622     }
00623   }
00624   return(false);
00625 }
00626 
00627 // protected ===================================================================
00628 bool Epetra_CrsGraph::FindGlobalIndexLoc(int NumIndices,
00629            const int* Indices,
00630            int Index,
00631            int Start,
00632            int& Loc) const
00633 {
00634   // If we have transformed the column indices, we must map this global Index to local
00635   if(CrsGraphData_->IndicesAreLocal_) {
00636     Index = LCID(Index);
00637   }
00638 
00639   if (CrsGraphData_->Sorted_) {
00640     int insertPoint;
00641     Loc = Epetra_Util_binary_search(Index, Indices, NumIndices, insertPoint);
00642     return( Loc > -1 );
00643   }
00644   else {
00645     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00646     for(j = 0; j < NumIndices; j++) {
00647       if(j0 >= NumIndices) 
00648   j0 = 0; // wrap around
00649       if(Indices[j0] == Index) {
00650   Loc = j0;
00651   return(true);
00652       }
00653       j0++;
00654     }
00655   }
00656   return(false);
00657 }
00658 
00659 // protected ===================================================================
00660 bool Epetra_CrsGraph::FindMyIndexLoc(int LocalRow,
00661              int Index,
00662              int Start,
00663              int& Loc) const
00664 {
00665   int NumIndices = NumMyIndices(LocalRow);
00666   int* locIndices = Indices(LocalRow);
00667 
00668   if(!CrsGraphData_->IndicesAreLocal_) {
00669     throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
00670   }
00671 
00672   if (CrsGraphData_->Sorted_) {
00673     int insertPoint;
00674     Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
00675     return( Loc > -1 );
00676   }
00677   else {
00678     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00679     for(j = 0; j < NumIndices; j++) {
00680       if(j0 >= NumIndices) 
00681   j0 = 0; // wrap around
00682       if(locIndices[j0] == Index) {
00683   Loc = j0;
00684   return(true);
00685       }
00686       j0++;
00687     }
00688   }
00689   return(false);
00690 }
00691 
00692 // protected ===================================================================
00693 bool Epetra_CrsGraph::FindMyIndexLoc(int NumIndices,
00694              const int* Indices,
00695              int Index,
00696              int Start,
00697              int& Loc) const
00698 {
00699   if(!CrsGraphData_->IndicesAreLocal_) {
00700     throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
00701   }
00702 
00703   if (CrsGraphData_->Sorted_) {
00704     int insertPoint;
00705     Loc = Epetra_Util_binary_search(Index, Indices, NumIndices, insertPoint);
00706     return( Loc > -1 );
00707   }
00708   else {
00709     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00710     for(j = 0; j < NumIndices; j++) {
00711       if(j0 >= NumIndices) 
00712   j0 = 0; // wrap around
00713       if(Indices[j0] == Index) {
00714   Loc = j0;
00715   return(true);
00716       }
00717       j0++;
00718     }
00719   }
00720   return(false);
00721 }
00722 
00723 //==============================================================================
00724 int Epetra_CrsGraph::FillComplete() {
00725   EPETRA_CHK_ERR(FillComplete(RowMap(), RowMap()));
00726   return(0); // uses FillComplete(...)'s shared-ownership test. 
00727 }
00728 
00729 //==============================================================================
00730 int Epetra_CrsGraph::FillComplete(const Epetra_BlockMap& DomainMap, const Epetra_BlockMap& RangeMap) {
00731   CrsGraphData_->DomainMap_ = DomainMap;
00732   CrsGraphData_->RangeMap_ = RangeMap;
00733 
00734   MakeIndicesLocal(DomainMap, RangeMap); // Convert indices to zero based on each processor
00735   SortIndices();  // Sort column entries from smallest to largest
00736   RemoveRedundantIndices(); // Get rid of any redundant index values
00737   DetermineTriangular(); //determine if lower/upper triangular
00738   CrsGraphData_->MakeImportExport(); // Build Import or Export objects
00739   ComputeGlobalConstants(); // Compute constants that require communication
00740   SetFilled(true);
00741 
00742   if(CrsGraphData_->ReferenceCount() > 1)
00743     return(1);
00744   else
00745     return(0);
00746 }
00747 
00748 //==============================================================================
00749 int Epetra_CrsGraph::TransformToLocal() {
00750   return(FillComplete());
00751 }
00752 
00753 //==============================================================================
00754 int Epetra_CrsGraph::TransformToLocal(const Epetra_BlockMap* DomainMap, const Epetra_BlockMap* RangeMap) {
00755   return(FillComplete(*DomainMap, *RangeMap));
00756 }
00757 
00758 // private =====================================================================
00759 int Epetra_CrsGraph::ComputeGlobalConstants()
00760 {
00761   if(GlobalConstantsComputed()) 
00762     return(0);
00763 
00764   Epetra_IntSerialDenseVector tempvec(8); // Temp space
00765 
00766   const int numMyBlockRows = NumMyBlockRows();
00767 
00768   CrsGraphData_->NumMyEntries_ = 0; // Compute Number of Nonzero entries and max
00769   CrsGraphData_->MaxNumIndices_ = 0;
00770   {for(int i = 0; i < numMyBlockRows; i++) {
00771     CrsGraphData_->NumMyEntries_ += CrsGraphData_->NumIndicesPerRow_[i];
00772     CrsGraphData_->MaxNumIndices_ = EPETRA_MAX(CrsGraphData_->MaxNumIndices_,CrsGraphData_->NumIndicesPerRow_[i]);
00773   }}
00774   
00775   // Case 1:  Constant block size (including blocksize = 1)
00776   if(RowMap().ConstantElementSize()) {
00777     tempvec[0] = CrsGraphData_->NumMyEntries_;
00778     tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
00779 
00780     Comm().SumAll(&tempvec[0], &tempvec[2], 2);
00781     Comm().MaxAll(&CrsGraphData_->MaxNumIndices_, &CrsGraphData_->GlobalMaxNumIndices_, 1);
00782     
00783     CrsGraphData_->NumGlobalEntries_ = tempvec[2];
00784     CrsGraphData_->NumGlobalBlockDiagonals_ = tempvec[3];
00785 
00786     int RowElementSize = RowMap().MaxElementSize();
00787     int ColElementSize = RowElementSize;
00788     CrsGraphData_->NumGlobalDiagonals_ = tempvec[3] * RowElementSize;
00789     CrsGraphData_->NumMyNonzeros_ = CrsGraphData_->NumMyEntries_ * RowElementSize * ColElementSize;
00790     CrsGraphData_->NumGlobalNonzeros_ = CrsGraphData_->NumGlobalEntries_ * RowElementSize * ColElementSize;
00791     CrsGraphData_->MaxNumNonzeros_ = CrsGraphData_->MaxNumIndices_ * RowElementSize * ColElementSize;
00792     CrsGraphData_->GlobalMaxNumNonzeros_ = CrsGraphData_->GlobalMaxNumIndices_ * RowElementSize * ColElementSize;
00793   }
00794 
00795   // Case 2:  Variable block size (more work)
00796   else {
00797     CrsGraphData_->NumMyNonzeros_ = 0;  // We will count the number of nonzeros here
00798     CrsGraphData_->MaxNumNonzeros_ = 0;  // We will determine the max number of nonzeros in any one block row
00799     int* RowElementSizeList = RowMap().ElementSizeList();
00800     int* ColElementSizeList = RowElementSizeList;
00801     if(Importer() != 0) 
00802       ColElementSizeList = ColMap().ElementSizeList();
00803     for(int i = 0; i < numMyBlockRows; i++){
00804       int NumEntries = CrsGraphData_->NumIndicesPerRow_[i];
00805       int* Indices = CrsGraphData_->Indices_[i];
00806       if(NumEntries > 0) {
00807   int CurNumNonzeros = 0;
00808   int RowDim = RowElementSizeList[i];
00809   for(int j = 0; j < NumEntries; j++) {
00810     int ColDim = ColElementSizeList[Indices[j]];
00811     CurNumNonzeros += RowDim*ColDim;
00812     CrsGraphData_->MaxColDim_ = EPETRA_MAX(CrsGraphData_->MaxColDim_, ColDim);
00813   }
00814   CrsGraphData_->MaxNumNonzeros_ = EPETRA_MAX(CrsGraphData_->MaxNumNonzeros_, CurNumNonzeros);
00815   CrsGraphData_->NumMyNonzeros_ += CurNumNonzeros;
00816       }
00817     }
00818     
00819     // Sum Up all nonzeros
00820     
00821     tempvec[0] = CrsGraphData_->NumMyEntries_;
00822     tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
00823     tempvec[2] = CrsGraphData_->NumMyDiagonals_;
00824     tempvec[3] = CrsGraphData_->NumMyNonzeros_;
00825     
00826     Comm().SumAll(&tempvec[0], &tempvec[4], 4);
00827     
00828     CrsGraphData_->NumGlobalEntries_ = tempvec[4];
00829     CrsGraphData_->NumGlobalBlockDiagonals_ = tempvec[5];
00830     CrsGraphData_->NumGlobalDiagonals_ = tempvec[6];
00831     CrsGraphData_->NumGlobalNonzeros_ = tempvec[7];
00832 
00833     tempvec[0] = CrsGraphData_->MaxNumIndices_;
00834     tempvec[1] = CrsGraphData_->MaxNumNonzeros_;
00835 
00836     Comm().MaxAll(&tempvec[0], &tempvec[2], 2);
00837 
00838     CrsGraphData_->GlobalMaxNumIndices_ = tempvec[2];
00839     CrsGraphData_->GlobalMaxNumNonzeros_ = tempvec[3];
00840   }
00841   
00842   CrsGraphData_->NumGlobalRows_ = CrsGraphData_->RangeMap_.NumGlobalPoints();
00843   CrsGraphData_->NumGlobalCols_ = DomainMap().NumGlobalPoints();
00844 
00845   CrsGraphData_->GlobalConstantsComputed_ = true;
00846   
00847   return(0);
00848 }
00849 
00850 void epetra_shellsort(int* list, int length)
00851 {
00852   int i, j, j2, temp, istep;
00853   unsigned step;
00854 
00855   step = length/2;
00856   while (step > 0)
00857   {
00858     for (i=step; i < length; i++)
00859     {
00860       istep = step;
00861       j = i;
00862       j2 = j-istep;
00863       temp = list[i];
00864       if (list[j2] > temp) {
00865         while ((j >= istep) && (list[j2] > temp))
00866         {
00867           list[j] = list[j2];
00868           j = j2;
00869           j2 -= istep;
00870         }
00871         list[j] = temp;
00872       }
00873     }
00874 
00875     if (step == 2)
00876       step = 1;
00877     else
00878       step = (int) (step / 2.2);
00879   }
00880 }
00881 
00882 //==============================================================================
00883 int Epetra_CrsGraph::SortIndices() {
00884   if(IndicesAreGlobal()) 
00885     EPETRA_CHK_ERR(-1);
00886   if(Sorted())
00887     return(0);
00888 
00889   // For each row, sort column entries from smallest to largest.
00890   // Use shell sort, which is fast if indices are already sorted.
00891 
00892   const int numMyBlockRows = NumMyBlockRows();
00893   for(int i = 0; i < numMyBlockRows; i++){
00894     int n = CrsGraphData_->NumIndicesPerRow_[i];
00895     int* const list = CrsGraphData_->Indices_[i];
00896 
00897     epetra_shellsort(list, n);
00898 //    int m = n/2;
00899     
00900 //    while(m > 0) {
00901  //     int max = n - m;
00902 //      for(int j = 0; j < max; j++) {
00903 //        int k = j;
00904 //        while(k>-1) {
00905 //    if(list[k+m] >= list[k])
00906 //      break;
00907 //    int itemp = list[k+m];
00908 //    list[k+m] = list[k];
00909 //    list[k] = itemp;
00910 //          k-=m;
00911 //  }
00912 //      }
00913 //      m = m/2;
00914 //    }
00915   }
00916   SetSorted(true);
00917 
00918   if(CrsGraphData_->ReferenceCount() > 1)
00919     return(1);
00920   else
00921     return(0);
00922 }
00923 
00924 void epetra_crsgraph_compress_out_duplicates(int len, int* list, int& newlen)
00925 {
00926   //
00927   //This function runs the array ('list') checking for
00928   //duplicate entries. Any duplicates that are found are
00929   //removed by sliding subsequent data down in the array,
00930   //over-writing the duplicates. Finally, the new length
00931   //of the array (i.e., the number of unique entries) is
00932   //placed in the output parameter 'newlen'. The array is
00933   //**not** re-allocated.
00934   //
00936   //Important assumption: The array contents are assumed to
00937   //be sorted before this function is called. If the array
00938   //contents are not sorted, then the behavior of this
00939   //function is undefined.
00941   //
00942 
00943   if (len < 2) return;
00944 
00945   int* ptr0 = &list[0];
00946   int* ptr1 = &list[1];
00947 
00948   int* ptr_end = &list[len-1];
00949 
00950   while(*ptr0 != *ptr1 && ptr1 < ptr_end) {
00951     ++ptr0;
00952     ++ptr1;
00953   }
00954 
00955   if (ptr1 < ptr_end) {
00956     //if ptr1 < ptr_end we've found a duplicate...
00957 
00958     ++ptr0;
00959     ++ptr1;
00960 
00961     while(*ptr0 == *ptr1 && ptr1 < ptr_end) ++ptr1;
00962 
00963     while(ptr1 < ptr_end) {
00964 
00965       int val = *ptr1++;
00966 
00967       while(val == *ptr1 && ptr1 < ptr_end) {
00968         ++ptr1;
00969       }
00970 
00971       *ptr0++ = val;
00972     }
00973 
00974     if (*(ptr0-1) != *ptr1) *ptr0++ = *ptr1;
00975 
00976     int num_removed = (int)(ptr_end - ptr0 + 1);
00977     newlen = len - num_removed;
00978   }
00979   else {
00980     if (*ptr0 == *ptr1) newlen = len - 1;
00981     else newlen = len;
00982   }
00983 }
00984 
00985 //==============================================================================
00986 int Epetra_CrsGraph::RemoveRedundantIndices()
00987 {
00988   if(!Sorted())
00989     EPETRA_CHK_ERR(-1);  // Must have sorted index set
00990   if(IndicesAreGlobal()) 
00991     EPETRA_CHK_ERR(-2); // Indices must be local
00992 
00993   // Note:  This function assumes that SortIndices was already called.
00994   // For each row, remove column indices that are repeated.
00995 
00996   const int numMyBlockRows = NumMyBlockRows();
00997 
00998   if(NoRedundancies() == false) {
00999     int* numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
01000     for(int i=0; i<numMyBlockRows; ++i) {
01001       int NumIndices = numIndicesPerRow[i];
01002       int* col_indices = this->Indices(i);
01003   
01004       if(NumIndices > 1) {
01005         epetra_crsgraph_compress_out_duplicates(NumIndices, col_indices,
01006                                                 numIndicesPerRow[i]);
01007       }
01008       // update vector size and address in memory
01009       if (!CrsGraphData_->StaticProfile_) {
01010         CrsGraphData_->SortedEntries_[i].entries_.assign(col_indices, col_indices+numIndicesPerRow[i]);
01011         if (numIndicesPerRow[i] > 0) {
01012           CrsGraphData_->Indices_[i] = &(CrsGraphData_->SortedEntries_[i].entries_[0]);
01013         }
01014         else {
01015           CrsGraphData_->Indices_[i] = NULL;
01016         }
01017       }
01018     }
01019   }
01020 
01021   SetNoRedundancies(true);
01022   return 0;
01023 }
01024 
01025 int Epetra_CrsGraph::DetermineTriangular()
01026 {
01027   // determine if graph is upper or lower triangular or has no diagonal
01028   
01029   if(!Sorted())
01030     EPETRA_CHK_ERR(-1);  // Must have sorted index set
01031   if(IndicesAreGlobal()) 
01032     EPETRA_CHK_ERR(-2); // Indices must be local
01033 
01034   CrsGraphData_->NumMyDiagonals_ = 0;
01035   CrsGraphData_->NumMyBlockDiagonals_ = 0;
01036 
01037   const Epetra_BlockMap& rowMap = RowMap();
01038   const Epetra_BlockMap& colMap = ColMap();
01039 
01040   const int numMyBlockRows = NumMyBlockRows();
01041 
01042   for(int i = 0; i < numMyBlockRows; i++) {
01043     int NumIndices = NumMyIndices(i);
01044     if(NumIndices > 0) {
01045       int ig = rowMap.GID(i);
01046       int* col_indices = this->Indices(i);
01047 
01048       int jl_0 = col_indices[0];
01049       int jl_n = col_indices[NumIndices-1];
01050 
01051       if(jl_n > i) CrsGraphData_->LowerTriangular_ = false;
01052       if(jl_0 < i) CrsGraphData_->UpperTriangular_ = false;
01053 
01054       //jl will be the local-index for the diagonal that we
01055       //want to search for.
01056       int jl = colMap.LID(ig);
01057 
01058       int insertPoint = -1;
01059       if (Epetra_Util_binary_search(jl, col_indices, NumIndices, insertPoint)>-1) {
01060         CrsGraphData_->NumMyBlockDiagonals_++;
01061         CrsGraphData_->NumMyDiagonals_ += rowMap.ElementSize(i);
01062       }
01063     }
01064   }
01065 
01066   CrsGraphData_->NoDiagonal_ = (CrsGraphData_->NumMyBlockDiagonals_ == 0);
01067 
01068   if(CrsGraphData_->ReferenceCount() > 1)
01069     return(1);
01070   else
01071     return(0);
01072 }
01073 
01074 // private =====================================================================
01075 int Epetra_CrsGraph::MakeColMap(const Epetra_BlockMap& DomainMap,
01076         const Epetra_BlockMap& RangeMap)
01077 {
01078   (void)RangeMap;
01079   int i;
01080   int j;
01081 
01082   if(CrsGraphData_->HaveColMap_) 
01083     return(0); // Already have a Column Map
01084 
01085   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01086   if(IndicesAreLocal()) 
01087     EPETRA_CHK_ERR(-1); // Return error: Indices must be global
01088   
01089   // Scan all column indices and sort into two groups: 
01090   // Local:  those whose GID matches a GID of the domain map on this processor and
01091   // Remote: All others.
01092   int numDomainElements = DomainMap.NumMyElements();
01093   bool * LocalGIDs  = 0;
01094   if (numDomainElements>0) LocalGIDs  = new bool[numDomainElements];
01095   for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
01096 
01097   // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
01098   // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
01099   // and the number of block rows.
01100   const int numMyBlockRows = NumMyBlockRows();
01101   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
01102   //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << endl;
01103   Epetra_HashTable RemoteGIDs(hashsize); 
01104   Epetra_HashTable RemoteGIDList(hashsize);
01105 
01106   int NumLocalColGIDs = 0;
01107   int NumRemoteColGIDs = 0;
01108   for(i = 0; i < numMyBlockRows; i++) {
01109     const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01110     int* ColIndices = CrsGraphData_->Indices_[i];
01111     for(j = 0; j < NumIndices; j++) {
01112       int GID = ColIndices[j];
01113       // Check if GID matches a row GID
01114       int LID = DomainMap.LID(GID);
01115       if(LID != -1) {
01116   bool alreadyFound = LocalGIDs[LID];
01117   if (!alreadyFound) {
01118           LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
01119           NumLocalColGIDs++;
01120   }
01121       }
01122       else {
01123   if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
01124     RemoteGIDs.Add(GID, NumRemoteColGIDs);
01125     RemoteGIDList.Add(NumRemoteColGIDs++, GID);
01126   }
01127       }
01128     }
01129   }
01130 
01131   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=DomainMap and quit
01132   if (DomainMap.Comm().NumProc()==1) { 
01133     
01134     if (NumRemoteColGIDs!=0) {
01135       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
01136     }
01137     if (NumLocalColGIDs==numDomainElements) {
01138       CrsGraphData_->ColMap_ = DomainMap;
01139       CrsGraphData_->HaveColMap_ = true;
01140       if (LocalGIDs!=0) delete [] LocalGIDs; 
01141       return(0); 
01142     }
01143   }
01144       
01145   // Now build integer array containing column GIDs
01146   // Build back end, containing remote GIDs, first
01147   int NumMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
01148   Epetra_IntSerialDenseVector ColIndices;
01149   if(NumMyBlockCols > 0) 
01150     ColIndices.Size(NumMyBlockCols);
01151 
01152   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01153 
01154   for(i = 0; i < NumRemoteColGIDs; i++) 
01155     RemoteColIndices[i] = RemoteGIDList.Get(i); 
01156 
01157   int NLists = 1;
01158   Epetra_IntSerialDenseVector PIDList;
01159   Epetra_IntSerialDenseVector SizeList;
01160   int* RemoteSizeList = 0;
01161   bool DoSizes = !DomainMap.ConstantElementSize(); // If not constant element size, then we must exchange
01162       
01163   if(NumRemoteColGIDs > 0) 
01164     PIDList.Size(NumRemoteColGIDs);
01165 
01166   if(DoSizes) {
01167     if(NumMyBlockCols > 0) 
01168       SizeList.Size(NumMyBlockCols);
01169     RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
01170     NLists++;
01171   }
01172   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01173       
01174   // Sort External column indices so that all columns coming from a given remote processor are contiguous
01175 
01176   Epetra_Util Util;
01177   int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
01178   SortLists[0] = RemoteColIndices;
01179   SortLists[1] = RemoteSizeList;
01180   Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
01181   if (CrsGraphData_->SortGhostsAssociatedWithEachProcessor_) {
01182     // Sort external column indices so that columns from a given remote processor are not only contiguous
01183     // but also in ascending order. NOTE: I don't know if the number of externals associated
01184     // with a given remote processor is known at this point ... so I count them here.
01185 
01186     NLists--;
01187     int StartCurrent, StartNext;
01188     StartCurrent = 0; StartNext = 1;
01189     while ( StartNext < NumRemoteColGIDs ) {
01190       if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
01191       else {
01192         if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01193         Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists);
01194         StartCurrent = StartNext; StartNext++;
01195       }
01196     }
01197     if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01198     Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists);
01199   }
01200 
01201   // Now fill front end. Two cases:
01202   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
01203   //     can simply read the domain GIDs into the front part of ColIndices, otherwise 
01204   // (2) We step through the GIDs of the DomainMap, checking to see if each domain GID is a column GID.
01205   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
01206 
01207   if(NumLocalColGIDs == DomainMap.NumMyElements()) {
01208     DomainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first NumMyBlockCols elements column GID list
01209     if(DoSizes) 
01210       DomainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
01211   }
01212   else {
01213     int NumMyElements = DomainMap.NumMyElements();
01214     int* MyGlobalElements = DomainMap.MyGlobalElements();
01215     int* ElementSizeList = 0;
01216     if(DoSizes) 
01217       ElementSizeList = DomainMap.ElementSizeList();
01218     int NumLocalAgain = 0;
01219     for(i = 0; i < NumMyElements; i++) {
01220       if(LocalGIDs[i]) {
01221   if(DoSizes) 
01222     SizeList[NumLocalAgain] = ElementSizeList[i];
01223   ColIndices[NumLocalAgain++] = MyGlobalElements[i];
01224       }
01225     }
01226     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
01227   }
01228 
01229   // Done with this array
01230   if (LocalGIDs!=0) delete [] LocalGIDs; 
01231 
01232 
01233   // Make Column map with same element sizes as Domain map
01234 
01235   if(DomainMap.MaxElementSize() == 1) { // Simple map
01236     Epetra_Map temp(-1, NumMyBlockCols, ColIndices.Values(), DomainMap.IndexBase(), DomainMap.Comm());
01237     CrsGraphData_->ColMap_ = temp;
01238   }
01239   else if(DomainMap.ConstantElementSize()) { // Constant Block size map
01240     Epetra_BlockMap temp(-1, NumMyBlockCols, ColIndices.Values(), DomainMap.MaxElementSize(),DomainMap.IndexBase(), DomainMap.Comm());
01241     CrsGraphData_->ColMap_ = temp;
01242   }
01243   else { // Most general case where block size is variable.
01244     Epetra_BlockMap temp(-1, NumMyBlockCols, ColIndices.Values(), SizeList.Values(), DomainMap.IndexBase(), DomainMap.Comm());
01245     CrsGraphData_->ColMap_ = temp;
01246   }
01247   CrsGraphData_->HaveColMap_ = true;
01248 
01249   return(0);
01250 }
01251 
01252 // protected ===================================================================
01253 int Epetra_CrsGraph::MakeIndicesLocal(const Epetra_BlockMap& DomainMap, const Epetra_BlockMap& RangeMap) {
01254   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01255   if(IndicesAreLocal() && IndicesAreGlobal()) 
01256     EPETRA_CHK_ERR(-1); // Return error: Indices must not be both local and global
01257 
01258   MakeColMap(DomainMap, RangeMap); // If user has not prescribed column map, create one from indices
01259   const Epetra_BlockMap& colmap = ColMap();
01260 
01261   // Store number of local columns
01262   CrsGraphData_->NumMyCols_ = ColMap().NumMyPoints();
01263   CrsGraphData_->NumMyBlockCols_ = ColMap().NumMyElements();
01264 
01265   // Transform indices to local index space
01266   const int numMyBlockRows = NumMyBlockRows();
01267 
01268   if(IndicesAreGlobal()) {
01269     // Check if ColMap is monotone. If not, the list will get unsorted.
01270     bool mapMonotone = true;
01271     {
01272       int oldGID = colmap.GID(0);
01273       for (int i=1; i<colmap.NumMyElements(); ++i) {
01274         if (oldGID > colmap.GID(i)) {
01275           mapMonotone = false;
01276           break;
01277         }
01278         oldGID = colmap.GID(i);
01279       }
01280     }
01281     if (Sorted())
01282       SetSorted(mapMonotone);
01283 
01284     // now comes the actual transformation
01285     for(int i = 0; i < numMyBlockRows; i++) {
01286       const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01287       int* ColIndices = CrsGraphData_->Indices_[i];
01288       for(int j = 0; j < NumIndices; j++) {
01289         int GID = ColIndices[j];
01290         int LID = colmap.LID(GID);
01291         if(LID != -1) 
01292           ColIndices[j] = LID;
01293         else 
01294           throw ReportError("Internal error in FillComplete ",-1); 
01295       }
01296     }
01297   }
01298   
01299   SetIndicesAreLocal(true);
01300   SetIndicesAreGlobal(false);
01301   
01302   if(CrsGraphData_->ReferenceCount() > 1)
01303     return(1); // return 1 if data is shared
01304   else
01305     return(0);
01306 }
01307 
01308 //==============================================================================
01309 int Epetra_CrsGraph::OptimizeStorage() {
01310   int NumIndices;
01311   const int numMyBlockRows = NumMyBlockRows();
01312 
01313   if(StorageOptimized()) 
01314     return(0); // Have we been here before?
01315   if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
01316 
01317   bool Contiguous = true; // Assume contiguous is true
01318   for(int i = 1; i < numMyBlockRows; i++) {
01319     NumIndices = CrsGraphData_->NumIndicesPerRow_[i-1];
01320     int NumAllocateIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[i-1];
01321 
01322     // Check if NumIndices is same as NumAllocatedIndices and 
01323     // check if end of beginning of current row starts immediately after end of previous row.
01324     if((NumIndices != NumAllocateIndices) || 
01325        (CrsGraphData_->Indices_[i] != CrsGraphData_->Indices_[i-1] + NumIndices)) {
01326       Contiguous = false;
01327       break;
01328     }
01329   }
01330 
01331   // NOTE:  At the end of the above loop set, there is a possibility that NumIndices and NumAllocatedIndices
01332   //        for the last row could be different, but I don't think it matters.
01333 
01334 
01335   if((CrsGraphData_->CV_ == View) && !Contiguous) 
01336     return(3);  // This is user data, it's not contiguous and we can't make it so.
01337 
01338   if(CrsGraphData_->IndexOffset_ .Values() != CrsGraphData_->NumIndicesPerRow_.Values())
01339     CrsGraphData_->IndexOffset_.MakeViewOf(CrsGraphData_->NumIndicesPerRow_);
01340 
01341   // This next step constructs the scan sum of the number of indices per row.  Note that the
01342   // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
01343   // careful with this sum operation
01344   int * NumIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
01345   int curNumIndices = NumIndicesPerRow[0];
01346   NumIndicesPerRow[0] = 0;
01347   for (int i=0; i<numMyBlockRows; ++i) {
01348     int nextNumIndices = NumIndicesPerRow[i+1];
01349     NumIndicesPerRow[i+1] = NumIndicesPerRow[i]+curNumIndices;
01350     curNumIndices = nextNumIndices;
01351   }
01352 
01353   if(!Contiguous) { // Must pack indices if not already contiguous
01354 
01355     // Allocate one big integer array for all index values
01356     if (!(StaticProfile())) { // If static profile, All_Indices_ is already allocated, only need to pack data
01357       int errorcode = CrsGraphData_->All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
01358       if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
01359     }
01360       // Pack indices into All_Indices_
01361 
01362     int* All_Indices = CrsGraphData_->All_Indices_.Values();
01363     int * IndexOffset = CrsGraphData_->IndexOffset_.Values();
01364     int ** Indices = CrsGraphData_->Indices_;
01365     
01366     if (!(StaticProfile())) {
01367 #ifdef Epetra_HAVE_OMP
01368 #pragma omp parallel for default(none) shared(IndexOffset,All_Indices,Indices)
01369 #endif   
01370       for(int i = 0; i < numMyBlockRows; i++) {
01371         int curNumIndices = IndexOffset[i+1] - IndexOffset[i];
01372         int* ColIndices = Indices[i];
01373         int *newColIndices = All_Indices+IndexOffset[i];
01374         for(int j = 0; j < curNumIndices; j++) newColIndices[j] = ColIndices[j];
01375       }
01376       for(int i = 0; i < numMyBlockRows; i++) {
01377         if (Indices[i]!=0) {
01378           CrsGraphData_->SortedEntries_[i].entries_.clear();
01379           Indices[i] = 0;
01380         }
01381      }
01382    } // End of contiguous non-static section   
01383    else {
01384 
01385      for(int i = 0; i < numMyBlockRows; i++) {
01386        int curNumIndices = IndexOffset[i+1] - IndexOffset[i];
01387        int* ColIndices = Indices[i];
01388        int *newColIndices = All_Indices+IndexOffset[i];
01389        if (ColIndices!=newColIndices) // No need to copy if pointing to same space
01390        for(int j = 0; j < curNumIndices; j++) newColIndices[j] = ColIndices[j];
01391        Indices[i] = 0;
01392      }
01393    } // End of Contiguous static section
01394   } // End of !Contiguous section
01395   else {
01396     //if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
01397     const int numMyNonzeros = NumMyNonzeros();
01398     if (numMyBlockRows > 0 && !(StaticProfile())) {
01399       int errorcode = CrsGraphData_->All_Indices_.Size(numMyNonzeros);
01400       if(errorcode != 0)  throw ReportError("Error with All_Indices_ allocation.", -99);
01401       int* all_indices_values = CrsGraphData_->All_Indices_.Values();
01402       int* indices_values = CrsGraphData_->Indices_[0];
01403 #ifdef Epetra_HAVE_OMP
01404 #pragma omp parallel for default(none) shared(all_indices_values,indices_values)
01405 #endif   
01406       for(int ii=0; ii<numMyNonzeros; ++ii) {
01407         all_indices_values[ii] = indices_values[ii];
01408       }
01409     }
01410   }
01411 
01412   // Delete unneeded storage
01413   CrsGraphData_->NumAllocatedIndicesPerRow_.Resize(0);
01414   delete [] CrsGraphData_->Indices_; CrsGraphData_->Indices_=0;
01415   CrsGraphData_->SortedEntries_.clear();
01416 
01417   SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
01418   CrsGraphData_->StorageOptimized_ = true;
01419 
01420   return(0);
01421 }
01422 
01423 //==============================================================================
01424 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01425 {
01426   int j;
01427 
01428   Row = LRID(Row); // Normalize row range
01429 
01430   if(Row < 0 || Row >= NumMyBlockRows()) 
01431     EPETRA_CHK_ERR(-1); // Not in Row range
01432 
01433   NumIndices = NumMyIndices(Row);
01434   if(LenOfIndices < NumIndices) 
01435     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01436 
01437   int * srcIndices = Indices(Row);
01438   if(IndicesAreLocal())  
01439     for(j = 0; j < NumIndices; j++) 
01440       targIndices[j] = GCID(srcIndices[j]);
01441   else 
01442     for(j = 0; j < NumIndices; j++)
01443       targIndices[j] = srcIndices[j];
01444   
01445   return(0);
01446 }
01447 
01448 //==============================================================================
01449 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01450 {
01451   int j;
01452 
01453   if(Row < 0 || Row >= NumMyBlockRows()) 
01454     EPETRA_CHK_ERR(-1); // Not in Row range
01455 
01456   NumIndices = NumMyIndices(Row);
01457   if(LenOfIndices < NumIndices) 
01458     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01459 
01460   if(IndicesAreGlobal()) 
01461     EPETRA_CHK_ERR(-3); // There are no local indices yet
01462 
01463   int * srcIndices = Indices(Row);
01464   for(j = 0; j < NumIndices; j++)
01465     targIndices[j] = srcIndices[j];
01466   
01467   return(0);
01468 }
01469 
01470 //==============================================================================
01471 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const 
01472 {
01473   Row = LRID(Row); // Normalize row range
01474 
01475   if(Row < 0 || Row >= NumMyBlockRows()) 
01476     EPETRA_CHK_ERR(-1); // Not in Row range
01477 
01478   if(IndicesAreLocal()) 
01479     EPETRA_CHK_ERR(-2); // There are no global indices
01480 
01481   NumIndices = NumMyIndices(Row);
01482 
01483   targIndices = Indices(Row);
01484   
01485   return(0);
01486 }
01487 
01488 //==============================================================================
01489 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const 
01490 {
01491   if(Row < 0 || Row >= NumMyBlockRows()) 
01492     EPETRA_CHK_ERR(-1); // Not in Row range
01493 
01494   if(IndicesAreGlobal()) 
01495     EPETRA_CHK_ERR(-2); // There are no local indices
01496 
01497   NumIndices = NumMyIndices(Row);
01498 
01499   targIndices = Indices(Row);
01500   
01501   return(0);
01502 }
01503 
01504 //==============================================================================
01505 int Epetra_CrsGraph::NumGlobalIndices(int Row) const {
01506   Row = LRID(Row);
01507   if(Row != -1) 
01508     return(NumMyIndices(Row));
01509   else 
01510     return(0); // No indices for this row on this processor
01511 }
01512 
01513 //==============================================================================
01514 int Epetra_CrsGraph::NumAllocatedGlobalIndices(int Row) const {
01515   Row = LRID(Row);
01516   if(Row != -1) 
01517     return(NumAllocatedMyIndices(Row));
01518   else 
01519     return(0); // No indices allocated for this row on this processor
01520 }
01521 
01522 //==============================================================================
01523 int Epetra_CrsGraph::ReplaceRowMap(const Epetra_BlockMap& newmap)
01524 {
01525   if (RowMap().PointSameAs(newmap)) {
01526     Epetra_DistObject::Map_ = newmap;
01527     CrsGraphData_->RowMap_ = newmap;
01528     CrsGraphData_->MakeImportExport();
01529     return(0);
01530   }
01531 
01532   return(-1);
01533 }
01534 
01535 //==============================================================================
01536 int Epetra_CrsGraph::ReplaceColMap(const Epetra_BlockMap& newmap)
01537 {
01538   if (ColMap().PointSameAs(newmap)) {
01539     CrsGraphData_->ColMap_ = newmap;
01540     CrsGraphData_->MakeImportExport();
01541     return(0);
01542   }
01543 
01544   return(-1);
01545 }
01546 
01547 // private =====================================================================
01548 int Epetra_CrsGraph::CheckSizes(const Epetra_SrcDistObject& Source) {
01549   try {
01550     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
01551     if(!A.GlobalConstantsComputed()) 
01552       EPETRA_CHK_ERR(-1); // Must have global constants to proceed
01553   }
01554   catch(...) {
01555     return(0); // No error at this point, object could be a RowMatrix
01556   }
01557   return(0);
01558 }
01559 
01560 // private =====================================================================
01561 int Epetra_CrsGraph::CopyAndPermute(const Epetra_SrcDistObject& Source,
01562             int NumSameIDs, 
01563             int NumPermuteIDs, 
01564             int* PermuteToLIDs,
01565             int* PermuteFromLIDs,
01566                                     const Epetra_OffsetIndex * Indexor)
01567 { 
01568   try {
01569     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01570     EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01571             PermuteFromLIDs,Indexor));
01572   }
01573   catch(...) {
01574     try {
01575       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01576       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01577                PermuteFromLIDs,Indexor));
01578     }
01579     catch(...) {
01580       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
01581     }
01582   }
01583   
01584   return(0);
01585 }
01586 
01587 // private =====================================================================
01588 int Epetra_CrsGraph::CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
01589                int NumSameIDs, 
01590                int NumPermuteIDs, 
01591                int* PermuteToLIDs,
01592                int* PermuteFromLIDs,
01593                                              const Epetra_OffsetIndex * Indexor)
01594 {
01595   (void)Indexor;
01596   int i;
01597   int j;
01598   int NumIndices;
01599   int FromRow;
01600   int ToRow;
01601   int MaxNumIndices = A.MaxNumEntries();
01602   Epetra_IntSerialDenseVector Indices;
01603   Epetra_SerialDenseVector Values;
01604 
01605   if(MaxNumIndices > 0) {
01606     Indices.Size(MaxNumIndices);
01607     Values.Size(MaxNumIndices); // Must extract values even though we discard them
01608   }
01609 
01610   const Epetra_Map& RowMap = A.RowMatrixRowMap();
01611   const Epetra_Map& ColMap = A.RowMatrixColMap();
01612   
01613   // Do copy first
01614   for(i = 0; i < NumSameIDs; i++) {
01615     ToRow = RowMap.GID(i);
01616     EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, MaxNumIndices, NumIndices, Values.Values(), Indices.Values()));
01617     for(j = 0; j < NumIndices; j++) 
01618       Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01619     // Place into target graph.  
01620     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices.Values());
01621     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01622   }
01623   
01624   // Do local permutation next
01625   for(i = 0; i < NumPermuteIDs; i++) {
01626     FromRow = PermuteFromLIDs[i];
01627     ToRow = GRID(PermuteToLIDs[i]);
01628     EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, MaxNumIndices, NumIndices, Values.Values(), Indices.Values()));
01629     for(j = 0; j < NumIndices; j++) 
01630       Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01631     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices.Values()); // Place into target graph.
01632     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01633   }
01634   
01635   return(0);
01636 }
01637 
01638 // private =====================================================================
01639 int Epetra_CrsGraph::CopyAndPermuteCrsGraph(const Epetra_CrsGraph& A,
01640               int NumSameIDs, 
01641               int NumPermuteIDs, 
01642               int* PermuteToLIDs,
01643               int* PermuteFromLIDs,
01644                                             const Epetra_OffsetIndex * Indexor)
01645 {
01646   (void)Indexor;
01647   int i;
01648   int Row;
01649   int NumIndices;
01650   int* Indices = 0;
01651   int FromRow, ToRow;
01652   int MaxNumIndices = A.MaxNumIndices();
01653   Epetra_IntSerialDenseVector IndicesVector;
01654 
01655   if(MaxNumIndices > 0 && A.IndicesAreLocal()) {
01656     IndicesVector.Size(MaxNumIndices);
01657     Indices = IndicesVector.Values();
01658   }
01659   
01660   // Do copy first
01661   if(NumSameIDs > 0) {
01662     if(A.IndicesAreLocal()) {
01663       for(i = 0; i < NumSameIDs; i++) {
01664         Row = GRID(i);
01665         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumIndices, Indices));
01666         // Place into target graph.  
01667         int ierr = InsertGlobalIndices(Row, NumIndices, Indices); 
01668         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01669       }
01670     }
01671     else { // A.IndiceAreGlobal()
01672       for(i = 0; i < NumSameIDs; i++) {
01673         Row = GRID(i);
01674         EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, Indices));
01675         // Place into target graph.  
01676         int ierr = InsertGlobalIndices(Row, NumIndices, Indices); 
01677         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01678       }
01679     } 
01680   }
01681 
01682   // Do local permutation next
01683   if(NumPermuteIDs > 0) {
01684     if(A.IndicesAreLocal()) {
01685       for(i = 0; i < NumPermuteIDs; i++) {
01686         FromRow = A.GRID(PermuteFromLIDs[i]);
01687         ToRow = GRID(PermuteToLIDs[i]);
01688         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, MaxNumIndices, NumIndices, Indices));
01689         // Place into target graph.
01690         int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices); 
01691         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01692       }
01693     }
01694     else { // A.IndiceAreGlobal()
01695       for(i = 0; i < NumPermuteIDs; i++) {
01696         FromRow = A.GRID(PermuteFromLIDs[i]);
01697         ToRow = GRID(PermuteToLIDs[i]);
01698         EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, Indices));
01699         // Place into target graph.
01700         int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices); 
01701         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01702       }
01703     }
01704   } 
01705   
01706   return(0);
01707 }
01708 
01709 // private =====================================================================
01710 int Epetra_CrsGraph::PackAndPrepare(const Epetra_SrcDistObject& Source, 
01711             int NumExportIDs, 
01712             int* ExportLIDs,
01713             int& LenExports, 
01714             char*& Exports, 
01715             int& SizeOfPacket, 
01716                                     int * Sizes,
01717                                     bool& VarSizes,
01718             Epetra_Distributor& Distor) 
01719 {
01720   int GlobalMaxNumIndices = 0;
01721   int TotalSendSize = 0;
01722 
01723   VarSizes = true;
01724 
01725   SizeOfPacket = (int)sizeof(int); 
01726 
01727   if(NumExportIDs <= 0) return(0);
01728 
01729   try {
01730     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01731     GlobalMaxNumIndices = A.GlobalMaxNumIndices();
01732     for( int i = 0; i < NumExportIDs; ++i )
01733     {
01734       Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
01735       TotalSendSize += Sizes[i];
01736     }
01737   }
01738   catch(...) {
01739     try {
01740       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01741       int MaxNumIndices = A.MaxNumEntries();
01742       A.Comm().MaxAll(&MaxNumIndices, &GlobalMaxNumIndices, 1);
01743       for( int i = 0; i < NumExportIDs; ++i )
01744       {
01745         int NumEntries;
01746         A.NumMyRowEntries( ExportLIDs[i], NumEntries );
01747         Sizes[i] = (NumEntries + 2);
01748         TotalSendSize += Sizes[i];
01749       }
01750     }
01751     catch(...) {
01752       EPETRA_CHK_ERR(-1); // Bad cast
01753     }
01754   }
01755 
01756   CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
01757 
01758   try {
01759     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01760     EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01761             SizeOfPacket, Sizes, VarSizes, Distor));
01762   }
01763   catch(...) {
01764     const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01765     EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01766              SizeOfPacket, Sizes, VarSizes, Distor));
01767   }
01768   return(0);
01769 }
01770 
01771 // private =====================================================================
01772 int Epetra_CrsGraph::PackAndPrepareCrsGraph(const Epetra_CrsGraph& A, 
01773               int NumExportIDs, 
01774               int* ExportLIDs,
01775               int& LenExports, 
01776               char*& Exports, 
01777               int& SizeOfPacket, 
01778                                             int* Sizes,
01779                                             bool& VarSizes,
01780               Epetra_Distributor& Distor)
01781 {
01782   (void)LenExports;
01783   (void)SizeOfPacket;
01784   (void)Sizes;
01785   (void)VarSizes;
01786   (void)Distor;
01787   int i;
01788   int NumIndices;
01789   int* Indices = 0;
01790   int FromRow;
01791   int* intptr;
01792   
01793   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01794   // 1st int: GRID of row where GRID is the global row ID for the source graph
01795   // next int:  NumIndices, Number of indices in row.
01796   // next NumIndices: The actual indices for the row.
01797   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01798   //   sized segments for current communication routines.
01799   int MaxNumIndices = A.MaxNumIndices();
01800   //if( MaxNumIndices ) Indices = new int[MaxNumIndices];
01801 
01802   intptr = (int*) Exports;
01803   for(i = 0; i < NumExportIDs; i++) {
01804     FromRow = A.GRID(ExportLIDs[i]);
01805     *intptr = FromRow;
01806     Indices = intptr + 2;
01807     EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, MaxNumIndices, NumIndices, Indices));
01808     intptr[1] = NumIndices; // Load second slot of segment
01809     intptr += (NumIndices+2); // Point to next segment
01810   }
01811 
01812   //if( Indices ) delete [] Indices;
01813     
01814   return(0);
01815 }
01816 
01817 // private =====================================================================
01818 int Epetra_CrsGraph::PackAndPrepareRowMatrix(const Epetra_RowMatrix& A, 
01819                int NumExportIDs, 
01820                int* ExportLIDs,
01821                int& LenExports, 
01822                char*& Exports, 
01823                int& SizeOfPacket, 
01824                                              int* Sizes,
01825                                              bool& VarSizes,
01826                Epetra_Distributor& Distor)
01827 {
01828   (void)LenExports;
01829   (void)SizeOfPacket;
01830   (void)Sizes;
01831   (void)VarSizes;
01832   (void)Distor;
01833   int i;
01834   int j;
01835   int NumIndices;
01836   int* Indices = 0;
01837   int FromRow;
01838   int* intptr;
01839   Epetra_SerialDenseVector Values;
01840   
01841   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01842   // 1st int: GRID of row where GRID is the global row ID for the source graph
01843   // next int:  NumIndices, Number of indices in row.
01844   // next NumIndices: The actual indices for the row.
01845   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01846   //   sized segments for current communication routines.
01847   int MaxNumIndices = A.MaxNumEntries();
01848   if(MaxNumIndices > 0) {
01849     Values.Size(MaxNumIndices);
01850 //    Indices = new int[MaxNumIndices];
01851   }
01852   const Epetra_Map& RowMap = A.RowMatrixRowMap();
01853   const Epetra_Map& ColMap = A.RowMatrixColMap();
01854 
01855   intptr = (int*) Exports;
01856   for(i = 0; i < NumExportIDs; i++) {
01857     FromRow = RowMap.GID(ExportLIDs[i]);
01858     *intptr = FromRow;
01859     Indices = intptr + 2;
01860     EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], MaxNumIndices, NumIndices, Values.Values(), Indices));
01861     for(j = 0; j < NumIndices; j++) Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01862     intptr[1] = NumIndices; // Load second slot of segment
01863     intptr += (NumIndices+2); // Point to next segment
01864   }
01865 
01866 //  if( Indices ) delete [] Indices;
01867  
01868   return(0);
01869 }
01870 
01871 // private =====================================================================
01872 int Epetra_CrsGraph::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
01873                                       int NumImportIDs, 
01874                                       int* ImportLIDs, 
01875                                       int LenImports, 
01876                                       char* Imports, 
01877                                       int& SizeOfPacket, 
01878                                       Epetra_Distributor& Distor,
01879                                       Epetra_CombineMode CombineMode,
01880                                       const Epetra_OffsetIndex * Indexor) 
01881 {
01882   (void)Source;
01883   (void)LenImports;
01884   (void)SizeOfPacket;
01885   (void)Distor;
01886   (void)CombineMode;
01887   (void)Indexor;
01888   if(NumImportIDs <= 0) 
01889     return(0);
01890 
01891   int NumIndices;
01892   int* Indices;
01893   int ToRow;
01894   int i;
01895   
01896   int* intptr;
01897   // Unpack it...
01898 
01899   // Each segment of Sends will be filled by a packed row of information for each row as follows:
01900   // 1st int: GRID of row where GRID is the global row ID for the source graph
01901   // next int:  NumIndices, Number of indices in row.
01902   // next NumIndices: The actual indices for the row.
01903 
01904   intptr = (int*) Imports;
01905     
01906   for(i = 0; i < NumImportIDs; i++) {
01907     ToRow = GRID(ImportLIDs[i]);
01908     assert((intptr[0])==ToRow); // Sanity check
01909     NumIndices = intptr[1];
01910     Indices = intptr + 2; 
01911     // Insert indices
01912     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices);
01913     if(ierr < 0) 
01914       EPETRA_CHK_ERR(ierr);
01915     intptr += (NumIndices+2); // Point to next segment
01916   }
01917 
01918   //destroy buffers since this operation is usually only done once
01919   if( LenExports_ ) {
01920     delete [] Exports_;
01921     Exports_ = 0;
01922     LenExports_ = 0;
01923   }
01924   if( LenImports_ ) {
01925     delete [] Imports_;
01926     Imports_ = 0;
01927     LenImports_ = 0;
01928   }
01929   
01930   return(0);
01931 }
01932 
01933 // protected ===================================================================
01934 bool Epetra_CrsGraph::GlobalConstantsComputed() const {
01935   int mineComputed = 0;
01936   int allComputed;
01937   if(CrsGraphData_->GlobalConstantsComputed_) 
01938     mineComputed = 1;
01939   RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
01940   // If allComputed comes back zero then at least one PE need global constants recomputed.
01941   return(allComputed==1);
01942 }
01943 
01944 // private =====================================================================
01945 void Epetra_CrsGraph::ComputeIndexState() {
01946   int myIndicesAreLocal = 0;
01947   int myIndicesAreGlobal = 0;
01948   if(CrsGraphData_->IndicesAreLocal_) 
01949     myIndicesAreLocal = 1;
01950   if(CrsGraphData_->IndicesAreGlobal_) 
01951     myIndicesAreGlobal = 1;
01952   int allIndicesAreLocal;
01953   int allIndicesAreGlobal;
01954   RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
01955   RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
01956   CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
01957   CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1);  // If indices are global on one PE should be local on all
01958 }
01959 
01960 //==============================================================================
01961 void Epetra_CrsGraph::Print (ostream& os) const {
01962   int MyPID = RowMap().Comm().MyPID();
01963   int NumProc = RowMap().Comm().NumProc();
01964 
01965   for(int iproc = 0; iproc < NumProc; iproc++) {
01966     if(MyPID == iproc) {
01967       if(MyPID == 0) {
01968   os << "\nNumber of Global Block Rows  = " << NumGlobalBlockRows()      << endl;
01969   os <<   "Number of Global Block Cols  = " << NumGlobalBlockCols()      << endl;
01970   os <<   "Number of Global Block Diags = " << NumGlobalBlockDiagonals() << endl;
01971   os <<   "Number of Global Entries     = " << NumGlobalEntries()        << endl;
01972   os << "\nNumber of Global Rows        = " << NumGlobalRows()           << endl;
01973   os <<   "Number of Global Cols        = " << NumGlobalCols()           << endl;
01974   os <<   "Number of Global Diagonals   = " << NumGlobalDiagonals()      << endl;
01975   os <<   "Number of Global Nonzeros    = " << NumGlobalNonzeros()       << endl;
01976   os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim()         << endl;
01977   os <<   "Global Maximum Block Col Dim = " << GlobalMaxColDim()         << endl;
01978   os <<   "Global Maximum Num Indices   = " << GlobalMaxNumIndices()     << endl;
01979   if(LowerTriangular()) os << " ** Matrix is Lower Triangular **"        << endl;
01980   if(UpperTriangular()) os << " ** Matrix is Upper Triangular **"        << endl;
01981   if(NoDiagonal())      os << " ** Matrix has no diagonal     **"        << endl << endl;
01982       }
01983       os << "\nNumber of My Block Rows  = " << NumMyBlockRows()      << endl;
01984       os <<   "Number of My Block Cols  = " << NumMyBlockCols()      << endl;
01985       os <<   "Number of My Block Diags = " << NumMyBlockDiagonals() << endl;
01986       os <<   "Number of My Entries     = " << NumMyEntries()        << endl;
01987       os << "\nNumber of My Rows        = " << NumMyRows()           << endl;
01988       os <<   "Number of My Cols        = " << NumMyCols()           << endl;
01989       os <<   "Number of My Diagonals   = " << NumMyDiagonals()      << endl;
01990       os <<   "Number of My Nonzeros    = " << NumMyNonzeros()       << endl;
01991       os << "\nMy Maximum Block Row Dim = " << MaxRowDim()           << endl;
01992       os <<   "My Maximum Block Col Dim = " << MaxColDim()           << endl;
01993       os <<   "My Maximum Num Indices   = " << MaxNumIndices()       << endl << endl;
01994 
01995       int NumMyBlockRows1 = NumMyBlockRows();
01996       int MaxNumIndices1 = MaxNumIndices();
01997       Epetra_IntSerialDenseVector Indices1(MaxNumIndices1);
01998       int NumIndices1;
01999       int i;
02000       int j;
02001       
02002       os.width(14);
02003       os <<  "       Row Index "; os << " ";
02004       for(j = 0; j < MaxNumIndices(); j++) {   
02005   os.width(12);
02006   os <<  "Col Index"; os << "      ";
02007       }
02008       os << endl;
02009       for(i = 0; i < NumMyBlockRows1; i++) {
02010   int Row = GRID(i); // Get global row number
02011   ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1.Values());
02012         
02013   os.width(14);
02014   os <<  Row ; os << "    ";  
02015   for(j = 0; j < NumIndices1 ; j++) {   
02016     os.width(12);
02017     os <<  Indices1[j]; os << "    ";
02018   }
02019   os << endl;
02020       }      
02021       os << flush;
02022     }
02023     // Do a few global ops to give I/O a chance to complete
02024     RowMap().Comm().Barrier();
02025     RowMap().Comm().Barrier();
02026     RowMap().Comm().Barrier();
02027   }
02028 }
02029 
02030 //==============================================================================
02031 Epetra_CrsGraph& Epetra_CrsGraph::operator = (const Epetra_CrsGraph& Source) {
02032   if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
02033     return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
02034 
02035   CleanupData();
02036   CrsGraphData_ = Source.CrsGraphData_;
02037   CrsGraphData_->IncrementReferenceCount();
02038 
02039   return(*this);
02040 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines