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

Generated on Thu Sep 18 12:37:57 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1