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 the max of 100
00967   // and the number of block rows.
00968   const int numMyBlockRows = NumMyBlockRows();
00969   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
00970   //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << endl;
00971   Epetra_HashTable RemoteGIDs(hashsize); 
00972   Epetra_HashTable RemoteGIDList(hashsize);
00973 
00974   int NumLocalColGIDs = 0;
00975   int NumRemoteColGIDs = 0;
00976   for(i = 0; i < numMyBlockRows; i++) {
00977     const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
00978     int* ColIndices = CrsGraphData_->Indices_[i];
00979     for(j = 0; j < NumIndices; j++) {
00980       int GID = ColIndices[j];
00981       // Check if GID matches a row GID
00982       int LID = DomainMap.LID(GID);
00983       if(LID != -1) {
00984   bool alreadyFound = LocalGIDs[LID];
00985   if (!alreadyFound) {
00986           LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
00987           NumLocalColGIDs++;
00988   }
00989       }
00990       else {
00991   if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
00992     RemoteGIDs.Add(GID, NumRemoteColGIDs);
00993     RemoteGIDList.Add(NumRemoteColGIDs++, GID);
00994   }
00995       }
00996     }
00997   }
00998 
00999   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=DomainMap and quit
01000   if (DomainMap.Comm().NumProc()==1) { 
01001     
01002     if (NumRemoteColGIDs!=0) {
01003       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
01004     }
01005     if (NumLocalColGIDs==numDomainElements) {
01006       CrsGraphData_->ColMap_ = DomainMap;
01007       CrsGraphData_->HaveColMap_ = true;
01008       if (LocalGIDs!=0) delete [] LocalGIDs; 
01009       return(0); 
01010     }
01011   }
01012       
01013   // Now build integer array containing column GIDs
01014   // Build back end, containing remote GIDs, first
01015   int NumMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
01016   Epetra_IntSerialDenseVector ColIndices;
01017   if(NumMyBlockCols > 0) 
01018     ColIndices.Size(NumMyBlockCols);
01019 
01020   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01021 
01022   for(i = 0; i < NumRemoteColGIDs; i++) 
01023     RemoteColIndices[i] = RemoteGIDList.Get(i); 
01024 
01025   int NLists = 1;
01026   Epetra_IntSerialDenseVector PIDList;
01027   Epetra_IntSerialDenseVector SizeList;
01028   int* RemoteSizeList = 0;
01029   bool DoSizes = !DomainMap.ConstantElementSize(); // If not constant element size, then we must exchange
01030       
01031   if(NumRemoteColGIDs > 0) 
01032     PIDList.Size(NumRemoteColGIDs);
01033 
01034   if(DoSizes) {
01035     if(NumMyBlockCols > 0) 
01036       SizeList.Size(NumMyBlockCols);
01037     RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
01038     NLists++;
01039   }
01040   EPETRA_CHK_ERR(DomainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01041       
01042   // Sort External column indices so that all columns coming from a given remote processor are contiguous
01043 
01044   Epetra_Util Util;
01045   int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
01046   SortLists[0] = RemoteColIndices;
01047   SortLists[1] = RemoteSizeList;
01048   Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists);
01049   if (CrsGraphData_->SortGhostsAssociatedWithEachProcessor_) {
01050     // Sort external column indices so that columns from a given remote processor are not only contiguous
01051     // but also in ascending order. NOTE: I don't know if the number of externals associated
01052     // with a given remote processor is known at this point ... so I count them here.
01053 
01054     NLists--;
01055     int StartCurrent, StartNext;
01056     StartCurrent = 0; StartNext = 1;
01057     while ( StartNext < NumRemoteColGIDs ) {
01058       if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
01059       else {
01060         if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01061         Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists);
01062         StartCurrent = StartNext; StartNext++;
01063       }
01064     }
01065     if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01066     Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists);
01067   }
01068 
01069   // Now fill front end. Two cases:
01070   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
01071   //     can simply read the domain GIDs into the front part of ColIndices, otherwise 
01072   // (2) We step through the GIDs of the DomainMap, checking to see if each domain GID is a column GID.
01073   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
01074 
01075   if(NumLocalColGIDs == DomainMap.NumMyElements()) {
01076     DomainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first NumMyBlockCols elements column GID list
01077     if(DoSizes) 
01078       DomainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
01079   }
01080   else {
01081     int NumMyElements = DomainMap.NumMyElements();
01082     int* MyGlobalElements = DomainMap.MyGlobalElements();
01083     int* ElementSizeList = 0;
01084     if(DoSizes) 
01085       ElementSizeList = DomainMap.ElementSizeList();
01086     int NumLocalAgain = 0;
01087     for(i = 0; i < NumMyElements; i++) {
01088       if(LocalGIDs[i]) {
01089   if(DoSizes) 
01090     SizeList[NumLocalAgain] = ElementSizeList[i];
01091   ColIndices[NumLocalAgain++] = MyGlobalElements[i];
01092       }
01093     }
01094     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
01095   }
01096 
01097   // Done with this array
01098   if (LocalGIDs!=0) delete [] LocalGIDs; 
01099 
01100 
01101   // Make Column map with same element sizes as Domain map
01102 
01103   if(DomainMap.MaxElementSize() == 1) { // Simple map
01104     Epetra_Map temp(-1, NumMyBlockCols, ColIndices.Values(), DomainMap.IndexBase(), DomainMap.Comm());
01105     CrsGraphData_->ColMap_ = temp;
01106   }
01107   else if(DomainMap.ConstantElementSize()) { // Constant Block size map
01108     Epetra_BlockMap temp(-1, NumMyBlockCols, ColIndices.Values(), DomainMap.MaxElementSize(),DomainMap.IndexBase(), DomainMap.Comm());
01109     CrsGraphData_->ColMap_ = temp;
01110   }
01111   else { // Most general case where block size is variable.
01112     Epetra_BlockMap temp(-1, NumMyBlockCols, ColIndices.Values(), SizeList.Values(), DomainMap.IndexBase(), DomainMap.Comm());
01113     CrsGraphData_->ColMap_ = temp;
01114   }
01115   CrsGraphData_->HaveColMap_ = true;
01116 
01117   return(0);
01118 }
01119 
01120 // protected ===================================================================
01121 int Epetra_CrsGraph::MakeIndicesLocal(const Epetra_BlockMap& DomainMap, const Epetra_BlockMap& RangeMap) {
01122   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01123   if(IndicesAreLocal() && IndicesAreGlobal()) 
01124     EPETRA_CHK_ERR(-1); // Return error: Indices must not be both local and global
01125 
01126   MakeColMap(DomainMap, RangeMap); // If user has not prescribed column map, create one from indices
01127   const Epetra_BlockMap& colmap = ColMap();
01128 
01129   // Store number of local columns
01130   CrsGraphData_->NumMyCols_ = ColMap().NumMyPoints();
01131   CrsGraphData_->NumMyBlockCols_ = ColMap().NumMyElements();
01132   // Transform indices to local index space
01133 
01134   const int numMyBlockRows = NumMyBlockRows();
01135 
01136   if(IndicesAreGlobal()) {
01137     for(int i = 0; i < numMyBlockRows; i++) {
01138       const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01139       int* ColIndices = CrsGraphData_->Indices_[i];
01140       for(int j = 0; j < NumIndices; j++) {
01141   int GID = ColIndices[j];
01142   int LID = colmap.LID(GID);
01143   if(LID != -1) 
01144     ColIndices[j] = LID;
01145   else 
01146     throw ReportError("Internal error in FillComplete ",-1); 
01147       }
01148     }
01149   }
01150   
01151   SetIndicesAreLocal(true);
01152   SetIndicesAreGlobal(false);
01153   
01154   if(CrsGraphData_->ReferenceCount() > 1)
01155     return(1); // return 1 if data is shared
01156   else
01157     return(0);
01158 }
01159 
01160 //==============================================================================
01161 int Epetra_CrsGraph::OptimizeStorage() {
01162   int i;
01163   int j;
01164   int NumIndices;
01165   const int numMyBlockRows = NumMyBlockRows();
01166 
01167   if(StorageOptimized()) 
01168     return(0); // Have we been here before?
01169   if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
01170 
01171   bool Contiguous = true; // Assume contiguous is true
01172   for(i = 1; i < numMyBlockRows; i++) {
01173     NumIndices = CrsGraphData_->NumIndicesPerRow_[i-1];
01174     int NumAllocateIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[i-1];
01175 
01176     // Check if NumIndices is same as NumAllocatedIndices and 
01177     // check if end of beginning of current row starts immediately after end of previous row.
01178     if((NumIndices != NumAllocateIndices) || 
01179        (CrsGraphData_->Indices_[i] != CrsGraphData_->Indices_[i-1] + NumIndices)) {
01180       Contiguous = false;
01181       break;
01182     }
01183   }
01184 
01185   // NOTE:  At the end of the above loop set, there is a possibility that NumIndices and NumAllocatedIndices
01186   //        for the last row could be different, but I don't think it matters.
01187 
01188 
01189   if((CrsGraphData_->CV_ == View) && !Contiguous) 
01190     return(3);  // This is user data, it's not contiguous and we can't make it so.
01191 
01192   if(CrsGraphData_->IndexOffset_ .Values() != CrsGraphData_->NumIndicesPerRow_.Values())
01193     CrsGraphData_->IndexOffset_.MakeViewOf(CrsGraphData_->NumIndicesPerRow_);
01194 
01195   // This next step constructs the scan sum of the number of indices per row.  Note that the
01196   // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
01197   // careful with this sum operation
01198   int * indices = CrsGraphData_->NumIndicesPerRow_.Values();
01199   int curNumIndices = indices[0];
01200   indices[0] = 0;
01201   for (i=0; i<numMyBlockRows; ++i) {
01202     int nextNumIndices = indices[i+1];
01203     indices[i+1] = indices[i]+curNumIndices;
01204     curNumIndices = nextNumIndices;
01205   }
01206 
01207   if(!Contiguous) { // Must pack indices if not already contiguous
01208 
01209     // Allocate one big integer array for all index values
01210     if (!(CrsGraphData_->StaticProfile_)) { // If static profile, All_Indices_ is already allocated, only need to pack data
01211       int errorcode = CrsGraphData_->All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
01212       if(errorcode != 0) 
01213   throw ReportError("Error with All_Indices_ allocation.", -99);
01214     }
01215     
01216     // Pack indices into All_Indices_
01217     
01218     int* tmp = CrsGraphData_->All_Indices_.Values();
01219     int * IndexOffset = CrsGraphData_->IndexOffset_ .Values();
01220     for(i = 0; i < numMyBlockRows; i++) {
01221       NumIndices = IndexOffset[i+1] - IndexOffset[i];
01222       int*& ColIndices = CrsGraphData_->Indices_[i];
01223       if (tmp!=ColIndices) { // No need to copy if pointing to same space
01224   for(j = 0; j < NumIndices; j++) 
01225     tmp[j] = ColIndices[j];
01226       }
01227       if (!(CrsGraphData_->StaticProfile_) && ColIndices!=0) delete [] ColIndices;
01228       CrsGraphData_->Indices_[i] = 0;
01229       tmp += NumIndices;  // tmp points to the offset in All_Indices_ where Indices_[i] starts.
01230     }
01231   } // End of !Contiguous section
01232   else {
01233     //if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
01234     if (numMyBlockRows > 0 && !(CrsGraphData_->StaticProfile_)) {
01235       int errorcode = CrsGraphData_->All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
01236       if(errorcode != 0) 
01237   throw ReportError("Error with All_Indices_ allocation.", -99);
01238       int* all_indices_values = CrsGraphData_->All_Indices_.Values();
01239       int* indices_values = CrsGraphData_->Indices_[0];
01240       for(int ii=0; ii<CrsGraphData_->NumMyNonzeros_; ++ii) {
01241         all_indices_values[ii] = indices_values[ii];
01242       }
01243     }
01244   }
01245 
01246   // Delete unneeded storage
01247   CrsGraphData_->NumAllocatedIndicesPerRow_.Resize(0);
01248   delete [] CrsGraphData_->Indices_; CrsGraphData_->Indices_=0;
01249 
01250   SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
01251   CrsGraphData_->StorageOptimized_ = true;
01252 
01253   return(0);
01254 }
01255 
01256 //==============================================================================
01257 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01258 {
01259   int j;
01260 
01261   Row = LRID(Row); // Normalize row range
01262 
01263   if(Row < 0 || Row >= NumMyBlockRows()) 
01264     EPETRA_CHK_ERR(-1); // Not in Row range
01265 
01266   NumIndices = NumMyIndices(Row);
01267   if(LenOfIndices < NumIndices) 
01268     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01269 
01270   int * srcIndices = Indices(Row);
01271   if(IndicesAreLocal())  
01272     for(j = 0; j < NumIndices; j++) 
01273       targIndices[j] = GCID(srcIndices[j]);
01274   else 
01275     for(j = 0; j < NumIndices; j++)
01276       targIndices[j] = srcIndices[j];
01277   
01278   return(0);
01279 }
01280 
01281 //==============================================================================
01282 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const 
01283 {
01284   int j;
01285 
01286   if(Row < 0 || Row >= NumMyBlockRows()) 
01287     EPETRA_CHK_ERR(-1); // Not in Row range
01288 
01289   NumIndices = NumMyIndices(Row);
01290   if(LenOfIndices < NumIndices) 
01291     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
01292 
01293   if(IndicesAreGlobal()) 
01294     EPETRA_CHK_ERR(-3); // There are no local indices yet
01295 
01296   int * srcIndices = Indices(Row);
01297   for(j = 0; j < NumIndices; j++)
01298     targIndices[j] = srcIndices[j];
01299   
01300   return(0);
01301 }
01302 
01303 //==============================================================================
01304 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const 
01305 {
01306   Row = LRID(Row); // Normalize row range
01307 
01308   if(Row < 0 || Row >= NumMyBlockRows()) 
01309     EPETRA_CHK_ERR(-1); // Not in Row range
01310 
01311   if(IndicesAreLocal()) 
01312     EPETRA_CHK_ERR(-2); // There are no global indices
01313 
01314   NumIndices = NumMyIndices(Row);
01315 
01316   targIndices = Indices(Row);
01317   
01318   return(0);
01319 }
01320 
01321 //==============================================================================
01322 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const 
01323 {
01324   if(Row < 0 || Row >= NumMyBlockRows()) 
01325     EPETRA_CHK_ERR(-1); // Not in Row range
01326 
01327   if(IndicesAreGlobal()) 
01328     EPETRA_CHK_ERR(-2); // There are no local indices
01329 
01330   NumIndices = NumMyIndices(Row);
01331 
01332   targIndices = Indices(Row);
01333   
01334   return(0);
01335 }
01336 
01337 //==============================================================================
01338 int Epetra_CrsGraph::NumGlobalIndices(int Row) const {
01339   Row = LRID(Row);
01340   if(Row != -1) 
01341     return(NumMyIndices(Row));
01342   else 
01343     return(0); // No indices for this row on this processor
01344 }
01345 
01346 //==============================================================================
01347 int Epetra_CrsGraph::NumAllocatedGlobalIndices(int Row) const {
01348   Row = LRID(Row);
01349   if(Row != -1) 
01350     return(NumAllocatedMyIndices(Row));
01351   else 
01352     return(0); // No indices allocated for this row on this processor
01353 }
01354 
01355 //==============================================================================
01356 int Epetra_CrsGraph::ReplaceRowMap(const Epetra_BlockMap& newmap)
01357 {
01358   if (RowMap().PointSameAs(newmap)) {
01359     Epetra_DistObject::Map_ = newmap;
01360     CrsGraphData_->RowMap_ = newmap;
01361     CrsGraphData_->MakeImportExport();
01362     return(0);
01363   }
01364 
01365   return(-1);
01366 }
01367 
01368 //==============================================================================
01369 int Epetra_CrsGraph::ReplaceColMap(const Epetra_BlockMap& newmap)
01370 {
01371   if (ColMap().PointSameAs(newmap)) {
01372     CrsGraphData_->ColMap_ = newmap;
01373     CrsGraphData_->MakeImportExport();
01374     return(0);
01375   }
01376 
01377   return(-1);
01378 }
01379 
01380 // private =====================================================================
01381 int Epetra_CrsGraph::CheckSizes(const Epetra_SrcDistObject& Source) {
01382   try {
01383     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
01384     if(!A.GlobalConstantsComputed()) 
01385       EPETRA_CHK_ERR(-1); // Must have global constants to proceed
01386   }
01387   catch(...) {
01388     return(0); // No error at this point, object could be a RowMatrix
01389   }
01390   return(0);
01391 }
01392 
01393 // private =====================================================================
01394 int Epetra_CrsGraph::CopyAndPermute(const Epetra_SrcDistObject& Source,
01395             int NumSameIDs, 
01396             int NumPermuteIDs, 
01397             int* PermuteToLIDs,
01398             int* PermuteFromLIDs,
01399                                     const Epetra_OffsetIndex * Indexor)
01400 { 
01401   try {
01402     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01403     EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01404             PermuteFromLIDs,Indexor));
01405   }
01406   catch(...) {
01407     try {
01408       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01409       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01410                PermuteFromLIDs,Indexor));
01411     }
01412     catch(...) {
01413       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
01414     }
01415   }
01416   
01417   return(0);
01418 }
01419 
01420 // private =====================================================================
01421 int Epetra_CrsGraph::CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
01422                int NumSameIDs, 
01423                int NumPermuteIDs, 
01424                int* PermuteToLIDs,
01425                int* PermuteFromLIDs,
01426                                              const Epetra_OffsetIndex * Indexor)
01427 {
01428   (void)Indexor;
01429   int i;
01430   int j;
01431   int NumIndices;
01432   int FromRow;
01433   int ToRow;
01434   int MaxNumIndices = A.MaxNumEntries();
01435   Epetra_IntSerialDenseVector Indices;
01436   Epetra_SerialDenseVector Values;
01437 
01438   if(MaxNumIndices > 0) {
01439     Indices.Size(MaxNumIndices);
01440     Values.Size(MaxNumIndices); // Must extract values even though we discard them
01441   }
01442 
01443   const Epetra_Map& RowMap = A.RowMatrixRowMap();
01444   const Epetra_Map& ColMap = A.RowMatrixColMap();
01445   
01446   // Do copy first
01447   for(i = 0; i < NumSameIDs; i++) {
01448     ToRow = RowMap.GID(i);
01449     EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, MaxNumIndices, NumIndices, Values.Values(), Indices.Values()));
01450     for(j = 0; j < NumIndices; j++) 
01451       Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01452     // Place into target graph.  
01453     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices.Values());
01454     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01455   }
01456   
01457   // Do local permutation next
01458   for(i = 0; i < NumPermuteIDs; i++) {
01459     FromRow = PermuteFromLIDs[i];
01460     ToRow = GRID(PermuteToLIDs[i]);
01461     EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, MaxNumIndices, NumIndices, Values.Values(), Indices.Values()));
01462     for(j = 0; j < NumIndices; j++) 
01463       Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01464     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices.Values()); // Place into target graph.
01465     if(ierr < 0) EPETRA_CHK_ERR(ierr);
01466   }
01467   
01468   return(0);
01469 }
01470 
01471 // private =====================================================================
01472 int Epetra_CrsGraph::CopyAndPermuteCrsGraph(const Epetra_CrsGraph& A,
01473               int NumSameIDs, 
01474               int NumPermuteIDs, 
01475               int* PermuteToLIDs,
01476               int* PermuteFromLIDs,
01477                                             const Epetra_OffsetIndex * Indexor)
01478 {
01479   (void)Indexor;
01480   int i;
01481   int Row;
01482   int NumIndices;
01483   int* Indices = 0;
01484   int FromRow, ToRow;
01485   int MaxNumIndices = A.MaxNumIndices();
01486   Epetra_IntSerialDenseVector IndicesVector;
01487 
01488   if(MaxNumIndices > 0 && A.IndicesAreLocal()) {
01489     IndicesVector.Size(MaxNumIndices);
01490     Indices = IndicesVector.Values();
01491   }
01492   
01493   // Do copy first
01494   if(NumSameIDs > 0) {
01495     if(A.IndicesAreLocal()) {
01496       for(i = 0; i < NumSameIDs; i++) {
01497         Row = GRID(i);
01498         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumIndices, Indices));
01499         // Place into target graph.  
01500         int ierr = InsertGlobalIndices(Row, NumIndices, Indices); 
01501         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01502       }
01503     }
01504     else { // A.IndiceAreGlobal()
01505       for(i = 0; i < NumSameIDs; i++) {
01506         Row = GRID(i);
01507         EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, Indices));
01508         // Place into target graph.  
01509         int ierr = InsertGlobalIndices(Row, NumIndices, Indices); 
01510         if(ierr < 0) EPETRA_CHK_ERR(ierr); 
01511       }
01512     } 
01513   }
01514 
01515   // Do local permutation next
01516   if(NumPermuteIDs > 0) {
01517     if(A.IndicesAreLocal()) {
01518       for(i = 0; i < NumPermuteIDs; i++) {
01519         FromRow = A.GRID(PermuteFromLIDs[i]);
01520         ToRow = GRID(PermuteToLIDs[i]);
01521         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, MaxNumIndices, NumIndices, Indices));
01522         // Place into target graph.
01523         int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices); 
01524         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01525       }
01526     }
01527     else { // A.IndiceAreGlobal()
01528       for(i = 0; i < NumPermuteIDs; i++) {
01529         FromRow = A.GRID(PermuteFromLIDs[i]);
01530         ToRow = GRID(PermuteToLIDs[i]);
01531         EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, Indices));
01532         // Place into target graph.
01533         int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices); 
01534         if (ierr < 0) EPETRA_CHK_ERR(ierr); 
01535       }
01536     }
01537   } 
01538   
01539   return(0);
01540 }
01541 
01542 // private =====================================================================
01543 int Epetra_CrsGraph::PackAndPrepare(const Epetra_SrcDistObject& Source, 
01544             int NumExportIDs, 
01545             int* ExportLIDs,
01546             int& LenExports, 
01547             char*& Exports, 
01548             int& SizeOfPacket, 
01549                                     int * Sizes,
01550                                     bool& VarSizes,
01551             Epetra_Distributor& Distor) 
01552 {
01553   int GlobalMaxNumIndices = 0;
01554   int TotalSendSize = 0;
01555 
01556   VarSizes = true;
01557 
01558   SizeOfPacket = (int)sizeof(int); 
01559 
01560   if(NumExportIDs <= 0) return(0);
01561 
01562   try {
01563     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01564     GlobalMaxNumIndices = A.GlobalMaxNumIndices();
01565     for( int i = 0; i < NumExportIDs; ++i )
01566     {
01567       Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
01568       TotalSendSize += Sizes[i];
01569     }
01570   }
01571   catch(...) {
01572     try {
01573       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01574       int MaxNumIndices = A.MaxNumEntries();
01575       A.Comm().MaxAll(&MaxNumIndices, &GlobalMaxNumIndices, 1);
01576       for( int i = 0; i < NumExportIDs; ++i )
01577       {
01578         int NumEntries;
01579         A.NumMyRowEntries( ExportLIDs[i], NumEntries );
01580         Sizes[i] = (NumEntries + 2);
01581         TotalSendSize += Sizes[i];
01582       }
01583     }
01584     catch(...) {
01585       EPETRA_CHK_ERR(-1); // Bad cast
01586     }
01587   }
01588 
01589   CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
01590 
01591   try {
01592     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
01593     EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01594             SizeOfPacket, Sizes, VarSizes, Distor));
01595   }
01596   catch(...) {
01597     const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
01598     EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
01599              SizeOfPacket, Sizes, VarSizes, Distor));
01600   }
01601   return(0);
01602 }
01603 
01604 // private =====================================================================
01605 int Epetra_CrsGraph::PackAndPrepareCrsGraph(const Epetra_CrsGraph& A, 
01606               int NumExportIDs, 
01607               int* ExportLIDs,
01608               int& LenExports, 
01609               char*& Exports, 
01610               int& SizeOfPacket, 
01611                                             int* Sizes,
01612                                             bool& VarSizes,
01613               Epetra_Distributor& Distor)
01614 {
01615   (void)LenExports;
01616   (void)SizeOfPacket;
01617   (void)Sizes;
01618   (void)VarSizes;
01619   (void)Distor;
01620   int i;
01621   int NumIndices;
01622   int* Indices = 0;
01623   int FromRow;
01624   int* intptr;
01625   
01626   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01627   // 1st int: GRID of row where GRID is the global row ID for the source graph
01628   // next int:  NumIndices, Number of indices in row.
01629   // next NumIndices: The actual indices for the row.
01630   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01631   //   sized segments for current communication routines.
01632   int MaxNumIndices = A.MaxNumIndices();
01633   //if( MaxNumIndices ) Indices = new int[MaxNumIndices];
01634 
01635   intptr = (int*) Exports;
01636   for(i = 0; i < NumExportIDs; i++) {
01637     FromRow = A.GRID(ExportLIDs[i]);
01638     *intptr = FromRow;
01639     Indices = intptr + 2;
01640     EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, MaxNumIndices, NumIndices, Indices));
01641     intptr[1] = NumIndices; // Load second slot of segment
01642     intptr += (NumIndices+2); // Point to next segment
01643   }
01644 
01645   //if( Indices ) delete [] Indices;
01646     
01647   return(0);
01648 }
01649 
01650 // private =====================================================================
01651 int Epetra_CrsGraph::PackAndPrepareRowMatrix(const Epetra_RowMatrix& A, 
01652                int NumExportIDs, 
01653                int* ExportLIDs,
01654                int& LenExports, 
01655                char*& Exports, 
01656                int& SizeOfPacket, 
01657                                              int* Sizes,
01658                                              bool& VarSizes,
01659                Epetra_Distributor& Distor)
01660 {
01661   (void)LenExports;
01662   (void)SizeOfPacket;
01663   (void)Sizes;
01664   (void)VarSizes;
01665   (void)Distor;
01666   int i;
01667   int j;
01668   int NumIndices;
01669   int* Indices = 0;
01670   int FromRow;
01671   int* intptr;
01672   Epetra_SerialDenseVector Values;
01673   
01674   // Each segment of Exports will be filled by a packed row of information for each row as follows:
01675   // 1st int: GRID of row where GRID is the global row ID for the source graph
01676   // next int:  NumIndices, Number of indices in row.
01677   // next NumIndices: The actual indices for the row.
01678   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
01679   //   sized segments for current communication routines.
01680   int MaxNumIndices = A.MaxNumEntries();
01681   if(MaxNumIndices > 0) {
01682     Values.Size(MaxNumIndices);
01683 //    Indices = new int[MaxNumIndices];
01684   }
01685   const Epetra_Map& RowMap = A.RowMatrixRowMap();
01686   const Epetra_Map& ColMap = A.RowMatrixColMap();
01687 
01688   intptr = (int*) Exports;
01689   for(i = 0; i < NumExportIDs; i++) {
01690     FromRow = RowMap.GID(ExportLIDs[i]);
01691     *intptr = FromRow;
01692     Indices = intptr + 2;
01693     EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], MaxNumIndices, NumIndices, Values.Values(), Indices));
01694     for(j = 0; j < NumIndices; j++) Indices[j] = ColMap.GID(Indices[j]); // convert to GIDs
01695     intptr[1] = NumIndices; // Load second slot of segment
01696     intptr += (NumIndices+2); // Point to next segment
01697   }
01698 
01699 //  if( Indices ) delete [] Indices;
01700  
01701   return(0);
01702 }
01703 
01704 // private =====================================================================
01705 int Epetra_CrsGraph::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
01706                                       int NumImportIDs, 
01707                                       int* ImportLIDs, 
01708                                       int LenImports, 
01709                                       char* Imports, 
01710                                       int& SizeOfPacket, 
01711                                       Epetra_Distributor& Distor,
01712                                       Epetra_CombineMode CombineMode,
01713                                       const Epetra_OffsetIndex * Indexor) 
01714 {
01715   (void)Source;
01716   (void)LenImports;
01717   (void)SizeOfPacket;
01718   (void)Distor;
01719   (void)CombineMode;
01720   (void)Indexor;
01721   if(NumImportIDs <= 0) 
01722     return(0);
01723 
01724   int NumIndices;
01725   int* Indices;
01726   int ToRow;
01727   int i;
01728   
01729   int* intptr;
01730   // Unpack it...
01731 
01732   // Each segment of Sends will be filled by a packed row of information for each row as follows:
01733   // 1st int: GRID of row where GRID is the global row ID for the source graph
01734   // next int:  NumIndices, Number of indices in row.
01735   // next NumIndices: The actual indices for the row.
01736 
01737   intptr = (int*) Imports;
01738     
01739   for(i = 0; i < NumImportIDs; i++) {
01740     ToRow = GRID(ImportLIDs[i]);
01741     assert((intptr[0])==ToRow); // Sanity check
01742     NumIndices = intptr[1];
01743     Indices = intptr + 2; 
01744     // Insert indices
01745     int ierr = InsertGlobalIndices(ToRow, NumIndices, Indices);
01746     if(ierr < 0) 
01747       EPETRA_CHK_ERR(ierr);
01748     intptr += (NumIndices+2); // Point to next segment
01749   }
01750 
01751   //destroy buffers since this operation is usually only done once
01752   if( LenExports_ ) {
01753     delete [] Exports_;
01754     Exports_ = 0;
01755     LenExports_ = 0;
01756   }
01757   if( LenImports_ ) {
01758     delete [] Imports_;
01759     Imports_ = 0;
01760     LenImports_ = 0;
01761   }
01762   
01763   return(0);
01764 }
01765 
01766 // protected ===================================================================
01767 bool Epetra_CrsGraph::GlobalConstantsComputed() const {
01768   int mineComputed = 0;
01769   int allComputed;
01770   if(CrsGraphData_->GlobalConstantsComputed_) 
01771     mineComputed = 1;
01772   RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
01773   // If allComputed comes back zero then at least one PE need global constants recomputed.
01774   return(allComputed==1);
01775 }
01776 
01777 // private =====================================================================
01778 void Epetra_CrsGraph::ComputeIndexState() {
01779   int myIndicesAreLocal = 0;
01780   int myIndicesAreGlobal = 0;
01781   if(CrsGraphData_->IndicesAreLocal_) 
01782     myIndicesAreLocal = 1;
01783   if(CrsGraphData_->IndicesAreGlobal_) 
01784     myIndicesAreGlobal = 1;
01785   int allIndicesAreLocal;
01786   int allIndicesAreGlobal;
01787   RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
01788   RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
01789   CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
01790   CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1);  // If indices are global on one PE should be local on all
01791 }
01792 
01793 //==============================================================================
01794 void Epetra_CrsGraph::Print (ostream& os) const {
01795   int MyPID = RowMap().Comm().MyPID();
01796   int NumProc = RowMap().Comm().NumProc();
01797 
01798   for(int iproc = 0; iproc < NumProc; iproc++) {
01799     if(MyPID == iproc) {
01800       if(MyPID == 0) {
01801   os << "\nNumber of Global Block Rows  = " << NumGlobalBlockRows()      << endl;
01802   os <<   "Number of Global Block Cols  = " << NumGlobalBlockCols()      << endl;
01803   os <<   "Number of Global Block Diags = " << NumGlobalBlockDiagonals() << endl;
01804   os <<   "Number of Global Entries     = " << NumGlobalEntries()        << endl;
01805   os << "\nNumber of Global Rows        = " << NumGlobalRows()           << endl;
01806   os <<   "Number of Global Cols        = " << NumGlobalCols()           << endl;
01807   os <<   "Number of Global Diagonals   = " << NumGlobalDiagonals()      << endl;
01808   os <<   "Number of Global Nonzeros    = " << NumGlobalNonzeros()       << endl;
01809   os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim()         << endl;
01810   os <<   "Global Maximum Block Col Dim = " << GlobalMaxColDim()         << endl;
01811   os <<   "Global Maximum Num Indices   = " << GlobalMaxNumIndices()     << endl;
01812   if(LowerTriangular()) os << " ** Matrix is Lower Triangular **"        << endl;
01813   if(UpperTriangular()) os << " ** Matrix is Upper Triangular **"        << endl;
01814   if(NoDiagonal())      os << " ** Matrix has no diagonal     **"        << endl << endl;
01815       }
01816       os << "\nNumber of My Block Rows  = " << NumMyBlockRows()      << endl;
01817       os <<   "Number of My Block Cols  = " << NumMyBlockCols()      << endl;
01818       os <<   "Number of My Block Diags = " << NumMyBlockDiagonals() << endl;
01819       os <<   "Number of My Entries     = " << NumMyEntries()        << endl;
01820       os << "\nNumber of My Rows        = " << NumMyRows()           << endl;
01821       os <<   "Number of My Cols        = " << NumMyCols()           << endl;
01822       os <<   "Number of My Diagonals   = " << NumMyDiagonals()      << endl;
01823       os <<   "Number of My Nonzeros    = " << NumMyNonzeros()       << endl;
01824       os << "\nMy Maximum Block Row Dim = " << MaxRowDim()           << endl;
01825       os <<   "My Maximum Block Col Dim = " << MaxColDim()           << endl;
01826       os <<   "My Maximum Num Indices   = " << MaxNumIndices()       << endl << endl;
01827 
01828       int NumMyBlockRows1 = NumMyBlockRows();
01829       int MaxNumIndices1 = MaxNumIndices();
01830       Epetra_IntSerialDenseVector Indices1(MaxNumIndices1);
01831       int NumIndices1;
01832       int i;
01833       int j;
01834       
01835       os.width(14);
01836       os <<  "       Row Index "; os << " ";
01837       for(j = 0; j < MaxNumIndices(); j++) {   
01838   os.width(12);
01839   os <<  "Col Index"; os << "      ";
01840       }
01841       os << endl;
01842       for(i = 0; i < NumMyBlockRows1; i++) {
01843   int Row = GRID(i); // Get global row number
01844   ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1.Values());
01845         
01846   os.width(14);
01847   os <<  Row ; os << "    ";  
01848   for(j = 0; j < NumIndices1 ; j++) {   
01849     os.width(12);
01850     os <<  Indices1[j]; os << "    ";
01851   }
01852   os << endl;
01853       }      
01854       os << flush;
01855     }
01856     // Do a few global ops to give I/O a chance to complete
01857     RowMap().Comm().Barrier();
01858     RowMap().Comm().Barrier();
01859     RowMap().Comm().Barrier();
01860   }
01861 }
01862 
01863 //==============================================================================
01864 Epetra_CrsGraph& Epetra_CrsGraph::operator = (const Epetra_CrsGraph& Source) {
01865   if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
01866     return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
01867 
01868   CleanupData();
01869   CrsGraphData_ = Source.CrsGraphData_;
01870   CrsGraphData_->IncrementReferenceCount();
01871 
01872   return(*this);
01873 }

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7