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 2011 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 global indices to local indices to 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   // This next step constructs the scan sum of the number of indices per row.  Note that the
01339   // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
01340   // careful with this sum operation
01341 
01342   if(CrsGraphData_->IndexOffset_.Values() != CrsGraphData_->NumIndicesPerRow_.Values())
01343     CrsGraphData_->IndexOffset_.MakeViewOf(CrsGraphData_->NumIndicesPerRow_);
01344 
01345   int * numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
01346   int curNumIndices = numIndicesPerRow[0];
01347   numIndicesPerRow[0] = 0;
01348   for (int i=0; i<numMyBlockRows; ++i) {
01349     int nextNumIndices = numIndicesPerRow[i+1];
01350     numIndicesPerRow[i+1] = numIndicesPerRow[i]+curNumIndices;
01351     curNumIndices = nextNumIndices;
01352   }
01353 
01354 // *******************************
01355 // Data NOT contiguous, make it so
01356 // *******************************
01357   if(!Contiguous) { // Must pack indices if not already contiguous
01358 
01359     // Allocate one big integer array for all index values
01360     if (!(StaticProfile())) { // If static profile, All_Indices_ is already allocated, only need to pack data
01361       int errorcode = CrsGraphData_->All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
01362       if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
01363     }
01364       // Pack indices into All_Indices_
01365 
01366     int* all_indices = CrsGraphData_->All_Indices_.Values();
01367     int * indexOffset = CrsGraphData_->IndexOffset_.Values();
01368     int ** indices = CrsGraphData_->Indices_;
01369     
01370     if (!(StaticProfile())) {
01371 #ifdef EPETRA_HAVE_OMP
01372 #pragma omp parallel for default(none) shared(indexOffset,all_indices,indices)
01373 #endif   
01374       for(int i = 0; i < numMyBlockRows; i++) {
01375   int numColIndices = indexOffset[i+1] - indexOffset[i];
01376         int* ColIndices = indices[i];
01377         int *newColIndices = all_indices+indexOffset[i];
01378         for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
01379       }
01380       for(int i = 0; i < numMyBlockRows; i++) {
01381         if (indices[i]!=0) {
01382           CrsGraphData_->SortedEntries_[i].entries_.clear();
01383           indices[i] = 0;
01384         }
01385      }
01386    } // End of non-contiguous non-static section   
01387    else {
01388 
01389      for(int i = 0; i < numMyBlockRows; i++) {
01390        int numColIndices = indexOffset[i+1] - indexOffset[i];
01391        int* ColIndices = indices[i];
01392        int *newColIndices = all_indices+indexOffset[i];
01393        if (ColIndices!=newColIndices) // No need to copy if pointing to same space
01394        for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
01395        indices[i] = 0;
01396      }
01397    } // End of non-Contiguous static section
01398   } // End of non-Contiguous section
01399   else { // Start of Contiguous section
01400     // if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
01401     // Execute the assignment block in parallel using the same pattern as SpMV
01402     // in order to improve page placement
01403     if (numMyBlockRows > 0 && !(StaticProfile())) {
01404       const int numMyNonzeros = NumMyNonzeros();
01405       int errorcode = CrsGraphData_->All_Indices_.Size(numMyNonzeros);
01406       if(errorcode != 0)  throw ReportError("Error with All_Indices_ allocation.", -99);
01407       int* new_all_indices = CrsGraphData_->All_Indices_.Values();
01408       int* old_all_indices = CrsGraphData_->Indices_[0];
01409       int * indexOffset = CrsGraphData_->IndexOffset_.Values();
01410 
01411 #ifdef EPETRA_HAVE_OMP
01412 #pragma omp parallel for default(none) shared(indexOffset,old_all_indices,new_all_indices)
01413 #endif
01414      for(int i = 0; i < numMyBlockRows; i++) {
01415        int numColIndices = indexOffset[i+1] - indexOffset[i];
01416        int *oldColIndices = old_all_indices+indexOffset[i];
01417        int *newColIndices = new_all_indices+indexOffset[i];
01418        for(int j = 0; j < numColIndices; j++) newColIndices[j] = oldColIndices[j];
01419      }
01420 
01421 //#ifdef EPETRA_HAVE_OMP
01422 //#pragma omp parallel for default(none) shared(all_indices_values,indices_values)
01423 //#endif   
01424 //      for(int ii=0; ii<numMyNonzeros; ++ii) {
01425 //        all_indices_values[ii] = indices_values[ii];
01426 //      }
01427     }
01428   }
01429 
01430   // Delete unneeded storage
01431   CrsGraphData_->NumAllocatedIndicesPerRow_.Resize(0);
01432   delete [] CrsGraphData_->Indices_; CrsGraphData_->Indices_=0;
01433   CrsGraphData_->SortedEntries_.clear();
01434 
01435   SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
01436   CrsGraphData_->StorageOptimized_ = true;
01437 
01438   return(0);
01439 }
01440 
01441 //==============================================================================
01442 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01443 {
01444   int j;
01445 
01446   Row = LRID(Row); // Normalize row range
01447 
01448   if(Row < 0 || Row >= NumMyBlockRows()) 
01449     EPETRA_CHK_ERR(-1); // Not in Row range
01450 
01451   NumIndices = NumMyIndices(Row);
01452   if(LenOfIndices < NumIndices) 
01453     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01454 
01455   int * srcIndices = Indices(Row);
01456   if(IndicesAreLocal())  
01457     for(j = 0; j < NumIndices; j++) 
01458       targIndices[j] = GCID(srcIndices[j]);
01459   else 
01460     for(j = 0; j < NumIndices; j++)
01461       targIndices[j] = srcIndices[j];
01462   
01463   return(0);
01464 }
01465 
01466 //==============================================================================
01467 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01468 {
01469   int j;
01470 
01471   if(Row < 0 || Row >= NumMyBlockRows()) 
01472     EPETRA_CHK_ERR(-1); // Not in Row range
01473 
01474   NumIndices = NumMyIndices(Row);
01475   if(LenOfIndices < NumIndices) 
01476     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01477 
01478   if(IndicesAreGlobal()) 
01479     EPETRA_CHK_ERR(-3); // There are no local indices yet
01480 
01481   int * srcIndices = Indices(Row);
01482   for(j = 0; j < NumIndices; j++)
01483     targIndices[j] = srcIndices[j];
01484   
01485   return(0);
01486 }
01487 
01488 //==============================================================================
01489 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const 
01490 {
01491   Row = LRID(Row); // Normalize row range
01492 
01493   if(Row < 0 || Row >= NumMyBlockRows()) 
01494     EPETRA_CHK_ERR(-1); // Not in Row range
01495 
01496   if(IndicesAreLocal()) 
01497     EPETRA_CHK_ERR(-2); // There are no global indices
01498 
01499   NumIndices = NumMyIndices(Row);
01500 
01501   targIndices = Indices(Row);
01502   
01503   return(0);
01504 }
01505 
01506 //==============================================================================
01507 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const 
01508 {
01509   if(Row < 0 || Row >= NumMyBlockRows()) 
01510     EPETRA_CHK_ERR(-1); // Not in Row range
01511 
01512   if(IndicesAreGlobal()) 
01513     EPETRA_CHK_ERR(-2); // There are no local indices
01514 
01515   NumIndices = NumMyIndices(Row);
01516 
01517   targIndices = Indices(Row);
01518   
01519   return(0);
01520 }
01521 
01522 //==============================================================================
01523 int Epetra_CrsGraph::NumGlobalIndices(int Row) const {
01524   Row = LRID(Row);
01525   if(Row != -1) 
01526     return(NumMyIndices(Row));
01527   else 
01528     return(0); // No indices for this row on this processor
01529 }
01530 
01531 //==============================================================================
01532 int Epetra_CrsGraph::NumAllocatedGlobalIndices(int Row) const {
01533   Row = LRID(Row);
01534   if(Row != -1) 
01535     return(NumAllocatedMyIndices(Row));
01536   else 
01537     return(0); // No indices allocated for this row on this processor
01538 }
01539 
01540 //==============================================================================
01541 int Epetra_CrsGraph::ReplaceRowMap(const Epetra_BlockMap& newmap)
01542 {
01543   if (RowMap().PointSameAs(newmap)) {
01544     Epetra_DistObject::Map_ = newmap;
01545     CrsGraphData_->RowMap_ = newmap;
01546     CrsGraphData_->MakeImportExport();
01547     return(0);
01548   }
01549 
01550   return(-1);
01551 }
01552 
01553 //==============================================================================
01554 int Epetra_CrsGraph::ReplaceColMap(const Epetra_BlockMap& newmap)
01555 {
01556   if (!HaveColMap() && !IndicesAreLocal() && !IndicesAreGlobal()) {
01557     CrsGraphData_->ColMap_            = newmap;
01558     CrsGraphData_->NumGlobalBlockCols_= newmap.NumGlobalElements();
01559     CrsGraphData_->NumMyBlockCols_    = newmap.NumMyElements();
01560     CrsGraphData_->MaxColDim_         = newmap.MaxElementSize();
01561     CrsGraphData_->GlobalMaxColDim_   = newmap.MaxElementSize();
01562     CrsGraphData_->NumGlobalCols_     = newmap.NumGlobalPoints();
01563     CrsGraphData_->NumMyCols_         = newmap.NumMyPoints();
01564     CrsGraphData_->HaveColMap_        = true;
01565     return(0);
01566   }
01567 
01568   if(ColMap().PointSameAs(newmap)) {
01569     CrsGraphData_->ColMap_ = newmap;
01570     CrsGraphData_->MakeImportExport();
01571     return(0);
01572   }
01573   
01574   return(-1);
01575 }
01576 
01577 // private =====================================================================
01578 int Epetra_CrsGraph::CheckSizes(const Epetra_SrcDistObject& Source) {
01579   try {
01580     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
01581     if(!A.GlobalConstantsComputed()) 
01582       EPETRA_CHK_ERR(-1); // Must have global constants to proceed
01583   }
01584   catch(...) {
01585     return(0); // No error at this point, object could be a RowMatrix
01586   }
01587   return(0);
01588 }
01589 
01590 // private =====================================================================
01591 int Epetra_CrsGraph::CopyAndPermute(const Epetra_SrcDistObject& Source,
01592             int NumSameIDs, 
01593             int NumPermuteIDs, 
01594             int* PermuteToLIDs,
01595             int* PermuteFromLIDs,
01596                                     const Epetra_OffsetIndex * Indexor)
01597 { 
01598   try {
01599     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01600     EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01601             PermuteFromLIDs,Indexor));
01602   }
01603   catch(...) {
01604     try {
01605       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01606       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01607                PermuteFromLIDs,Indexor));
01608     }
01609     catch(...) {
01610       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
01611     }
01612   }
01613   
01614   return(0);
01615 }
01616 
01617 // private =====================================================================
01618 int Epetra_CrsGraph::CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
01619                int NumSameIDs, 
01620                int NumPermuteIDs, 
01621                int* PermuteToLIDs,
01622                int* PermuteFromLIDs,
01623                                              const Epetra_OffsetIndex * Indexor)
01624 {
01625   (void)Indexor;
01626   int i;
01627   int j;
01628   int NumIndices;
01629   int FromRow;
01630   int ToRow;
01631   int maxNumIndices = A.MaxNumEntries();
01632   Epetra_IntSerialDenseVector indices;
01633   Epetra_SerialDenseVector Values;
01634 
01635   if(maxNumIndices > 0) {
01636     indices.Size(maxNumIndices);
01637     Values.Size(maxNumIndices); // Must extract values even though we discard them
01638   }
01639 
01640   const Epetra_Map& rowMap = A.RowMatrixRowMap();
01641   const Epetra_Map& colMap = A.RowMatrixColMap();
01642   
01643   // Do copy first
01644   for(i = 0; i < NumSameIDs; i++) {
01645     ToRow = rowMap.GID(i);
01646     EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumIndices, NumIndices, Values.Values(), indices.Values()));
01647     for(j = 0; j < NumIndices; j++) 
01648       indices[j] = colMap.GID(indices[j]); // convert to GIDs
01649     // Place into target graph.  
01650     int ierr = InsertGlobalIndices(ToRow, NumIndices, indices.Values());
01651     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01652   }
01653   
01654   // Do local permutation next
01655   for(i = 0; i < NumPermuteIDs; i++) {
01656     FromRow = PermuteFromLIDs[i];
01657     ToRow = GRID(PermuteToLIDs[i]);
01658     EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumIndices, NumIndices, Values.Values(), indices.Values()));
01659     for(j = 0; j < NumIndices; j++) 
01660       indices[j] = colMap.GID(indices[j]); // convert to GIDs
01661     int ierr = InsertGlobalIndices(ToRow, NumIndices, indices.Values()); // Place into target graph.
01662     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01663   }
01664   
01665   return(0);
01666 }
01667 
01668 // private =====================================================================
01669 int Epetra_CrsGraph::CopyAndPermuteCrsGraph(const Epetra_CrsGraph& A,
01670               int NumSameIDs, 
01671               int NumPermuteIDs, 
01672               int* PermuteToLIDs,
01673               int* PermuteFromLIDs,
01674                                             const Epetra_OffsetIndex * Indexor)
01675 {
01676   (void)Indexor;
01677   int i;
01678   int Row;
01679   int NumIndices;
01680   int* indices = 0;
01681   int FromRow, ToRow;
01682   int maxNumIndices = A.MaxNumIndices();
01683   Epetra_IntSerialDenseVector IndicesVector;
01684 
01685   if(maxNumIndices > 0 && A.IndicesAreLocal()) {
01686     IndicesVector.Size(maxNumIndices);
01687     indices = IndicesVector.Values();
01688   }
01689   
01690   // Do copy first
01691   if(NumSameIDs > 0) {
01692     if(A.IndicesAreLocal()) {
01693       for(i = 0; i < NumSameIDs; i++) {
01694         Row = GRID(i);
01695         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumIndices, NumIndices, indices));
01696         // Place into target graph.  
01697         int ierr = InsertGlobalIndices(Row, NumIndices, indices); 
01698         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01699       }
01700     }
01701     else { // A.IndiceAreGlobal()
01702       for(i = 0; i < NumSameIDs; i++) {
01703         Row = GRID(i);
01704         EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, indices));
01705         // Place into target graph.  
01706         int ierr = InsertGlobalIndices(Row, NumIndices, indices); 
01707         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01708       }
01709     } 
01710   }
01711 
01712   // Do local permutation next
01713   if(NumPermuteIDs > 0) {
01714     if(A.IndicesAreLocal()) {
01715       for(i = 0; i < NumPermuteIDs; i++) {
01716         FromRow = A.GRID(PermuteFromLIDs[i]);
01717         ToRow = GRID(PermuteToLIDs[i]);
01718         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
01719         // Place into target graph.
01720         int ierr = InsertGlobalIndices(ToRow, NumIndices, indices); 
01721         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01722       }
01723     }
01724     else { // A.IndiceAreGlobal()
01725       for(i = 0; i < NumPermuteIDs; i++) {
01726         FromRow = A.GRID(PermuteFromLIDs[i]);
01727         ToRow = GRID(PermuteToLIDs[i]);
01728         EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, indices));
01729         // Place into target graph.
01730         int ierr = InsertGlobalIndices(ToRow, NumIndices, indices); 
01731         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01732       }
01733     }
01734   } 
01735   
01736   return(0);
01737 }
01738 
01739 // private =====================================================================
01740 int Epetra_CrsGraph::PackAndPrepare(const Epetra_SrcDistObject& Source, 
01741             int NumExportIDs, 
01742             int* ExportLIDs,
01743             int& LenExports, 
01744             char*& Exports, 
01745             int& SizeOfPacket, 
01746                                     int * Sizes,
01747                                     bool& VarSizes,
01748             Epetra_Distributor& Distor) 
01749 {
01750   int globalMaxNumIndices = 0;
01751   int TotalSendSize = 0;
01752 
01753   VarSizes = true;
01754 
01755   SizeOfPacket = (int)sizeof(int); 
01756 
01757   if(NumExportIDs <= 0) return(0);
01758 
01759   try {
01760     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01761     globalMaxNumIndices = A.GlobalMaxNumIndices();
01762     for( int i = 0; i < NumExportIDs; ++i )
01763     {
01764       Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
01765       TotalSendSize += Sizes[i];
01766     }
01767   }
01768   catch(...) {
01769     try {
01770       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01771       int maxNumIndices = A.MaxNumEntries();
01772       A.Comm().MaxAll(&maxNumIndices, &globalMaxNumIndices, 1);
01773       for( int i = 0; i < NumExportIDs; ++i )
01774       {
01775         int NumEntries;
01776         A.NumMyRowEntries( ExportLIDs[i], NumEntries );
01777         Sizes[i] = (NumEntries + 2);
01778         TotalSendSize += Sizes[i];
01779       }
01780     }
01781     catch(...) {
01782       EPETRA_CHK_ERR(-1); // Bad cast
01783     }
01784   }
01785 
01786   CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
01787 
01788   try {
01789     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01790     EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01791             SizeOfPacket, Sizes, VarSizes, Distor));
01792   }
01793   catch(...) {
01794     const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01795     EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01796              SizeOfPacket, Sizes, VarSizes, Distor));
01797   }
01798   return(0);
01799 }
01800 
01801 // private =====================================================================
01802 int Epetra_CrsGraph::PackAndPrepareCrsGraph(const Epetra_CrsGraph& A, 
01803               int NumExportIDs, 
01804               int* ExportLIDs,
01805               int& LenExports, 
01806               char*& Exports, 
01807               int& SizeOfPacket, 
01808                                             int* Sizes,
01809                                             bool& VarSizes,
01810               Epetra_Distributor& Distor)
01811 {
01812   (void)LenExports;
01813   (void)SizeOfPacket;
01814   (void)Sizes;
01815   (void)VarSizes;
01816   (void)Distor;
01817   int i;
01818   int NumIndices;
01819   int* indices = 0;
01820   int FromRow;
01821   int* intptr;
01822   
01823   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01824   // 1st int: GRID of row where GRID is the global row ID for the source graph
01825   // next int:  NumIndices, Number of indices in row.
01826   // next NumIndices: The actual indices for the row.
01827   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01828   //   sized segments for current communication routines.
01829   int maxNumIndices = A.MaxNumIndices();
01830   //if( maxNumIndices ) indices = new int[maxNumIndices];
01831 
01832   intptr = (int*) Exports;
01833   for(i = 0; i < NumExportIDs; i++) {
01834     FromRow = A.GRID(ExportLIDs[i]);
01835     *intptr = FromRow;
01836     indices = intptr + 2;
01837     EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
01838     intptr[1] = NumIndices; // Load second slot of segment
01839     intptr += (NumIndices+2); // Point to next segment
01840   }
01841 
01842   //if( indices ) delete [] indices;
01843     
01844   return(0);
01845 }
01846 
01847 // private =====================================================================
01848 int Epetra_CrsGraph::PackAndPrepareRowMatrix(const Epetra_RowMatrix& A, 
01849                int NumExportIDs, 
01850                int* ExportLIDs,
01851                int& LenExports, 
01852                char*& Exports, 
01853                int& SizeOfPacket, 
01854                                              int* Sizes,
01855                                              bool& VarSizes,
01856                Epetra_Distributor& Distor)
01857 {
01858   (void)LenExports;
01859   (void)SizeOfPacket;
01860   (void)Sizes;
01861   (void)VarSizes;
01862   (void)Distor;
01863   int i;
01864   int j;
01865   int NumIndices;
01866   int* indices = 0;
01867   int FromRow;
01868   int* intptr;
01869   Epetra_SerialDenseVector Values;
01870   
01871   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01872   // 1st int: GRID of row where GRID is the global row ID for the source graph
01873   // next int:  NumIndices, Number of indices in row.
01874   // next NumIndices: The actual indices for the row.
01875   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01876   //   sized segments for current communication routines.
01877   int maxNumIndices = A.MaxNumEntries();
01878   if(maxNumIndices > 0) {
01879     Values.Size(maxNumIndices);
01880 //    indices = new int[maxNumIndices];
01881   }
01882   const Epetra_Map& rowMap = A.RowMatrixRowMap();
01883   const Epetra_Map& colMap = A.RowMatrixColMap();
01884 
01885   intptr = (int*) Exports;
01886   for(i = 0; i < NumExportIDs; i++) {
01887     FromRow = rowMap.GID(ExportLIDs[i]);
01888     *intptr = FromRow;
01889     indices = intptr + 2;
01890     EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), indices));
01891     for(j = 0; j < NumIndices; j++) indices[j] = colMap.GID(indices[j]); // convert to GIDs
01892     intptr[1] = NumIndices; // Load second slot of segment
01893     intptr += (NumIndices+2); // Point to next segment
01894   }
01895 
01896 //  if( indices ) delete [] indices;
01897  
01898   return(0);
01899 }
01900 
01901 // private =====================================================================
01902 int Epetra_CrsGraph::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
01903                                       int NumImportIDs, 
01904                                       int* ImportLIDs, 
01905                                       int LenImports, 
01906                                       char* Imports, 
01907                                       int& SizeOfPacket, 
01908                                       Epetra_Distributor& Distor,
01909                                       Epetra_CombineMode CombineMode,
01910                                       const Epetra_OffsetIndex * Indexor) 
01911 {
01912   (void)Source;
01913   (void)LenImports;
01914   (void)SizeOfPacket;
01915   (void)Distor;
01916   (void)CombineMode;
01917   (void)Indexor;
01918   if(NumImportIDs <= 0) 
01919     return(0);
01920 
01921   int NumIndices;
01922   int* indices;
01923   int ToRow;
01924   int i;
01925   
01926   int* intptr;
01927   // Unpack it...
01928 
01929   // Each segment of Sends will be filled by a packed row of information for each row as follows:
01930   // 1st int: GRID of row where GRID is the global row ID for the source graph
01931   // next int:  NumIndices, Number of indices in row.
01932   // next NumIndices: The actual indices for the row.
01933 
01934   intptr = (int*) Imports;
01935     
01936   for(i = 0; i < NumImportIDs; i++) {
01937     ToRow = GRID(ImportLIDs[i]);
01938     assert((intptr[0])==ToRow); // Sanity check
01939     NumIndices = intptr[1];
01940     indices = intptr + 2; 
01941     // Insert indices
01942     int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
01943     if(ierr < 0) 
01944       EPETRA_CHK_ERR(ierr);
01945     intptr += (NumIndices+2); // Point to next segment
01946   }
01947 
01948   //destroy buffers since this operation is usually only done once
01949   if( LenExports_ ) {
01950     delete [] Exports_;
01951     Exports_ = 0;
01952     LenExports_ = 0;
01953   }
01954   if( LenImports_ ) {
01955     delete [] Imports_;
01956     Imports_ = 0;
01957     LenImports_ = 0;
01958   }
01959   
01960   return(0);
01961 }
01962 
01963 // protected ===================================================================
01964 bool Epetra_CrsGraph::GlobalConstantsComputed() const {
01965   int mineComputed = 0;
01966   int allComputed;
01967   if(CrsGraphData_->GlobalConstantsComputed_) 
01968     mineComputed = 1;
01969   RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
01970   // If allComputed comes back zero then at least one PE need global constants recomputed.
01971   return(allComputed==1);
01972 }
01973 
01974 // private =====================================================================
01975 void Epetra_CrsGraph::ComputeIndexState() {
01976   int myIndicesAreLocal = 0;
01977   int myIndicesAreGlobal = 0;
01978   if(CrsGraphData_->IndicesAreLocal_) 
01979     myIndicesAreLocal = 1;
01980   if(CrsGraphData_->IndicesAreGlobal_) 
01981     myIndicesAreGlobal = 1;
01982   int allIndicesAreLocal;
01983   int allIndicesAreGlobal;
01984   RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
01985   RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
01986   CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
01987   CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1);  // If indices are global on one PE should be local on all
01988 }
01989 
01990 //==============================================================================
01991 void Epetra_CrsGraph::Print (ostream& os) const {
01992   int MyPID = RowMap().Comm().MyPID();
01993   int NumProc = RowMap().Comm().NumProc();
01994 
01995   for(int iproc = 0; iproc < NumProc; iproc++) {
01996     if(MyPID == iproc) {
01997       if(MyPID == 0) {
01998   os << "\nNumber of Global Block Rows  = " << NumGlobalBlockRows()      << endl;
01999   os <<   "Number of Global Block Cols  = " << NumGlobalBlockCols()      << endl;
02000   os <<   "Number of Global Block Diags = " << NumGlobalBlockDiagonals() << endl;
02001   os <<   "Number of Global Entries     = " << NumGlobalEntries()        << endl;
02002   os << "\nNumber of Global Rows        = " << NumGlobalRows()           << endl;
02003   os <<   "Number of Global Cols        = " << NumGlobalCols()           << endl;
02004   os <<   "Number of Global Diagonals   = " << NumGlobalDiagonals()      << endl;
02005   os <<   "Number of Global Nonzeros    = " << NumGlobalNonzeros()       << endl;
02006   os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim()         << endl;
02007   os <<   "Global Maximum Block Col Dim = " << GlobalMaxColDim()         << endl;
02008   os <<   "Global Maximum Num Indices   = " << GlobalMaxNumIndices()     << endl;
02009   if(LowerTriangular()) os << " ** Matrix is Lower Triangular **"        << endl;
02010   if(UpperTriangular()) os << " ** Matrix is Upper Triangular **"        << endl;
02011   if(NoDiagonal())      os << " ** Matrix has no diagonal     **"        << endl << endl;
02012       }
02013       os << "\nNumber of My Block Rows  = " << NumMyBlockRows()      << endl;
02014       os <<   "Number of My Block Cols  = " << NumMyBlockCols()      << endl;
02015       os <<   "Number of My Block Diags = " << NumMyBlockDiagonals() << endl;
02016       os <<   "Number of My Entries     = " << NumMyEntries()        << endl;
02017       os << "\nNumber of My Rows        = " << NumMyRows()           << endl;
02018       os <<   "Number of My Cols        = " << NumMyCols()           << endl;
02019       os <<   "Number of My Diagonals   = " << NumMyDiagonals()      << endl;
02020       os <<   "Number of My Nonzeros    = " << NumMyNonzeros()       << endl;
02021       os << "\nMy Maximum Block Row Dim = " << MaxRowDim()           << endl;
02022       os <<   "My Maximum Block Col Dim = " << MaxColDim()           << endl;
02023       os <<   "My Maximum Num Indices   = " << MaxNumIndices()       << endl << endl;
02024 
02025       int NumMyBlockRows1 = NumMyBlockRows();
02026       int MaxNumIndices1 = MaxNumIndices();
02027       Epetra_IntSerialDenseVector Indices1(MaxNumIndices1);
02028       int NumIndices1;
02029       int i;
02030       int j;
02031       
02032       os.width(14);
02033       os <<  "       Row Index "; os << " ";
02034       for(j = 0; j < MaxNumIndices(); j++) {   
02035   os.width(12);
02036   os <<  "Col Index"; os << "      ";
02037       }
02038       os << endl;
02039       for(i = 0; i < NumMyBlockRows1; i++) {
02040   int Row = GRID(i); // Get global row number
02041   ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1.Values());
02042         
02043   os.width(14);
02044   os <<  Row ; os << "    ";  
02045   for(j = 0; j < NumIndices1 ; j++) {   
02046     os.width(12);
02047     os <<  Indices1[j]; os << "    ";
02048   }
02049   os << endl;
02050       }      
02051       os << flush;
02052     }
02053     // Do a few global ops to give I/O a chance to complete
02054     RowMap().Comm().Barrier();
02055     RowMap().Comm().Barrier();
02056     RowMap().Comm().Barrier();
02057   }
02058 }
02059 
02060 //==============================================================================
02061 Epetra_CrsGraph& Epetra_CrsGraph::operator = (const Epetra_CrsGraph& Source) {
02062   if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
02063     return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
02064 
02065   CleanupData();
02066   CrsGraphData_ = Source.CrsGraphData_;
02067   CrsGraphData_->IncrementReferenceCount();
02068 
02069   return(*this);
02070 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines