Epetra Package Browser (Single Doxygen Collection) Development
Epetra_CrsGraph.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //               Epetra: Linear Algebra Services Package
00006 //                 Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_ConfigDefs.h"
00045 #include "Epetra_CrsGraph.h"
00046 #include "Epetra_Import.h"
00047 #include "Epetra_Export.h"
00048 #include "Epetra_Distributor.h"
00049 #include "Epetra_Util.h"
00050 #include "Epetra_Comm.h"
00051 #include "Epetra_HashTable.h"
00052 #include "Epetra_BlockMap.h"
00053 #include "Epetra_Map.h"
00054 #include "Epetra_RowMatrix.h"
00055 #include "Epetra_IntSerialDenseVector.h"
00056 #include "Epetra_SerialComm.h"
00057 
00058 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00059 #include "Epetra_LongLongSerialDenseVector.h"
00060 #endif
00061 
00062 #include "Epetra_SerialDenseVector.h"
00063 #include "Epetra_OffsetIndex.h"
00064 
00065 //==============================================================================
00066 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV,
00067         const Epetra_BlockMap& rowMap,
00068         const int* numIndicesPerRow, bool staticProfile)
00069   : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
00070     CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, staticProfile))
00071 {
00072   Allocate(numIndicesPerRow, 1, staticProfile);
00073 }
00074 
00075 //==============================================================================
00076 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV,
00077         const Epetra_BlockMap& rowMap,
00078         int numIndicesPerRow, bool staticProfile)
00079   : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
00080     CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, staticProfile))
00081 {
00082   Allocate(&numIndicesPerRow, 0, staticProfile);
00083 }
00084 
00085 //==============================================================================
00086 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV,
00087          const Epetra_BlockMap& rowMap,
00088          const Epetra_BlockMap& colMap,
00089          const int* numIndicesPerRow, bool staticProfile)
00090   : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
00091     CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, colMap, staticProfile))
00092 {
00093   if(!rowMap.GlobalIndicesTypeMatch(colMap))
00094      throw ReportError("Epetra_CrsGraph::Epetra_CrsGraph: cannot be called with different indices types for rowMap and colMap", -1);
00095 
00096   Allocate(numIndicesPerRow, 1, staticProfile);
00097 }
00098 
00099 //==============================================================================
00100 Epetra_CrsGraph::Epetra_CrsGraph(Epetra_DataAccess CV,
00101          const Epetra_BlockMap& rowMap,
00102          const Epetra_BlockMap& colMap,
00103          int numIndicesPerRow, bool staticProfile)
00104   : Epetra_DistObject(rowMap, "Epetra::CrsGraph"),
00105     CrsGraphData_(new Epetra_CrsGraphData(CV, rowMap, colMap, staticProfile))
00106 {
00107   if(!rowMap.GlobalIndicesTypeMatch(colMap))
00108      throw ReportError("Epetra_CrsGraph::Epetra_CrsGraph: cannot be called with different indices types for rowMap and colMap", -1);
00109 
00110   Allocate(&numIndicesPerRow, 0, staticProfile);
00111 }
00112 
00113 //==============================================================================
00114 Epetra_CrsGraph::Epetra_CrsGraph(const Epetra_CrsGraph& Graph)
00115   : Epetra_DistObject(Graph),
00116     CrsGraphData_(Graph.CrsGraphData_)
00117 {
00118   CrsGraphData_->IncrementReferenceCount();
00119 }
00120 
00121 // private =====================================================================
00122 template<typename int_type>
00123 int Epetra_CrsGraph::TAllocate(const int* numIndicesPerRow, int Inc, bool staticProfile) {
00124   int i;
00125   const int numMyBlockRows = CrsGraphData_->NumMyBlockRows_;
00126 
00127   // Portions specific to ISDVs
00128   // Note that NumIndicesPerRow_ will be 1 element longer than needed.
00129   // This is because, if OptimizeStorage() is called, the storage associated with
00130   // NumIndicesPerRow_ will be reused implicitly for the IndexOffset_ array.
00131   // Although a bit fragile, for users who care about efficient memory allocation,
00132   // the time and memory fragmentation are important: Mike Heroux Feb 2005.
00133   int errorcode = CrsGraphData_->NumIndicesPerRow_.Size(numMyBlockRows+1);
00134   if(errorcode != 0)
00135     throw ReportError("Error with NumIndicesPerRow_ allocation.", -99);
00136 
00137   errorcode = CrsGraphData_->NumAllocatedIndicesPerRow_.Size(numMyBlockRows);
00138   if(errorcode != 0)
00139     throw ReportError("Error with NumAllocatedIndicesPerRow_ allocation.", -99);
00140 
00141   int nnz = 0;
00142   if(CrsGraphData_->CV_ == Copy) {
00143     if (numIndicesPerRow != 0) {
00144       for(i = 0; i < numMyBlockRows; i++) {
00145   int nnzr = numIndicesPerRow[i*Inc];
00146   nnz += nnzr;
00147       }
00148     }
00149   }
00150   CrsGraphData_->NumMyEntries_ = nnz; // Define this now since we have it and might want to use it
00151   //***
00152 
00153   CrsGraphData_->MaxNumIndices_ = 0;
00154 
00155   Epetra_CrsGraphData::IndexData<int_type>& Data = CrsGraphData_->Data<int_type>();
00156 
00157   // Allocate and initialize entries if we are copying data
00158   if(CrsGraphData_->CV_ == Copy) {
00159     if (staticProfile) Data.All_Indices_.Size(nnz);
00160     int_type * all_indices = Data.All_Indices_.Values(); // First address of contiguous buffer
00161     for(i = 0; i < numMyBlockRows; i++) {
00162       const int NumIndices = numIndicesPerRow==0 ? 0 :numIndicesPerRow[i*Inc];
00163       const int_type indexBaseMinusOne = (int_type) IndexBase64() - 1;
00164 
00165       if(NumIndices > 0) {
00166   if (staticProfile) {
00167     Data.Indices_[i] = all_indices;
00168     all_indices += NumIndices;
00169     int_type* ColIndices = Data.Indices_[i];
00170     for(int j = 0; j < NumIndices; j++)
00171       ColIndices[j] = indexBaseMinusOne; // Fill column indices with out-of-range values
00172   }
00173   else {
00174     // reserve memory in the STL vector, and then resize it to zero
00175     // again in order to signal the program that no data is in there
00176     // yet.
00177     Data.SortedEntries_[i].entries_.resize(NumIndices,
00178                  indexBaseMinusOne);
00179     Data.Indices_[i] = NumIndices > 0 ? &Data.SortedEntries_[i].entries_[0]: NULL;
00180     Data.SortedEntries_[i].entries_.resize(0);
00181   }
00182       }
00183       else {
00184   Data.Indices_[i] = 0;
00185       }
00186 
00187       CrsGraphData_->NumAllocatedIndicesPerRow_[i] = NumIndices;
00188     }
00189     if (staticProfile) assert(Data.All_Indices_.Values()+nnz==all_indices); // Sanity check
00190   }
00191   else { // CV_ == View
00192     for(i = 0; i < numMyBlockRows; i++) {
00193       Data.Indices_[i] = 0;
00194     }
00195   }
00196 
00197   SetAllocated(true);
00198 
00199   return(0);
00200 }
00201 
00202 int Epetra_CrsGraph::Allocate(const int* numIndicesPerRow, int Inc, bool staticProfile)
00203 {
00204   if(RowMap().GlobalIndicesInt()) {
00205     return TAllocate<int>(numIndicesPerRow, Inc, staticProfile);
00206   }
00207 
00208   if(RowMap().GlobalIndicesLongLong()) {
00209 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00210     TAllocate<int>(numIndicesPerRow, Inc, staticProfile);
00211     TAllocate<long long>(numIndicesPerRow, Inc, staticProfile);
00212     return 0;
00213 #else
00214     throw ReportError("Epetra_CrsGraph::Allocate: ERROR, GlobalIndicesLongLong but no API for it.",-1);
00215 #endif
00216   }
00217 
00218   throw ReportError("Epetra_CrsGraph::Allocate: Internal error.", -1);
00219 }
00220 
00221 
00222 // private =====================================================================
00223 /*
00224 int Epetra_CrsGraph::ReAllocate() {
00225   // Reallocate storage that was deleted in OptimizeStorage
00226 
00227   // Since NIPR is in Copy mode, NAIPR will become a Copy as well
00228   CrsGraphData_->NumAllocatedIndicesPerRow_ = CrsGraphData_->NumIndicesPerRow_;
00229 
00230   CrsGraphData_->StorageOptimized_ = false;
00231 
00232   return(0);
00233 }
00234 */
00235 //==============================================================================
00236 Epetra_CrsGraph::~Epetra_CrsGraph()
00237 {
00238   CleanupData();
00239 }
00240 
00241 // private =====================================================================
00242 void Epetra_CrsGraph::CleanupData() {
00243   if(CrsGraphData_ != 0) {
00244     CrsGraphData_->DecrementReferenceCount();
00245     if(CrsGraphData_->ReferenceCount() == 0) {
00246       delete CrsGraphData_;
00247       CrsGraphData_ = 0;
00248     }
00249   }
00250 }
00251 
00252 //==============================================================================
00253 template<typename int_type>
00254 int Epetra_CrsGraph::InsertGlobalIndices(int_type Row, int NumIndices, int_type* indices) {
00255   if(IndicesAreLocal())
00256     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00257   if(IndicesAreContiguous())
00258     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00259   SetIndicesAreGlobal(true);
00260   int locRow = LRID(Row); // Find local row number for this global row index
00261 
00262   EPETRA_CHK_ERR(InsertIndicesIntoSorted(locRow, NumIndices, indices));
00263 
00264   if(CrsGraphData_->ReferenceCount() > 1)
00265     return(1);
00266   else
00267     return(0);
00268 }
00269 
00270 //==============================================================================
00271 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00272 int Epetra_CrsGraph::InsertGlobalIndices(int Row, int NumIndices, int* indices) {
00273 
00274   if(RowMap().GlobalIndicesInt())
00275     return InsertGlobalIndices<int>(Row, NumIndices, indices);
00276   else
00277     throw ReportError("Epetra_CrsGraph::InsertGlobalIndices int version called for a graph that is not int.", -1);
00278 }
00279 #endif
00280 //==============================================================================
00281 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00282 int Epetra_CrsGraph::InsertGlobalIndices(long long Row, int NumIndices, long long* indices) {
00283 
00284   if(RowMap().GlobalIndicesLongLong())
00285     return InsertGlobalIndices<long long>(Row, NumIndices, indices);
00286   else
00287     throw ReportError("Epetra_CrsGraph::InsertGlobalIndices long long version called for a graph that is not long long.", -1);
00288 }
00289 #endif
00290 //==============================================================================
00291 int Epetra_CrsGraph::InsertMyIndices(int Row, int NumIndices, int* indices) {
00292 
00293   if(IndicesAreGlobal()) {
00294     EPETRA_CHK_ERR(-2); // Cannot insert local indices into a global graph
00295   }
00296   if(IndicesAreContiguous())
00297     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00298 
00299   if (CrsGraphData_->HaveColMap_) {
00300     SetIndicesAreLocal(true);
00301   }
00302   else {
00303      if (!IndicesAreLocal()) {
00304        EPETRA_CHK_ERR(-4);
00305      }
00306   }
00307 
00308 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00309   EPETRA_CHK_ERR(InsertIndicesIntoSorted(Row, NumIndices, indices));
00310 #else
00311   throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
00312 #endif
00313 
00314   if(CrsGraphData_->ReferenceCount() > 1)
00315     return(1);
00316   else
00317     return(0);
00318 }
00319 
00320 // protected ===================================================================
00321 template<typename int_type>
00322 int Epetra_CrsGraph::InsertIndices(int Row,
00323            int NumIndices,
00324            int_type* UserIndices)
00325 {
00326   if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
00327 
00328   SetSorted(false); // No longer in sorted state.
00329   CrsGraphData_->NoRedundancies_ = false; // Redundancies possible.
00330   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00331 
00332   int j;
00333   int ierr = 0;
00334 
00335   if(Row < 0 || Row >= NumMyBlockRows())
00336     EPETRA_CHK_ERR(-2); // Not in Row range
00337 
00338   int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
00339   int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00340 
00341   Epetra_CrsGraphData::IndexData<int_type>& Data = CrsGraphData_->Data<int_type>();
00342 
00343   if(CrsGraphData_->CV_ == View) {
00344     if(Data.Indices_[Row] != 0)
00345       ierr = 2; // This row has been defined already.  Issue warning.
00346     Data.Indices_[Row] = UserIndices;
00347     current_numAllocIndices = NumIndices;
00348     current_numIndices = NumIndices;
00349   }
00350   else {
00351     // if HaveColMap_ is true, UserIndices is copied into a new array,
00352     // and then modified. The UserIndices pointer is updated to point to this
00353     // new array. If HaveColMap_ is false, nothing is done. This way,
00354     // the same UserIndices pointer can be used later on regardless of whether
00355     // changes were made.
00356     if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
00357       if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
00358         delete [] Data.TempColIndices_;
00359         Data.TempColIndices_ = new int_type[NumIndices];
00360         CrsGraphData_->NumTempColIndices_ = NumIndices;
00361       }
00362       int_type * tempIndices = Data.TempColIndices_;
00363       int loc = 0;
00364       if(IndicesAreLocal()) {
00365         for(j = 0; j < NumIndices; ++j)
00366           if(CrsGraphData_->ColMap_.MyLID(static_cast<int>(UserIndices[j])))
00367             tempIndices[loc++] = UserIndices[j];
00368       }
00369       else {
00370         for(j = 0; j < NumIndices; ++j) {
00371           const int Index = CrsGraphData_->ColMap_.LID(UserIndices[j]);
00372           if (Index > -1)
00373             tempIndices[loc++] = UserIndices[j];
00374         }
00375 
00376       }
00377       if(loc != NumIndices)
00378         ierr = 2; //Some columns excluded
00379       NumIndices = loc;
00380       UserIndices = tempIndices;
00381     }
00382 
00383     int start = current_numIndices;
00384     int stop = start + NumIndices;
00385     if (CrsGraphData_->StaticProfile_) {
00386       if(stop > current_numAllocIndices)
00387         EPETRA_CHK_ERR(-2); // Cannot reallocate storage if graph created using StaticProfile
00388     }
00389     else {
00390       if (current_numAllocIndices > 0 && stop > current_numAllocIndices)
00391         ierr = 3;
00392       Data.SortedEntries_[Row].entries_.resize(stop, IndexBase64() - 1);
00393       Data.Indices_[Row] = stop>0 ? &Data.SortedEntries_[Row].entries_[0] : NULL;
00394 
00395       current_numAllocIndices =  (int) Data.SortedEntries_[Row].entries_.capacity();
00396     }
00397 
00398     current_numIndices = stop;
00399     int_type* RowIndices = Data.Indices_[Row]+start;
00400     for(j = 0; j < NumIndices; j++) {
00401       RowIndices[j] = UserIndices[j];
00402     }
00403   }
00404 
00405   if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
00406     CrsGraphData_->MaxNumIndices_ = current_numIndices;
00407   }
00408   EPETRA_CHK_ERR(ierr);
00409 
00410 
00411   if(CrsGraphData_->ReferenceCount() > 1)
00412     return(1); // return 1 if data is shared
00413   else
00414     return(0);
00415 }
00416 
00417 // =========================================================================
00418 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00419 int Epetra_CrsGraph::InsertIndices(int Row,
00420            int NumIndices,
00421            int* UserIndices)
00422 {
00423   if(RowMap().GlobalIndicesTypeValid())
00424     return InsertIndices<int>(Row, NumIndices, UserIndices);
00425   else
00426     throw ReportError("Epetra_CrsGraph::InsertIndices global index type unknown.", -1);
00427 }
00428 #endif
00429 // =========================================================================
00430 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00431 int Epetra_CrsGraph::InsertIndices(int Row,
00432            int NumIndices,
00433            long long* UserIndices)
00434 {
00435   if(RowMap().GlobalIndicesLongLong())
00436     return InsertIndices<long long>(Row, NumIndices, UserIndices);
00437   else
00438     throw ReportError("Epetra_CrsGraph::InsertIndices long long version called for a graph that is not long long.", -1);
00439 }
00440 #endif
00441 
00442 // =========================================================================
00443 template<typename int_type>
00444 int Epetra_CrsGraph::InsertIndicesIntoSorted(int Row,
00445               int NumIndices,
00446               int_type* UserIndices)
00447 {
00448   // This function is only valid for COPY mode with non-static profile and
00449   // sorted entries. Otherwise, go to the other function.
00450   if (!CrsGraphData_->NoRedundancies_ || !CrsGraphData_->Sorted_  ||
00451       CrsGraphData_->StaticProfile_ || CrsGraphData_->CV_ == View )
00452     return InsertIndices(Row, NumIndices, UserIndices);
00453 
00454   if (StorageOptimized()) EPETRA_CHK_ERR(-1); // Cannot insert into an optimized graph
00455 
00456   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00457 
00458   int ierr = 0;
00459 
00460   if(Row < 0 || Row >= NumMyBlockRows())
00461     EPETRA_CHK_ERR(-2); // Not in Row range
00462 
00463   int& current_numAllocIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[Row];
00464   int& current_numIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00465 
00466   Epetra_CrsGraphData::IndexData<int_type>& Data = CrsGraphData_->Data<int_type>();
00467 
00468   // if HaveColMap_ is true, UserIndices filters out excluded indices,
00469   // and then modified. The UserIndices pointer is updated to point to this
00470   // new array. If HaveColMap_ is false, nothing is done. This way,
00471   // the same UserIndices pointer can be used later on regardless of whether
00472   // changes were made.
00473   if(CrsGraphData_->HaveColMap_) { //only insert indices in col map if defined
00474     if (CrsGraphData_->NumTempColIndices_ < NumIndices) {
00475       delete [] Data.TempColIndices_;
00476       Data.TempColIndices_ = new int_type[NumIndices];
00477       CrsGraphData_->NumTempColIndices_ = NumIndices;
00478     }
00479     int_type * tempIndices = Data.TempColIndices_;
00480     int loc = 0;
00481     if(IndicesAreLocal()) {
00482       for(int j = 0; j < NumIndices; ++j)
00483         if(CrsGraphData_->ColMap_.MyLID(static_cast<int>(UserIndices[j])))
00484           tempIndices[loc++] = UserIndices[j];
00485     }
00486     else {
00487       for(int j = 0; j < NumIndices; ++j) {
00488         if (CrsGraphData_->ColMap_.MyGID(UserIndices[j])) {
00489           tempIndices[loc++] = UserIndices[j];
00490         }
00491       }
00492     }
00493     if(loc != NumIndices)
00494       ierr = 2; //Some columns excluded
00495     NumIndices = loc;
00496     UserIndices = tempIndices;
00497   }
00498 
00499   // for non-static profile, directly insert into a list that we always
00500   // keep sorted.
00501   Data.SortedEntries_[Row].AddEntries(NumIndices, UserIndices);
00502   current_numIndices = (int) Data.SortedEntries_[Row].entries_.size();
00503   current_numAllocIndices = (int) Data.SortedEntries_[Row].entries_.capacity();
00504   // reset the pointer to the respective data
00505   Data.Indices_[Row] = current_numIndices > 0 ? &Data.SortedEntries_[Row].entries_[0] : NULL;
00506 
00507   if (CrsGraphData_->MaxNumIndices_ < current_numIndices) {
00508     CrsGraphData_->MaxNumIndices_ = current_numIndices;
00509   }
00510   EPETRA_CHK_ERR(ierr);
00511 
00512   if(CrsGraphData_->ReferenceCount() > 1)
00513     return(1); // return 1 if data is shared
00514   else
00515     return(0);
00516 }
00517 
00518 //==============================================================================
00519 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00520 int Epetra_CrsGraph::InsertIndicesIntoSorted(int Row,
00521               int NumIndices,
00522               int* UserIndices)
00523 {
00524   if(RowMap().GlobalIndicesTypeValid())
00525     return InsertIndicesIntoSorted<int>(Row, NumIndices, UserIndices);
00526   else
00527     throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted global index type unknown.", -1);
00528 }
00529 #endif
00530 //==============================================================================
00531 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00532 int Epetra_CrsGraph::InsertIndicesIntoSorted(int Row,
00533               int NumIndices,
00534               long long* UserIndices)
00535 {
00536   if(RowMap().GlobalIndicesLongLong())
00537     return InsertIndicesIntoSorted<long long>(Row, NumIndices, UserIndices);
00538   else
00539     throw ReportError("Epetra_CrsGraph::InsertIndicesIntoSorted long long version called for a graph that is not long long.", -1);
00540 }
00541 #endif
00542 //==============================================================================
00543 template<typename int_type>
00544 int Epetra_CrsGraph::RemoveGlobalIndices(int_type Row, int NumIndices, int_type* indices) {
00545   int j;
00546   int k;
00547   int ierr = 0;
00548   int Loc;
00549 
00550   if(IndicesAreContiguous() || StorageOptimized())
00551     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00552 
00553   if(IndicesAreLocal())
00554     EPETRA_CHK_ERR(-2); // Cannot remove global indices from a filled graph
00555 
00556   if(CrsGraphData_->CV_ == View)
00557     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00558 
00559   int locRow = LRID(Row); // Normalize row range
00560 
00561   if(locRow < 0 || locRow >= NumMyBlockRows())
00562     EPETRA_CHK_ERR(-1); // Not in Row range
00563 
00564   int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[locRow];
00565 
00566   Epetra_CrsGraphData::IndexData<int_type>& Data = CrsGraphData_->Data<int_type>();
00567 
00568   for(j = 0; j < NumIndices; j++) {
00569     int_type Index = indices[j];
00570     if(FindGlobalIndexLoc(locRow,Index,j,Loc)) {
00571       for(k = Loc+1; k < NumCurrentIndices; k++)
00572         Data.Indices_[locRow][k-1] = Data.Indices_[locRow][k];
00573       NumCurrentIndices--;
00574       CrsGraphData_->NumIndicesPerRow_[locRow]--;
00575       if (!CrsGraphData_->StaticProfile_)
00576         Data.SortedEntries_[locRow].entries_.pop_back();
00577       else
00578         Data.Indices_[locRow][NumCurrentIndices-1] = IndexBase64() - 1;
00579     }
00580   }
00581   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00582 
00583   EPETRA_CHK_ERR(ierr);
00584 
00585   if(CrsGraphData_->ReferenceCount() > 1)
00586     return(1);
00587   else
00588     return(0);
00589 }
00590 
00591 //==============================================================================
00592 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00593 int Epetra_CrsGraph::RemoveGlobalIndices(int Row, int NumIndices, int* indices)
00594 {
00595   if(RowMap().GlobalIndicesInt())
00596     return RemoveGlobalIndices<int>(Row, NumIndices, indices);
00597   else
00598     throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices int version called for a graph that is not int.", -1);
00599 }
00600 #endif
00601 //==============================================================================
00602 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00603 int Epetra_CrsGraph::RemoveGlobalIndices(long long Row, int NumIndices, long long* indices)
00604 {
00605   if(RowMap().GlobalIndicesLongLong())
00606     return RemoveGlobalIndices<long long>(Row, NumIndices, indices);
00607   else
00608     throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices long long version called for a graph that is not long long.", -1);
00609 }
00610 #endif
00611 //==============================================================================
00612 int Epetra_CrsGraph::RemoveMyIndices(int Row, int NumIndices, int* indices) {
00613 
00614   if(IndicesAreContiguous() || StorageOptimized())
00615     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00616 
00617   if(IndicesAreGlobal())
00618     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00619 
00620   int j;
00621   int k;
00622   int ierr = 0;
00623   int Loc;
00624 
00625   if(CrsGraphData_->CV_ == View)
00626     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00627 
00628   if(Row < 0 || Row >= NumMyBlockRows())
00629     EPETRA_CHK_ERR(-1); // Not in Row range
00630 
00631   int NumCurrentIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00632 
00633   Epetra_CrsGraphData::IndexData<int>& Data = CrsGraphData_->Data<int>();
00634 
00635   for(j = 0; j < NumIndices; j++) {
00636     int Index = indices[j];
00637     if(FindMyIndexLoc(Row,Index,j,Loc)) {
00638       for(k = Loc + 1; k < NumCurrentIndices; k++)
00639         Data.Indices_[Row][k-1] = Data.Indices_[Row][k];
00640       NumCurrentIndices--;
00641       CrsGraphData_->NumIndicesPerRow_[Row]--;
00642       if (!CrsGraphData_->StaticProfile_)
00643         Data.SortedEntries_[Row].entries_.pop_back();
00644       else
00645         Data.Indices_[Row][NumCurrentIndices-1] = (int) IndexBase64() - 1;
00646     }
00647   }
00648   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00649 
00650   EPETRA_CHK_ERR(ierr);
00651 
00652   if(CrsGraphData_->ReferenceCount() > 1)
00653     return(1);
00654   else
00655     return(0);
00656 }
00657 
00658 //==============================================================================
00659 template<typename int_type>
00660 int Epetra_CrsGraph::TRemoveGlobalIndices(long long Row) {
00661   int j;
00662   int ierr = 0;
00663 
00664   if(IndicesAreContiguous() || StorageOptimized())
00665     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00666 
00667   if(IndicesAreLocal())
00668     EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
00669   if(CrsGraphData_->CV_ == View)
00670     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00671 
00672   // Normalize row range
00673 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
00674   int locRow = LRID((int) Row);
00675 #else
00676   int locRow = LRID(Row);
00677 #endif
00678 
00679   if(locRow < 0 || locRow >= NumMyBlockRows())
00680     EPETRA_CHK_ERR(-1); // Not in Row range
00681 
00682   Epetra_CrsGraphData::IndexData<int_type>& Data = CrsGraphData_->Data<int_type>();
00683 
00684   if (CrsGraphData_->StaticProfile_) {
00685     int NumIndices = CrsGraphData_->NumIndicesPerRow_[locRow];
00686 
00687     const int_type indexBaseMinusOne = (int_type) IndexBase64() - 1;
00688     for(j = 0; j < NumIndices; j++)
00689       Data.Indices_[locRow][j] = indexBaseMinusOne; // Set to invalid
00690   }
00691   else
00692     Data.SortedEntries_[locRow].entries_.resize(0);
00693 
00694   CrsGraphData_->NumIndicesPerRow_[locRow] = 0;
00695 
00696 
00697   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00698   EPETRA_CHK_ERR(ierr);
00699 
00700   if(CrsGraphData_->ReferenceCount() > 1)
00701     return(1);
00702   else
00703     return(0);
00704 }
00705 
00706 int Epetra_CrsGraph::RemoveGlobalIndices(long long Row)
00707 {
00708   if(RowMap().GlobalIndicesLongLong())
00709 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00710     return TRemoveGlobalIndices<long long>(Row);
00711 #else
00712     throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: ERROR, GlobalIndicesLongLong but no API for it.",-1);
00713 #endif
00714 
00715   if(RowMap().GlobalIndicesInt())
00716 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00717     return TRemoveGlobalIndices<int>(Row);
00718 #else
00719     throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: ERROR, GlobalIndicesInt but no API for it.",-1);
00720 #endif
00721 
00722   throw ReportError("Epetra_CrsGraph::RemoveGlobalIndices: Internal error.", -1);
00723 }
00724 
00725 //==============================================================================
00726 int Epetra_CrsGraph::RemoveMyIndices(int Row)
00727 {
00728   int ierr = 0;
00729 
00730   if(IndicesAreContiguous() || StorageOptimized())
00731     EPETRA_CHK_ERR(-1); // Indices cannot be individually deleted and newed
00732 
00733   if(IndicesAreGlobal())
00734     EPETRA_CHK_ERR(-2); // Cannot remove from a filled graph
00735 
00736   if(CrsGraphData_->CV_ == View)
00737     EPETRA_CHK_ERR(-3); // This is a view only.  Cannot remove entries.
00738 
00739   if(Row < 0 || Row >= NumMyBlockRows())
00740     EPETRA_CHK_ERR(-1); // Not in Row range
00741 
00742   Epetra_CrsGraphData::IndexData<int>& Data = CrsGraphData_->Data<int>();
00743 
00744   if (CrsGraphData_->StaticProfile_) {
00745     int NumIndices = CrsGraphData_->NumIndicesPerRow_[Row];
00746     for(int j = 0; j < NumIndices; j++)
00747       Data.Indices_[Row][j] = -1; // Set to invalid
00748   }
00749   else
00750     Data.SortedEntries_[Row].entries_.resize(0);
00751 
00752   CrsGraphData_->NumIndicesPerRow_[Row] = 0;
00753 
00754   SetGlobalConstantsComputed(false); // No longer have valid global constants.
00755   EPETRA_CHK_ERR(ierr);
00756 
00757   if(CrsGraphData_->ReferenceCount() > 1)
00758     return(1);
00759   else
00760     return(0);
00761 }
00762 
00763 // protected ===================================================================
00764 template<typename int_type>
00765 bool Epetra_CrsGraph::FindGlobalIndexLoc(int LocalRow,
00766            int_type Index,
00767            int Start,
00768            int& Loc) const
00769 {
00770   int NumIndices = NumMyIndices(LocalRow);
00771 
00772   // If we have transformed the column indices, we must map this global Index to local
00773   if(CrsGraphData_->IndicesAreLocal_) {
00774     int* locIndices = Indices(LocalRow);
00775     int locIndex = LCID(Index);
00776     if (CrsGraphData_->Sorted_) {
00777       int insertPoint;
00778       Loc = Epetra_Util_binary_search(locIndex, locIndices, NumIndices, insertPoint);
00779       return( Loc > -1 );
00780     }
00781     else {
00782       int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00783       for(j = 0; j < NumIndices; j++) {
00784         if(j0 >= NumIndices)
00785           j0 = 0; // wrap around
00786       if(locIndices[j0] == locIndex) {
00787           Loc = j0;
00788           return(true);
00789       }
00790       j0++;
00791      }
00792     }
00793   }
00794   else {
00795     int_type* locIndices = TIndices<int_type>(LocalRow);
00796     if (CrsGraphData_->Sorted_) {
00797       int insertPoint;
00798       Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
00799       return( Loc > -1 );
00800     }
00801     else {
00802       int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00803       for(j = 0; j < NumIndices; j++) {
00804         if(j0 >= NumIndices)
00805           j0 = 0; // wrap around
00806       if(locIndices[j0] == Index) {
00807           Loc = j0;
00808           return(true);
00809       }
00810       j0++;
00811      }
00812     }
00813   }
00814 
00815   return(false);
00816 }
00817 
00818 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00819 bool Epetra_CrsGraph::FindGlobalIndexLoc(int LocalRow,
00820            int Index,
00821            int Start,
00822            int& Loc) const
00823 {
00824   if(RowMap().GlobalIndicesInt())
00825   return FindGlobalIndexLoc<int>(LocalRow, Index, Start, Loc);
00826   else
00827   throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc int version called for a graph that is not int.", -1);
00828 }
00829 #endif
00830 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00831 bool Epetra_CrsGraph::FindGlobalIndexLoc(int LocalRow,
00832            long long Index,
00833            int Start,
00834            int& Loc) const
00835 {
00836   if(RowMap().GlobalIndicesLongLong())
00837   return FindGlobalIndexLoc<long long>(LocalRow, Index, Start, Loc);
00838   else
00839   throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc long long version called for a graph that is not long long.", -1);
00840 }
00841 #endif
00842 
00843 // protected ===================================================================
00844 template<typename int_type>
00845 bool Epetra_CrsGraph::FindGlobalIndexLoc(int NumIndices,
00846            const int_type* indices,
00847            int_type Index,
00848            int Start,
00849            int& Loc) const
00850 {
00851   // If we have transformed the column indices, we must map this global Index to local
00852   if(CrsGraphData_->IndicesAreLocal_) {
00853     Index = LCID(Index);
00854   }
00855 
00856   if (CrsGraphData_->Sorted_) {
00857     int insertPoint;
00858     Loc = Epetra_Util_binary_search(Index, indices, NumIndices, insertPoint);
00859     return( Loc > -1 );
00860   }
00861   else {
00862     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00863     for(j = 0; j < NumIndices; j++) {
00864       if(j0 >= NumIndices)
00865   j0 = 0; // wrap around
00866       if(indices[j0] == Index) {
00867   Loc = j0;
00868   return(true);
00869       }
00870       j0++;
00871     }
00872   }
00873   return(false);
00874 }
00875 
00876 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00877 bool Epetra_CrsGraph::FindGlobalIndexLoc(int NumIndices,
00878            const int* indices,
00879            int Index,
00880            int Start,
00881            int& Loc) const
00882 {
00883   if(RowMap().GlobalIndicesInt())
00884   return FindGlobalIndexLoc<int>(NumIndices, indices, Index, Start, Loc);
00885   else
00886   throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc int version called for a graph that is not int.", -1);
00887 }
00888 #endif
00889 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00890 bool Epetra_CrsGraph::FindGlobalIndexLoc(int NumIndices,
00891            const long long* indices,
00892            long long Index,
00893            int Start,
00894            int& Loc) const
00895 {
00896   if(RowMap().GlobalIndicesLongLong())
00897   return FindGlobalIndexLoc<long long>(NumIndices, indices, Index, Start, Loc);
00898   else
00899   throw ReportError("Epetra_CrsGraph::FindGlobalIndexLoc long long version called for a graph that is not long long.", -1);
00900 }
00901 #endif
00902 
00903 // protected ===================================================================
00904 bool Epetra_CrsGraph::FindMyIndexLoc(int LocalRow,
00905              int Index,
00906              int Start,
00907              int& Loc) const
00908 {
00909   int NumIndices = NumMyIndices(LocalRow);
00910   int* locIndices = Indices(LocalRow);
00911 
00912   if(!CrsGraphData_->IndicesAreLocal_) {
00913     throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
00914   }
00915 
00916   if (CrsGraphData_->Sorted_) {
00917     int insertPoint;
00918     Loc = Epetra_Util_binary_search(Index, locIndices, NumIndices, insertPoint);
00919     return( Loc > -1 );
00920   }
00921   else {
00922     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00923     for(j = 0; j < NumIndices; j++) {
00924       if(j0 >= NumIndices)
00925   j0 = 0; // wrap around
00926       if(locIndices[j0] == Index) {
00927   Loc = j0;
00928   return(true);
00929       }
00930       j0++;
00931     }
00932   }
00933   return(false);
00934 }
00935 
00936 // protected ===================================================================
00937 bool Epetra_CrsGraph::FindMyIndexLoc(int NumIndices,
00938              const int* indices,
00939              int Index,
00940              int Start,
00941              int& Loc) const
00942 {
00943   if(!CrsGraphData_->IndicesAreLocal_) {
00944     throw ReportError("Epetra_CrsGraph::FindMyIndexLoc", -1);// Indices must be local
00945   }
00946 
00947   if (CrsGraphData_->Sorted_) {
00948     int insertPoint;
00949     Loc = Epetra_Util_binary_search(Index, indices, NumIndices, insertPoint);
00950     return( Loc > -1 );
00951   }
00952   else {
00953     int j, j0 = Start; // Start search at index Start (must be >= 0 and < NumIndices)
00954     for(j = 0; j < NumIndices; j++) {
00955       if(j0 >= NumIndices)
00956   j0 = 0; // wrap around
00957       if(indices[j0] == Index) {
00958   Loc = j0;
00959   return(true);
00960       }
00961       j0++;
00962     }
00963   }
00964   return(false);
00965 }
00966 
00967 //==============================================================================
00968 int Epetra_CrsGraph::FillComplete() {
00969   EPETRA_CHK_ERR(FillComplete(RowMap(), RowMap()));
00970   return(0); // uses FillComplete(...)'s shared-ownership test.
00971 }
00972 
00973 //==============================================================================
00974 int Epetra_CrsGraph::FillComplete(const Epetra_BlockMap& domainMap, const Epetra_BlockMap& rangeMap) {
00975   if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
00976      throw ReportError("Epetra_CrsGraph::FillComplete: cannot be called with different indices types for domainMap and rangeMap", -1);
00977 
00978   if(!RowMap().GlobalIndicesTypeMatch(domainMap))
00979     throw ReportError("Epetra_CrsGraph::FillComplete: cannot be called with different indices types for row map and incoming rangeMap", -1);
00980 
00981   CrsGraphData_->DomainMap_ = domainMap;
00982   CrsGraphData_->RangeMap_ = rangeMap;
00983 
00984   MakeIndicesLocal(domainMap, rangeMap); // Convert global indices to local indices to on each processor
00985   SortIndices();  // Sort column entries from smallest to largest
00986   RemoveRedundantIndices(); // Get rid of any redundant index values
00987   DetermineTriangular(); //determine if lower/upper triangular
00988   CrsGraphData_->MakeImportExport(); // Build Import or Export objects
00989   ComputeGlobalConstants(); // Compute constants that require communication
00990   SetFilled(true);
00991 
00992   if(CrsGraphData_->ReferenceCount() > 1)
00993     return(1);
00994   else
00995     return(0);
00996 }
00997 
00998 //==============================================================================
00999 int Epetra_CrsGraph::TransformToLocal() {
01000   return(FillComplete());
01001 }
01002 
01003 //==============================================================================
01004 int Epetra_CrsGraph::TransformToLocal(const Epetra_BlockMap* domainMap, const Epetra_BlockMap* rangeMap) {
01005   return(FillComplete(*domainMap, *rangeMap));
01006 }
01007 
01008 // private =====================================================================
01009 int Epetra_CrsGraph::ComputeGlobalConstants()
01010 {
01011   if(GlobalConstantsComputed())
01012     return(0);
01013 
01014 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01015   Epetra_LongLongSerialDenseVector
01016 #else
01017   Epetra_IntSerialDenseVector
01018 #endif
01019      tempvec(8); // Temp space
01020 
01021   const int numMyBlockRows = NumMyBlockRows();
01022 
01023   CrsGraphData_->NumMyEntries_ = 0; // Compute Number of Nonzero entries and max
01024   CrsGraphData_->MaxNumIndices_ = 0;
01025   {for(int i = 0; i < numMyBlockRows; i++) {
01026     CrsGraphData_->NumMyEntries_ += CrsGraphData_->NumIndicesPerRow_[i];
01027     CrsGraphData_->MaxNumIndices_ = EPETRA_MAX(CrsGraphData_->MaxNumIndices_,CrsGraphData_->NumIndicesPerRow_[i]);
01028   }}
01029 
01030   // Case 1:  Constant block size (including blocksize = 1)
01031   if(RowMap().ConstantElementSize() && ColMap().ConstantElementSize() && RowMap().ElementSize() == ColMap().ElementSize()) {
01032     // Jim Westfall reported a fix on 22 June 2013 where the second two conditions
01033     // above are necessary.  The added conditions check for the case when the row map
01034     // and col map are constant but different as possible with VBR sub matrices used
01035     // in global assemble
01036 
01037     tempvec[0] = CrsGraphData_->NumMyEntries_;
01038     tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
01039 
01040     Comm().SumAll(&tempvec[0], &tempvec[2], 2);
01041     int tmp_MaxNumIndices = CrsGraphData_->MaxNumIndices_;
01042     Comm().MaxAll(&tmp_MaxNumIndices, &CrsGraphData_->GlobalMaxNumIndices_, 1);
01043 
01044     CrsGraphData_->NumGlobalEntries_ = tempvec[2];
01045     CrsGraphData_->NumGlobalBlockDiagonals_ = tempvec[3];
01046 
01047     int RowElementSize = RowMap().MaxElementSize();
01048     int ColElementSize = RowElementSize;
01049     CrsGraphData_->NumGlobalDiagonals_ = tempvec[3] * RowElementSize;
01050     CrsGraphData_->NumMyNonzeros_ = CrsGraphData_->NumMyEntries_ * RowElementSize * ColElementSize;
01051     CrsGraphData_->NumGlobalNonzeros_ = CrsGraphData_->NumGlobalEntries_ * RowElementSize * ColElementSize;
01052     CrsGraphData_->MaxNumNonzeros_ = CrsGraphData_->MaxNumIndices_ * RowElementSize * ColElementSize;
01053     CrsGraphData_->GlobalMaxNumNonzeros_ = CrsGraphData_->GlobalMaxNumIndices_ * RowElementSize * ColElementSize;
01054   }
01055 
01056   // Case 2:  Variable block size (more work)
01057   else {
01058     CrsGraphData_->NumMyNonzeros_ = 0;  // We will count the number of nonzeros here
01059     CrsGraphData_->MaxNumNonzeros_ = 0;  // We will determine the max number of nonzeros in any one block row
01060     int* RowElementSizeList = RowMap().ElementSizeList();
01061     int* ColElementSizeList = RowElementSizeList;
01062     if(Importer() != 0)
01063       ColElementSizeList = ColMap().ElementSizeList();
01064       const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
01065     for(int i = 0; i < numMyBlockRows; i++){
01066       int NumEntries = CrsGraphData_->NumIndicesPerRow_[i];
01067       int* indices = intData.Indices_[i];
01068       if(NumEntries > 0) {
01069   int CurNumNonzeros = 0;
01070   int RowDim = RowElementSizeList[i];
01071   for(int j = 0; j < NumEntries; j++) {
01072     int ColDim = ColElementSizeList[indices[j]];
01073     CurNumNonzeros += RowDim*ColDim;
01074     CrsGraphData_->MaxColDim_ = EPETRA_MAX(CrsGraphData_->MaxColDim_, ColDim);
01075   }
01076   CrsGraphData_->MaxNumNonzeros_ = EPETRA_MAX(CrsGraphData_->MaxNumNonzeros_, CurNumNonzeros);
01077   CrsGraphData_->NumMyNonzeros_ += CurNumNonzeros;
01078       }
01079     }
01080 
01081     // Sum Up all nonzeros
01082 
01083     tempvec[0] = CrsGraphData_->NumMyEntries_;
01084     tempvec[1] = CrsGraphData_->NumMyBlockDiagonals_;
01085     tempvec[2] = CrsGraphData_->NumMyDiagonals_;
01086     tempvec[3] = CrsGraphData_->NumMyNonzeros_;
01087 
01088     Comm().SumAll(&tempvec[0], &tempvec[4], 4);
01089 
01090     CrsGraphData_->NumGlobalEntries_ = tempvec[4];
01091     CrsGraphData_->NumGlobalBlockDiagonals_ = tempvec[5];
01092     CrsGraphData_->NumGlobalDiagonals_ = tempvec[6];
01093     CrsGraphData_->NumGlobalNonzeros_ = tempvec[7];
01094 
01095     tempvec[0] = CrsGraphData_->MaxNumIndices_;
01096     tempvec[1] = CrsGraphData_->MaxNumNonzeros_;
01097 
01098     Comm().MaxAll(&tempvec[0], &tempvec[2], 2);
01099 
01100     CrsGraphData_->GlobalMaxNumIndices_ = (int) tempvec[2];
01101     CrsGraphData_->GlobalMaxNumNonzeros_ = (int) tempvec[3];
01102   }
01103 
01104   CrsGraphData_->NumGlobalRows_ = CrsGraphData_->RangeMap_.NumGlobalPoints64();
01105   CrsGraphData_->NumGlobalCols_ = DomainMap().NumGlobalPoints64();
01106 
01107   CrsGraphData_->GlobalConstantsComputed_ = true;
01108 
01109   return(0);
01110 }
01111 
01112 void epetra_shellsort(int* list, int length)
01113 {
01114   int i, j, j2, temp, istep;
01115   unsigned step;
01116 
01117   step = length/2;
01118   while (step > 0)
01119   {
01120     for (i=step; i < length; i++)
01121     {
01122       istep = step;
01123       j = i;
01124       j2 = j-istep;
01125       temp = list[i];
01126       if (list[j2] > temp) {
01127         while ((j >= istep) && (list[j2] > temp))
01128         {
01129           list[j] = list[j2];
01130           j = j2;
01131           j2 -= istep;
01132         }
01133         list[j] = temp;
01134       }
01135     }
01136 
01137     if (step == 2)
01138       step = 1;
01139     else
01140       step = (int) (step / 2.2);
01141   }
01142 }
01143 
01144 //==============================================================================
01145 int Epetra_CrsGraph::SortIndices() {
01146   if(IndicesAreGlobal())
01147     EPETRA_CHK_ERR(-1);
01148   if(Sorted())
01149     return(0);
01150 
01151   // For each row, sort column entries from smallest to largest.
01152   // Use shell sort, which is fast if indices are already sorted.
01153 
01154   const int numMyBlockRows = NumMyBlockRows();
01155   const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
01156   for(int i = 0; i < numMyBlockRows; i++){
01157     int n = CrsGraphData_->NumIndicesPerRow_[i];
01158     int* const list = intData.Indices_[i];
01159 
01160     epetra_shellsort(list, n);
01161 //    int m = n/2;
01162 
01163 //    while(m > 0) {
01164  //     int max = n - m;
01165 //      for(int j = 0; j < max; j++) {
01166 //        int k = j;
01167 //        while(k>-1) {
01168 //    if(list[k+m] >= list[k])
01169 //      break;
01170 //    int itemp = list[k+m];
01171 //    list[k+m] = list[k];
01172 //    list[k] = itemp;
01173 //          k-=m;
01174 //  }
01175 //      }
01176 //      m = m/2;
01177 //    }
01178   }
01179   SetSorted(true);
01180 
01181   if(CrsGraphData_->ReferenceCount() > 1)
01182     return(1);
01183   else
01184     return(0);
01185 }
01186 
01187 void epetra_crsgraph_compress_out_duplicates(int len, int* list, int& newlen)
01188 {
01189   //
01190   //This function runs the array ('list') checking for
01191   //duplicate entries. Any duplicates that are found are
01192   //removed by sliding subsequent data down in the array,
01193   //over-writing the duplicates. Finally, the new length
01194   //of the array (i.e., the number of unique entries) is
01195   //placed in the output parameter 'newlen'. The array is
01196   //**not** re-allocated.
01197   //
01199   //Important assumption: The array contents are assumed to
01200   //be sorted before this function is called. If the array
01201   //contents are not sorted, then the behavior of this
01202   //function is undefined.
01204   //
01205 
01206   if (len < 2) return;
01207 
01208   int* ptr0 = &list[0];
01209   int* ptr1 = &list[1];
01210 
01211   int* ptr_end = &list[len-1];
01212 
01213   while(*ptr0 != *ptr1 && ptr1 < ptr_end) {
01214     ++ptr0;
01215     ++ptr1;
01216   }
01217 
01218   if (ptr1 < ptr_end) {
01219     //if ptr1 < ptr_end we've found a duplicate...
01220 
01221     ++ptr0;
01222     ++ptr1;
01223 
01224     while(*ptr0 == *ptr1 && ptr1 < ptr_end) ++ptr1;
01225 
01226     while(ptr1 < ptr_end) {
01227 
01228       int val = *ptr1++;
01229 
01230       while(val == *ptr1 && ptr1 < ptr_end) {
01231         ++ptr1;
01232       }
01233 
01234       *ptr0++ = val;
01235     }
01236 
01237     if (*(ptr0-1) != *ptr1) *ptr0++ = *ptr1;
01238 
01239     int num_removed = (int)(ptr_end - ptr0 + 1);
01240     newlen = len - num_removed;
01241   }
01242   else {
01243     if (*ptr0 == *ptr1) newlen = len - 1;
01244     else newlen = len;
01245   }
01246 }
01247 
01248 //==============================================================================
01249 int Epetra_CrsGraph::RemoveRedundantIndices()
01250 {
01251   if(!Sorted())
01252     EPETRA_CHK_ERR(-1);  // Must have sorted index set
01253   if(IndicesAreGlobal())
01254     EPETRA_CHK_ERR(-2); // Indices must be local
01255 
01256   // Note:  This function assumes that SortIndices was already called.
01257   // For each row, remove column indices that are repeated.
01258 
01259   const int numMyBlockRows = NumMyBlockRows();
01260   bool found_redundancies = false;
01261 
01262   if(NoRedundancies() == false) {
01263     int* numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
01264     Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
01265     for(int i=0; i<numMyBlockRows; ++i) {
01266       int NumIndices = numIndicesPerRow[i];
01267       int* col_indices = this->Indices(i);
01268 
01269       if(NumIndices > 1) {
01270         epetra_crsgraph_compress_out_duplicates(NumIndices, col_indices,
01271                                                 numIndicesPerRow[i]);
01272       }
01273       if (NumIndices != numIndicesPerRow[i]) {
01274           found_redundancies = true;
01275       }
01276     }
01277     if (found_redundancies && !CrsGraphData_->StaticProfile_)
01278     {
01279       for(int i=0; i<numMyBlockRows; ++i) {
01280         int* col_indices = this->Indices(i);
01281 
01282         // update vector size and address in memory
01283         intData.SortedEntries_[i].entries_.assign(col_indices, col_indices+numIndicesPerRow[i]);
01284         if (numIndicesPerRow[i] > 0) {
01285           intData.Indices_[i] = &(intData.SortedEntries_[i].entries_[0]);
01286         }
01287         else {
01288           intData.Indices_[i] = NULL;
01289         }
01290       }
01291     }
01292   }
01293 
01294   SetNoRedundancies(true);
01295   return 0;
01296 }
01297 
01298 //==============================================================================
01299 int Epetra_CrsGraph::DetermineTriangular()
01300 {
01301   // determine if graph is upper or lower triangular or has no diagonal
01302 
01303   if(!Sorted())
01304     EPETRA_CHK_ERR(-1);  // Must have sorted index set
01305   if(IndicesAreGlobal())
01306     EPETRA_CHK_ERR(-2); // Indices must be local
01307 
01308   CrsGraphData_->NumMyDiagonals_ = 0;
01309   CrsGraphData_->NumMyBlockDiagonals_ = 0;
01310 
01311   const Epetra_BlockMap& rowMap = RowMap();
01312   const Epetra_BlockMap& colMap = ColMap();
01313 
01314   const int numMyBlockRows = NumMyBlockRows();
01315 
01316   for(int i = 0; i < numMyBlockRows; i++) {
01317     int NumIndices = NumMyIndices(i);
01318     if(NumIndices > 0) {
01319 #if defined(EPETRA_NO_64BIT_GLOBAL_INDICES) && !defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
01320       int ig = rowMap.GID(i);
01321 #else
01322       long long ig = rowMap.GID64(i);
01323 #endif
01324       int* col_indices = this->Indices(i);
01325 
01326       int jl_0 = col_indices[0];
01327       int jl_n = col_indices[NumIndices-1];
01328 
01329       if(jl_n > i) CrsGraphData_->LowerTriangular_ = false;
01330       if(jl_0 < i) CrsGraphData_->UpperTriangular_ = false;
01331 
01332       //jl will be the local-index for the diagonal that we
01333       //want to search for.
01334       int jl = colMap.LID(ig);
01335 
01336       int insertPoint = -1;
01337       if (Epetra_Util_binary_search(jl, col_indices, NumIndices, insertPoint)>-1) {
01338         CrsGraphData_->NumMyBlockDiagonals_++;
01339         CrsGraphData_->NumMyDiagonals_ += rowMap.ElementSize(i);
01340       }
01341     }
01342   }
01343 
01344   CrsGraphData_->NoDiagonal_ = (CrsGraphData_->NumMyBlockDiagonals_ == 0);
01345 
01346   if(CrsGraphData_->ReferenceCount() > 1)
01347     return(1);
01348   else
01349     return(0);
01350 }
01351 
01352 // private =====================================================================
01353 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01354 int Epetra_CrsGraph::MakeColMap_int(const Epetra_BlockMap& domainMap,
01355         const Epetra_BlockMap& rangeMap)
01356 {
01357   (void)rangeMap;
01358   int i;
01359   int j;
01360 
01361   if(CrsGraphData_->HaveColMap_)
01362     return(0); // Already have a Column Map
01363 
01364   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01365   if(IndicesAreLocal())
01366     EPETRA_CHK_ERR(-1); // Return error: Indices must be global
01367 
01368   // Scan all column indices and sort into two groups:
01369   // Local:  those whose GID matches a GID of the domain map on this processor and
01370   // Remote: All others.
01371   int numDomainElements = domainMap.NumMyElements();
01372   bool * LocalGIDs  = 0;
01373   if (numDomainElements>0) LocalGIDs  = new bool[numDomainElements];
01374   for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
01375 
01376   // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
01377   // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
01378   // and the number of block rows.
01379   const int numMyBlockRows = NumMyBlockRows();
01380   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
01381   //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
01382   Epetra_HashTable<int> RemoteGIDs(hashsize);
01383   Epetra_HashTable<int> RemoteGIDList(hashsize);
01384 
01385   int NumLocalColGIDs = 0;
01386   int NumRemoteColGIDs = 0;
01387   const Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
01388   for(i = 0; i < numMyBlockRows; i++) {
01389     const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01390     int* ColIndices = intData.Indices_[i];
01391     for(j = 0; j < NumIndices; j++) {
01392       int GID = ColIndices[j];
01393       // Check if GID matches a row GID
01394       int LID = domainMap.LID(GID);
01395       if(LID != -1) {
01396   bool alreadyFound = LocalGIDs[LID];
01397   if (!alreadyFound) {
01398           LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
01399           NumLocalColGIDs++;
01400   }
01401       }
01402       else {
01403   if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
01404     RemoteGIDs.Add(GID, NumRemoteColGIDs);
01405     RemoteGIDList.Add(NumRemoteColGIDs++, GID);
01406   }
01407       }
01408     }
01409   }
01410 
01411   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
01412   if (domainMap.Comm().NumProc()==1) {
01413 
01414     if (NumRemoteColGIDs!=0) {
01415       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
01416     }
01417     if (NumLocalColGIDs==numDomainElements) {
01418       CrsGraphData_->ColMap_ = domainMap;
01419       CrsGraphData_->HaveColMap_ = true;
01420       if (LocalGIDs!=0) delete [] LocalGIDs;
01421       return(0);
01422     }
01423   }
01424 
01425   // Now build integer array containing column GIDs
01426   // Build back end, containing remote GIDs, first
01427   int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
01428   Epetra_IntSerialDenseVector ColIndices;
01429   if(numMyBlockCols > 0)
01430     ColIndices.Size(numMyBlockCols);
01431 
01432   int* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01433 
01434   for(i = 0; i < NumRemoteColGIDs; i++)
01435     RemoteColIndices[i] = RemoteGIDList.Get(i);
01436 
01437   int NLists = 1;
01438   Epetra_IntSerialDenseVector PIDList;
01439   Epetra_IntSerialDenseVector SizeList;
01440   int* RemoteSizeList = 0;
01441   bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
01442 
01443   if(NumRemoteColGIDs > 0)
01444     PIDList.Size(NumRemoteColGIDs);
01445 
01446   if(DoSizes) {
01447     if(numMyBlockCols > 0)
01448       SizeList.Size(numMyBlockCols);
01449     RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
01450     NLists++;
01451   }
01452   EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01453 
01454   // Sort External column indices so that all columns coming from a given remote processor are contiguous
01455 
01456   Epetra_Util Util;
01457   int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
01458   SortLists[0] = RemoteColIndices;
01459   SortLists[1] = RemoteSizeList;
01460   Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, SortLists, 0, 0);
01461   if (CrsGraphData_->SortGhostsAssociatedWithEachProcessor_) {
01462     // Sort external column indices so that columns from a given remote processor are not only contiguous
01463     // but also in ascending order. NOTE: I don't know if the number of externals associated
01464     // with a given remote processor is known at this point ... so I count them here.
01465 
01466     NLists--;
01467     int StartCurrent, StartNext;
01468     StartCurrent = 0; StartNext = 1;
01469     while ( StartNext < NumRemoteColGIDs ) {
01470       if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
01471       else {
01472         if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01473         Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
01474         StartCurrent = StartNext; StartNext++;
01475       }
01476     }
01477     if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01478     Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
01479   }
01480 
01481   // Now fill front end. Two cases:
01482   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
01483   //     can simply read the domain GIDs into the front part of ColIndices, otherwise
01484   // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
01485   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
01486 
01487   if(NumLocalColGIDs == domainMap.NumMyElements()) {
01488     domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
01489     if(DoSizes)
01490       domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
01491   }
01492   else {
01493     int NumMyElements = domainMap.NumMyElements();
01494     int* MyGlobalElements = domainMap.MyGlobalElements();
01495     int* ElementSizeList = 0;
01496     if(DoSizes)
01497       ElementSizeList = domainMap.ElementSizeList();
01498     int NumLocalAgain = 0;
01499     for(i = 0; i < NumMyElements; i++) {
01500       if(LocalGIDs[i]) {
01501   if(DoSizes)
01502     SizeList[NumLocalAgain] = ElementSizeList[i];
01503   ColIndices[NumLocalAgain++] = MyGlobalElements[i];
01504       }
01505     }
01506     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
01507   }
01508 
01509   // Done with this array
01510   if (LocalGIDs!=0) delete [] LocalGIDs;
01511 
01512 
01513   // Make Column map with same element sizes as Domain map
01514 
01515   if(domainMap.MaxElementSize() == 1) { // Simple map
01516     Epetra_Map temp((int) -1, numMyBlockCols, ColIndices.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
01517     CrsGraphData_->ColMap_ = temp;
01518   }
01519   else if(domainMap.ConstantElementSize()) { // Constant Block size map
01520     Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(), (int) domainMap.IndexBase64(), domainMap.Comm());
01521     CrsGraphData_->ColMap_ = temp;
01522   }
01523   else { // Most general case where block size is variable.
01524     Epetra_BlockMap temp((int) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), (int) domainMap.IndexBase64(), domainMap.Comm());
01525     CrsGraphData_->ColMap_ = temp;
01526   }
01527   CrsGraphData_->HaveColMap_ = true;
01528 
01529   return(0);
01530 }
01531 #endif
01532 //==============================================================================
01533 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01534 int Epetra_CrsGraph::MakeColMap_LL(const Epetra_BlockMap& domainMap,
01535         const Epetra_BlockMap& rangeMap)
01536 {
01537   (void)rangeMap;
01538   int i;
01539   int j;
01540 
01541   if(CrsGraphData_->HaveColMap_)
01542     return(0); // Already have a Column Map
01543 
01544   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01545   if(IndicesAreLocal())
01546     EPETRA_CHK_ERR(-1); // Return error: Indices must be global
01547 
01548   // Scan all column indices and sort into two groups:
01549   // Local:  those whose GID matches a GID of the domain map on this processor and
01550   // Remote: All others.
01551   int numDomainElements = domainMap.NumMyElements();
01552   bool * LocalGIDs  = 0;
01553   if (numDomainElements>0) LocalGIDs  = new bool[numDomainElements];
01554   for (i=0; i<numDomainElements; i++) LocalGIDs[i] = false; // Assume domain GIDs are not local
01555 
01556   // In principle it is good to have RemoteGIDs and RemotGIDList be as long as the number of remote GIDs
01557   // on this processor, but this would require two passes through the column IDs, so we make it the max of 100
01558   // and the number of block rows.
01559   const int numMyBlockRows = NumMyBlockRows();
01560   int  hashsize = numMyBlockRows; if (hashsize < 100) hashsize = 100;
01561   //cout << "numMyBlockRows = " << numMyBlockRows << " hashsize = " << hashsize << std::endl;
01562   Epetra_HashTable<int> RemoteGIDs(hashsize);
01563   Epetra_HashTable<long long> RemoteGIDList(hashsize);
01564 
01565   int NumLocalColGIDs = 0;
01566   int NumRemoteColGIDs = 0;
01567 
01568   if(IndicesAreLocal())
01569   {
01570     Epetra_CrsGraphData::IndexData<int>& intData = CrsGraphData_->Data<int>();
01571 
01572     for(i = 0; i < numMyBlockRows; i++) {
01573     const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01574     int* ColIndices = intData.Indices_[i];
01575     for(j = 0; j < NumIndices; j++) {
01576       int GID = ColIndices[j];
01577       // Check if GID matches a row GID
01578       int LID = domainMap.LID(GID);
01579       if(LID != -1) {
01580     bool alreadyFound = LocalGIDs[LID];
01581     if (!alreadyFound) {
01582         LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
01583         NumLocalColGIDs++;
01584     }
01585       }
01586       else {
01587     if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
01588       RemoteGIDs.Add(GID, NumRemoteColGIDs);
01589       RemoteGIDList.Add(NumRemoteColGIDs++, GID);
01590     }
01591       }
01592     }
01593     }
01594   }
01595   else if(IndicesAreGlobal())
01596   {
01597     Epetra_CrsGraphData::IndexData<long long>& LLData = CrsGraphData_->Data<long long>();
01598 
01599     for(i = 0; i < numMyBlockRows; i++) {
01600     const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01601     long long* ColIndices = LLData.Indices_[i];
01602     for(j = 0; j < NumIndices; j++) {
01603       long long GID = ColIndices[j];
01604       // Check if GID matches a row GID
01605       int LID = domainMap.LID(GID);
01606       if(LID != -1) {
01607     bool alreadyFound = LocalGIDs[LID];
01608     if (!alreadyFound) {
01609         LocalGIDs[LID] = true; // There is a column in the graph associated with this domain map GID
01610         NumLocalColGIDs++;
01611     }
01612       }
01613       else {
01614     if(RemoteGIDs.Get(GID) == -1) { // This means its a new remote GID
01615       RemoteGIDs.Add(GID, NumRemoteColGIDs);
01616       RemoteGIDList.Add(NumRemoteColGIDs++, GID);
01617     }
01618       }
01619     }
01620     }
01621   }
01622 
01623   // Possible short-circuit:  If all domain map GIDs are present as column indices, then set ColMap=domainMap and quit
01624   if (domainMap.Comm().NumProc()==1) {
01625 
01626     if (NumRemoteColGIDs!=0) {
01627       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
01628     }
01629     if (NumLocalColGIDs==numDomainElements) {
01630       CrsGraphData_->ColMap_ = domainMap;
01631       CrsGraphData_->HaveColMap_ = true;
01632       if (LocalGIDs!=0) delete [] LocalGIDs;
01633       return(0);
01634     }
01635   }
01636 
01637   // Now build integer array containing column GIDs
01638   // Build back end, containing remote GIDs, first
01639   int numMyBlockCols = NumLocalColGIDs + NumRemoteColGIDs;
01640   Epetra_LongLongSerialDenseVector ColIndices;
01641   if(numMyBlockCols > 0)
01642     ColIndices.Size(numMyBlockCols);
01643 
01644   long long* RemoteColIndices = ColIndices.Values() + NumLocalColGIDs; // Points to back end of ColIndices
01645 
01646   for(i = 0; i < NumRemoteColGIDs; i++)
01647     RemoteColIndices[i] = RemoteGIDList.Get(i);
01648 
01649   int NLists = 0;
01650   Epetra_IntSerialDenseVector PIDList;
01651   Epetra_IntSerialDenseVector SizeList;
01652   int* RemoteSizeList = 0;
01653   bool DoSizes = !domainMap.ConstantElementSize(); // If not constant element size, then we must exchange
01654 
01655   if(NumRemoteColGIDs > 0)
01656     PIDList.Size(NumRemoteColGIDs);
01657 
01658   if(DoSizes) {
01659     if(numMyBlockCols > 0)
01660       SizeList.Size(numMyBlockCols);
01661     RemoteSizeList = SizeList.Values() + NumLocalColGIDs;
01662     NLists++;
01663   }
01664   EPETRA_CHK_ERR(domainMap.RemoteIDList(NumRemoteColGIDs, RemoteColIndices, PIDList.Values(), 0, RemoteSizeList));
01665 
01666   // Sort External column indices so that all columns coming from a given remote processor are contiguous
01667 
01668   Epetra_Util Util;
01669   //int* SortLists[2]; // this array is allocated on the stack, and so we won't need to delete it.bb
01670   //SortLists[0] = RemoteColIndices;
01671   //SortLists[1] = RemoteSizeList;
01672   Util.Sort(true, NumRemoteColGIDs, PIDList.Values(), 0, 0, NLists, &RemoteSizeList, 1, &RemoteColIndices);
01673   if (CrsGraphData_->SortGhostsAssociatedWithEachProcessor_) {
01674     // Sort external column indices so that columns from a given remote processor are not only contiguous
01675     // but also in ascending order. NOTE: I don't know if the number of externals associated
01676     // with a given remote processor is known at this point ... so I count them here.
01677 
01678   int* SortLists[1] = {0};
01679 
01680     int StartCurrent, StartNext;
01681     StartCurrent = 0; StartNext = 1;
01682     while ( StartNext < NumRemoteColGIDs ) {
01683       if ((PIDList.Values())[StartNext]==(PIDList.Values())[StartNext-1]) StartNext++;
01684       else {
01685         if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01686         Util.Sort(true,StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]),0,0,NLists,SortLists, 0, 0);
01687         StartCurrent = StartNext; StartNext++;
01688       }
01689     }
01690     if(DoSizes) SortLists[0] = &(RemoteSizeList[StartCurrent]);
01691     Util.Sort(true, StartNext-StartCurrent, &(RemoteColIndices[StartCurrent]), 0, 0, NLists, SortLists, 0, 0);
01692   }
01693 
01694   // Now fill front end. Two cases:
01695   // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
01696   //     can simply read the domain GIDs into the front part of ColIndices, otherwise
01697   // (2) We step through the GIDs of the domainMap, checking to see if each domain GID is a column GID.
01698   //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
01699 
01700   if(NumLocalColGIDs == domainMap.NumMyElements()) {
01701     domainMap.MyGlobalElements(ColIndices.Values()); // Load Global Indices into first numMyBlockCols elements column GID list
01702     if(DoSizes)
01703       domainMap.ElementSizeList(SizeList.Values()); // Load ElementSizeList too
01704   }
01705   else {
01706     int NumMyElements = domainMap.NumMyElements();
01707     long long* MyGlobalElements = domainMap.MyGlobalElements64();
01708     int* ElementSizeList = 0;
01709     if(DoSizes)
01710       ElementSizeList = domainMap.ElementSizeList();
01711     int NumLocalAgain = 0;
01712     for(i = 0; i < NumMyElements; i++) {
01713       if(LocalGIDs[i]) {
01714   if(DoSizes)
01715     SizeList[NumLocalAgain] = ElementSizeList[i];
01716   ColIndices[NumLocalAgain++] = MyGlobalElements[i];
01717       }
01718     }
01719     assert(NumLocalAgain==NumLocalColGIDs); // Sanity test
01720   }
01721 
01722   // Done with this array
01723   if (LocalGIDs!=0) delete [] LocalGIDs;
01724 
01725 
01726   // Make Column map with same element sizes as Domain map
01727 
01728   if(domainMap.MaxElementSize() == 1) { // Simple map
01729   Epetra_Map temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.IndexBase64(), domainMap.Comm());
01730   CrsGraphData_->ColMap_ = temp;
01731   }
01732   else if(domainMap.ConstantElementSize()) { // Constant Block size map
01733   Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), domainMap.MaxElementSize(),domainMap.IndexBase64(), domainMap.Comm());
01734   CrsGraphData_->ColMap_ = temp;
01735   }
01736   else { // Most general case where block size is variable.
01737   Epetra_BlockMap temp((long long) -1, numMyBlockCols, ColIndices.Values(), SizeList.Values(), domainMap.IndexBase64(), domainMap.Comm());
01738   CrsGraphData_->ColMap_ = temp;
01739   }
01740   CrsGraphData_->HaveColMap_ = true;
01741 
01742   return(0);
01743 }
01744 #endif
01745 
01746 int Epetra_CrsGraph::MakeColMap(const Epetra_BlockMap& domainMap,
01747         const Epetra_BlockMap& rangeMap)
01748 {
01749   if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
01750     throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for domainMap and rangeMap", -1);
01751 
01752   if(!RowMap().GlobalIndicesTypeMatch(domainMap))
01753     throw ReportError("Epetra_CrsGraph::MakeColMap: cannot be called with different indices types for row map and incoming rangeMap", -1);
01754 
01755   if(RowMap().GlobalIndicesInt())
01756 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01757     return MakeColMap_int(domainMap, rangeMap);
01758 #else
01759     throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesInt but no API for it.",-1);
01760 #endif
01761 
01762   if(RowMap().GlobalIndicesLongLong())
01763 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01764     return MakeColMap_LL(domainMap, rangeMap);
01765 #else
01766     throw ReportError("Epetra_CrsGraph::MakeColMap: ERROR, GlobalIndicesLongLong but no API for it.",-1);
01767 #endif
01768 
01769   throw ReportError("Epetra_CrsGraph::MakeColMap: Internal error, unable to determine global index type of maps", -1);
01770 }
01771 
01772 // protected ===================================================================
01773 int Epetra_CrsGraph::MakeIndicesLocal(const Epetra_BlockMap& domainMap, const Epetra_BlockMap& rangeMap) {
01774   if(!domainMap.GlobalIndicesTypeMatch(rangeMap))
01775      throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for domainMap and rangeMap", -1);
01776 
01777   if(!RowMap().GlobalIndicesTypeMatch(domainMap))
01778     throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: cannot be called with different indices types for row map and incoming rangeMap", -1);
01779 
01780   ComputeIndexState(); // Update index state by checking IndicesAreLocal/Global on all PEs
01781   if(IndicesAreLocal() && IndicesAreGlobal())
01782     EPETRA_CHK_ERR(-1); // Return error: Indices must not be both local and global
01783 
01784   MakeColMap(domainMap, rangeMap); // If user has not prescribed column map, create one from indices
01785   const Epetra_BlockMap& colmap = ColMap();
01786 
01787   // Store number of local columns
01788   CrsGraphData_->NumMyCols_ = ColMap().NumMyPoints();
01789   CrsGraphData_->NumMyBlockCols_ = ColMap().NumMyElements();
01790 
01791   // Transform indices to local index space
01792   const int numMyBlockRows = NumMyBlockRows();
01793 
01794   if(IndicesAreGlobal()) {
01795     // Check if ColMap is monotone. If not, the list will get unsorted.
01796     bool mapMonotone = true;
01797     {
01798       long long oldGID = colmap.GID64(0);
01799       for (int i=1; i<colmap.NumMyElements(); ++i) {
01800         if (oldGID > colmap.GID64(i)) {
01801           mapMonotone = false;
01802           break;
01803         }
01804         oldGID = colmap.GID64(i);
01805       }
01806     }
01807     if (Sorted())
01808       SetSorted(mapMonotone);
01809 
01810   // We don't call Data<int>() here because that would not work (it will throw)
01811   // if GlobalIndicesLongLong() and IndicesAreGlobal().  This is because
01812   // the local flag is not set yet.  We are in the middle of the transaction here.
01813   // In all other cases, one should call the function Data<int> or Data<long long>
01814   // instead of obtaining the pointer directly.
01815   Epetra_CrsGraphData::IndexData<int>& intData = *(CrsGraphData_->data);
01816 
01817   if(RowMap().GlobalIndicesInt())
01818   {
01819     // now comes the actual transformation
01820     for(int i = 0; i < numMyBlockRows; i++) {
01821       const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01822       int* ColIndices = intData.Indices_[i];
01823       for(int j = 0; j < NumIndices; j++) {
01824       int GID = ColIndices[j];
01825       int LID = colmap.LID(GID);
01826       if(LID != -1)
01827         ColIndices[j] = LID;
01828       else
01829         throw ReportError("Internal error in FillComplete ",-1);
01830       }
01831     }
01832   }
01833   else if(RowMap().GlobalIndicesLongLong())
01834   {
01835 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01836     Epetra_CrsGraphData::IndexData<long long>& LL_Data = CrsGraphData_->Data<long long>();
01837 
01838     if (!CrsGraphData_->StaticProfile_) {
01839       // Use the resize trick used in TAllocate.
01840       const long long indexBaseMinusOne = IndexBase64() - 1;
01841       for(int i = 0; i < numMyBlockRows; i++) {
01842         const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01843         intData.SortedEntries_[i].entries_.resize(NumIndices, indexBaseMinusOne);
01844         intData.Indices_[i] = NumIndices > 0 ? &intData.SortedEntries_[i].entries_[0]: NULL;
01845         intData.SortedEntries_[i].entries_.resize(0);
01846       }
01847     }
01848 
01849     // now comes the actual transformation
01850     for(int i = 0; i < numMyBlockRows; i++) {
01851       const int NumIndices = CrsGraphData_->NumIndicesPerRow_[i];
01852       long long* ColIndices = LL_Data.Indices_[i];
01853       int* intColIndices = intData.Indices_[i];
01854       for(int j = 0; j < NumIndices; j++) {
01855       long long GID = ColIndices[j];
01856       int LID = colmap.LID(GID);
01857       if(LID != -1)
01858         intColIndices[j] = LID;
01859       else
01860         throw ReportError("Internal error in FillComplete ",-1);
01861       }
01862     }
01863 
01864     LL_Data.Deallocate(); // deallocate long long data since indices are local now.
01865 #else
01866   throw ReportError("Epetra_CrsGraph::MakeIndicesLocal: GlobalIndicesLongLong but no long long API", -1);
01867 #endif
01868   }
01869 
01870   }
01871 
01872   SetIndicesAreLocal(true);
01873   SetIndicesAreGlobal(false);
01874 
01875   if(CrsGraphData_->ReferenceCount() > 1)
01876     return(1); // return 1 if data is shared
01877   else
01878     return(0);
01879 }
01880 
01881 //==============================================================================
01882 int Epetra_CrsGraph::OptimizeStorage() {
01883   int NumIndices;
01884   const int numMyBlockRows = NumMyBlockRows();
01885 
01886   Epetra_CrsGraphData::IndexData<int>& Data = CrsGraphData_->Data<int>();
01887 
01888   if(StorageOptimized())
01889     return(0); // Have we been here before?
01890   if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
01891 
01892   bool Contiguous = true; // Assume contiguous is true
01893   for(int i = 1; i < numMyBlockRows; i++) {
01894     NumIndices = CrsGraphData_->NumIndicesPerRow_[i-1];
01895     int NumAllocateIndices = CrsGraphData_->NumAllocatedIndicesPerRow_[i-1];
01896 
01897     // Check if NumIndices is same as NumAllocatedIndices and
01898     // check if end of beginning of current row starts immediately after end of previous row.
01899     if((NumIndices != NumAllocateIndices) ||
01900        (Data.Indices_[i] != Data.Indices_[i-1] + NumIndices)) {
01901       Contiguous = false;
01902       break;
01903     }
01904   }
01905 
01906   // NOTE:  At the end of the above loop set, there is a possibility that NumIndices and NumAllocatedIndices
01907   //        for the last row could be different, but I don't think it matters.
01908 
01909 
01910   if((CrsGraphData_->CV_ == View) && !Contiguous)
01911     return(3);  // This is user data, it's not contiguous and we can't make it so.
01912 
01913   // This next step constructs the scan sum of the number of indices per row.  Note that the
01914   // storage associated with NumIndicesPerRow is used by IndexOffset, so we have to be
01915   // careful with this sum operation
01916 
01917   if(CrsGraphData_->IndexOffset_.Values() != CrsGraphData_->NumIndicesPerRow_.Values())
01918     CrsGraphData_->IndexOffset_.MakeViewOf(CrsGraphData_->NumIndicesPerRow_);
01919 
01920   int * numIndicesPerRow = CrsGraphData_->NumIndicesPerRow_.Values();
01921   int curNumIndices = numIndicesPerRow[0];
01922   numIndicesPerRow[0] = 0;
01923   for (int i=0; i<numMyBlockRows; ++i) {
01924     int nextNumIndices = numIndicesPerRow[i+1];
01925     numIndicesPerRow[i+1] = numIndicesPerRow[i]+curNumIndices;
01926     curNumIndices = nextNumIndices;
01927   }
01928 
01929 // *******************************
01930 // Data NOT contiguous, make it so
01931 // *******************************
01932   if(!Contiguous) { // Must pack indices if not already contiguous
01933 
01934     // Allocate one big integer array for all index values
01935     if (!(StaticProfile())) { // If static profile, All_Indices_ is already allocated, only need to pack data
01936       int errorcode = Data.All_Indices_.Size(CrsGraphData_->NumMyNonzeros_);
01937       if(errorcode != 0) throw ReportError("Error with All_Indices_ allocation.", -99);
01938     }
01939       // Pack indices into All_Indices_
01940 
01941     int* all_indices = Data.All_Indices_.Values();
01942     int * indexOffset = CrsGraphData_->IndexOffset_.Values();
01943     int ** indices = Data.Indices_;
01944 
01945     if (!(StaticProfile())) {
01946 #ifdef EPETRA_HAVE_OMP
01947 #pragma omp parallel for default(none) shared(indexOffset,all_indices,indices)
01948 #endif
01949       for(int i = 0; i < numMyBlockRows; i++) {
01950   int numColIndices = indexOffset[i+1] - indexOffset[i];
01951         int* ColIndices = indices[i];
01952         int *newColIndices = all_indices+indexOffset[i];
01953         for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
01954       }
01955       for(int i = 0; i < numMyBlockRows; i++) {
01956         if (indices[i]!=0) {
01957           Data.SortedEntries_[i].entries_.clear();
01958           indices[i] = 0;
01959         }
01960      }
01961    } // End of non-contiguous non-static section
01962    else {
01963 
01964      for(int i = 0; i < numMyBlockRows; i++) {
01965        int numColIndices = indexOffset[i+1] - indexOffset[i];
01966        int* ColIndices = indices[i];
01967        int *newColIndices = all_indices+indexOffset[i];
01968        if (ColIndices!=newColIndices) // No need to copy if pointing to same space
01969        for(int j = 0; j < numColIndices; j++) newColIndices[j] = ColIndices[j];
01970        indices[i] = 0;
01971      }
01972    } // End of non-Contiguous static section
01973   } // End of non-Contiguous section
01974   else { // Start of Contiguous section
01975     // if contiguous, set All_Indices_ from CrsGraphData_->Indices_[0].
01976     // Execute the assignment block in parallel using the same pattern as SpMV
01977     // in order to improve page placement
01978     if (numMyBlockRows > 0 && !(StaticProfile())) {
01979       const int numMyNonzeros = NumMyNonzeros();
01980       int errorcode = Data.All_Indices_.Size(numMyNonzeros);
01981       if(errorcode != 0)  throw ReportError("Error with All_Indices_ allocation.", -99);
01982       int* new_all_indices = Data.All_Indices_.Values();
01983       int* old_all_indices = Data.Indices_[0];
01984       int * indexOffset = CrsGraphData_->IndexOffset_.Values();
01985 
01986 #ifdef EPETRA_HAVE_OMP
01987 #pragma omp parallel for default(none) shared(indexOffset,old_all_indices,new_all_indices)
01988 #endif
01989      for(int i = 0; i < numMyBlockRows; i++) {
01990        int numColIndices = indexOffset[i+1] - indexOffset[i];
01991        int *oldColIndices = old_all_indices+indexOffset[i];
01992        int *newColIndices = new_all_indices+indexOffset[i];
01993        for(int j = 0; j < numColIndices; j++) newColIndices[j] = oldColIndices[j];
01994      }
01995 
01996 //#ifdef EPETRA_HAVE_OMP
01997 //#pragma omp parallel for default(none) shared(all_indices_values,indices_values)
01998 //#endif
01999 //      for(int ii=0; ii<numMyNonzeros; ++ii) {
02000 //        all_indices_values[ii] = indices_values[ii];
02001 //      }
02002     }
02003   }
02004 
02005   // Delete unneeded storage
02006   CrsGraphData_->NumAllocatedIndicesPerRow_.Resize(0);
02007   delete [] Data.Indices_; Data.Indices_=0;
02008   Data.SortedEntries_.clear();
02009 
02010   SetIndicesAreContiguous(true); // Can no longer dynamically add or remove indices
02011   CrsGraphData_->StorageOptimized_ = true;
02012 
02013 /*
02014 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
02015   All_IndicesPlus1(); // see if preemptively calling this improves Multiply timing.
02016 #endif
02017 */
02018 
02019   return(0);
02020 }
02021 
02022 //==============================================================================
02023 template<typename int_type>
02024 int Epetra_CrsGraph::ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
02025 {
02026   int j;
02027 
02028   int locRow = LRID(Row); // Normalize row range
02029 
02030   if(locRow < 0 || locRow >= NumMyBlockRows())
02031     EPETRA_CHK_ERR(-1); // Not in Row range
02032 
02033   NumIndices = NumMyIndices(locRow);
02034   if(LenOfIndices < NumIndices)
02035     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
02036 
02037   if(IndicesAreLocal())
02038   {
02039     int * srcIndices = TIndices<int>(locRow);
02040     // static_cast is ok because global indices were created from int values and hence must fit ints
02041     for(j = 0; j < NumIndices; j++)
02042       targIndices[j] = static_cast<int_type>(GCID64(srcIndices[j]));
02043   }
02044   else
02045   {
02046     int_type * srcIndices = TIndices<int_type>(locRow);
02047     for(j = 0; j < NumIndices; j++)
02048       targIndices[j] = srcIndices[j];
02049   }
02050 
02051   return(0);
02052 }
02053 
02054 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02055 int Epetra_CrsGraph::ExtractGlobalRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
02056 {
02057   if(RowMap().GlobalIndicesInt())
02058     return ExtractGlobalRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
02059   else
02060   throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy int version called for a graph that is not int.", -1);
02061 }
02062 #endif
02063 
02064 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02065 int Epetra_CrsGraph::ExtractGlobalRowCopy(long long Row, int LenOfIndices, int& NumIndices, long long* targIndices) const
02066 {
02067   if(RowMap().GlobalIndicesLongLong())
02068     return ExtractGlobalRowCopy<long long>(Row, LenOfIndices, NumIndices, targIndices);
02069   else
02070   throw ReportError("Epetra_CrsGraph::ExtractGlobalRowCopy long long version called for a graph that is not long long.", -1);
02071 }
02072 #endif
02073 
02074 //==============================================================================
02075 template<typename int_type>
02076 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int_type* targIndices) const
02077 {
02078   int j;
02079 
02080   if(Row < 0 || Row >= NumMyBlockRows())
02081     EPETRA_CHK_ERR(-1); // Not in Row range
02082 
02083   NumIndices = NumMyIndices(Row);
02084   if(LenOfIndices < NumIndices)
02085     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumIndices
02086 
02087   if(IndicesAreGlobal())
02088     EPETRA_CHK_ERR(-3); // There are no local indices yet
02089 
02090   int * srcIndices = TIndices<int>(Row);
02091   for(j = 0; j < NumIndices; j++)
02092     targIndices[j] = srcIndices[j];
02093 
02094   return(0);
02095 }
02096 
02097 int Epetra_CrsGraph::ExtractMyRowCopy(int Row, int LenOfIndices, int& NumIndices, int* targIndices) const
02098 {
02099   if(RowMap().GlobalIndicesTypeValid())
02100     return ExtractMyRowCopy<int>(Row, LenOfIndices, NumIndices, targIndices);
02101   else
02102     throw ReportError("Epetra_CrsGraph::ExtractMyRowCopy graph global index type unknown.", -1);
02103 }
02104 
02105 //==============================================================================
02106 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02107 int Epetra_CrsGraph::ExtractGlobalRowView(int Row, int& NumIndices, int*& targIndices) const
02108 {
02109   if(!RowMap().GlobalIndicesInt())
02110     throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView int version called for a graph that is not int.", -1);
02111 
02112   int locRow = LRID(Row); // Normalize row range
02113 
02114   if(locRow < 0 || locRow >= NumMyBlockRows())
02115     EPETRA_CHK_ERR(-1); // Not in Row range
02116 
02117   if(IndicesAreLocal())
02118     EPETRA_CHK_ERR(-2); // There are no global indices
02119 
02120   NumIndices = NumMyIndices(locRow);
02121 
02122   targIndices = TIndices<int>(locRow);
02123 
02124   return(0);
02125 }
02126 #endif
02127 //==============================================================================
02128 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02129 int Epetra_CrsGraph::ExtractGlobalRowView(long long Row, int& NumIndices, long long*& targIndices) const
02130 {
02131   if(!RowMap().GlobalIndicesLongLong())
02132     throw ReportError("Epetra_CrsGraph::ExtractGlobalRowView long long version called for a graph that is not long long.", -1);
02133 
02134   int locRow = LRID(Row); // Normalize row range
02135 
02136   if(locRow < 0 || locRow >= NumMyBlockRows())
02137     EPETRA_CHK_ERR(-1); // Not in Row range
02138 
02139   if(IndicesAreLocal())
02140     EPETRA_CHK_ERR(-2); // There are no global indices
02141 
02142   NumIndices = NumMyIndices(locRow);
02143 
02144   targIndices = TIndices<long long>(locRow);
02145 
02146   return(0);
02147 }
02148 #endif
02149 //==============================================================================
02150 int Epetra_CrsGraph::ExtractMyRowView(int Row, int& NumIndices, int*& targIndices) const
02151 {
02152   if(Row < 0 || Row >= NumMyBlockRows())
02153     EPETRA_CHK_ERR(-1); // Not in Row range
02154 
02155   if(IndicesAreGlobal())
02156     EPETRA_CHK_ERR(-2); // There are no local indices
02157 
02158   NumIndices = NumMyIndices(Row);
02159 
02160   targIndices = TIndices<int>(Row);
02161 
02162   return(0);
02163 }
02164 
02165 //==============================================================================
02166 int Epetra_CrsGraph::NumGlobalIndices(long long Row) const {
02167 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
02168   int locRow = LRID((int) Row);
02169 #else
02170   int locRow = LRID(Row);
02171 #endif
02172   if(locRow != -1)
02173     return(NumMyIndices(locRow));
02174   else
02175     return(0); // No indices for this row on this processor
02176 }
02177 
02178 //==============================================================================
02179 int Epetra_CrsGraph::NumAllocatedGlobalIndices(long long Row) const {
02180 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
02181   int locRow = LRID((int) Row);
02182 #else
02183   int locRow = LRID(Row);
02184 #endif
02185   if(locRow != -1)
02186     return(NumAllocatedMyIndices(locRow));
02187   else
02188     return(0); // No indices allocated for this row on this processor
02189 }
02190 
02191 //==============================================================================
02192 int Epetra_CrsGraph::ReplaceRowMap(const Epetra_BlockMap& newmap)
02193 {
02194   if (RowMap().PointSameAs(newmap)) {
02195     Epetra_DistObject::Map_ = newmap;
02196     CrsGraphData_->RowMap_ = newmap;
02197     CrsGraphData_->MakeImportExport();
02198     return(0);
02199   }
02200 
02201   return(-1);
02202 }
02203 
02204 //==============================================================================
02205 int Epetra_CrsGraph::ReplaceColMap(const Epetra_BlockMap& newmap)
02206 {
02207   if (!HaveColMap() && !IndicesAreLocal() && !IndicesAreGlobal() && newmap.GlobalIndicesTypeMatch(RowMap())) {
02208     CrsGraphData_->ColMap_            = newmap;
02209     CrsGraphData_->NumGlobalBlockCols_= newmap.NumGlobalElements64();
02210     CrsGraphData_->NumMyBlockCols_    = newmap.NumMyElements();
02211     CrsGraphData_->MaxColDim_         = newmap.MaxElementSize();
02212     CrsGraphData_->GlobalMaxColDim_   = newmap.MaxElementSize();
02213     CrsGraphData_->NumGlobalCols_     = newmap.NumGlobalPoints64();
02214     CrsGraphData_->NumMyCols_         = newmap.NumMyPoints();
02215     CrsGraphData_->HaveColMap_        = true;
02216     return(0);
02217   }
02218 
02219   if(ColMap().PointSameAs(newmap)) {
02220     CrsGraphData_->ColMap_ = newmap;
02221     CrsGraphData_->MakeImportExport();
02222     return(0);
02223   }
02224 
02225   return(-1);
02226 }
02227 
02228 
02229 //==============================================================================
02230 int Epetra_CrsGraph::ReplaceDomainMapAndImporter(const Epetra_BlockMap& NewDomainMap, const Epetra_Import * NewImporter) {
02231   int rv=0;
02232   if( !NewImporter && ColMap().SameAs(NewDomainMap)) {
02233     CrsGraphData_->DomainMap_ = NewDomainMap;
02234     delete CrsGraphData_->Importer_;
02235     CrsGraphData_->Importer_ = 0;
02236   }
02237   else if(NewImporter && ColMap().SameAs(NewImporter->TargetMap()) && NewDomainMap.SameAs(NewImporter->SourceMap())) {
02238     CrsGraphData_->DomainMap_ = NewDomainMap;
02239     delete CrsGraphData_->Importer_;
02240     CrsGraphData_->Importer_  = new Epetra_Import(*NewImporter);
02241   }
02242   else
02243     rv=-1;
02244   return rv;
02245 }
02246 
02247 //==============================================================================
02248 int Epetra_CrsGraph::RemoveEmptyProcessesInPlace(const Epetra_BlockMap * newMap) {
02249   const Epetra_BlockMap *newDomainMap=0, *newRangeMap=0, *newColMap=0;
02250   Epetra_Import * newImport=0;
02251   Epetra_Export * newExport=0;
02252 
02253   const Epetra_Comm *newComm = newMap ? &newMap->Comm() : 0;
02254 
02255   if(DomainMap().SameAs(RowMap())) {
02256     // Common case: original domain and row Maps are identical.
02257     // In that case, we need only replace the original domain Map
02258     // with the new Map.  This ensures that the new domain and row
02259     // Maps _stay_ identical.
02260     newDomainMap = newMap;
02261   }
02262   else
02263     newDomainMap = DomainMap().ReplaceCommWithSubset(newComm);
02264 
02265   if(RangeMap().SameAs(RowMap())){
02266     // Common case: original range and row Maps are identical.  In
02267     // that case, we need only replace the original range Map with
02268     // the new Map.  This ensures that the new range and row Maps
02269     // _stay_ identical.
02270     newRangeMap = newMap;
02271   }
02272   else
02273     newRangeMap = RangeMap().ReplaceCommWithSubset(newComm);
02274 
02275   newColMap=ColMap().ReplaceCommWithSubset(newComm);
02276 
02277   if(newComm) {
02278     // (Re)create the Export and / or Import if necessary.
02279     //
02280     // The operations below are collective on the new communicator.
02281     //
02282     if(RangeMap().DataPtr() != RowMap().DataPtr())
02283       newExport = new Epetra_Export(*newMap,*newRangeMap);
02284 
02285     if(DomainMap().DataPtr() != ColMap().DataPtr())
02286       newImport = new Epetra_Import(*newColMap,*newDomainMap);
02287   }
02288 
02289   // CrsGraphData things
02290   if(CrsGraphData_->ReferenceCount() !=1)
02291     throw ReportError("Epetra_CrsGraph::RemoveEmptyProcessesInPlace does not work for shared CrsGraphData_",-2);
02292 
02293   // Dummy map for the non-active procs
02294 #if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
02295   Epetra_Map dummy;
02296 #else
02297   Epetra_SerialComm SComm;
02298   Epetra_Map dummy(0,0,SComm);
02299 #endif
02300 
02301   delete CrsGraphData_->Importer_;
02302   delete CrsGraphData_->Exporter_;
02303 
02304 
02305   CrsGraphData_->RowMap_    = newMap       ? *newMap      : dummy;
02306   CrsGraphData_->ColMap_    = newColMap    ? *newColMap    : dummy;
02307   CrsGraphData_->DomainMap_ = newDomainMap ? *newDomainMap : dummy;
02308   CrsGraphData_->RangeMap_  = newRangeMap  ? *newRangeMap : dummy;
02309   CrsGraphData_->Importer_  = newImport;
02310   CrsGraphData_->Exporter_  = newExport;
02311 
02312   // Epetra_DistObject things
02313   if(newMap) {
02314     Map_  = *newMap;
02315     Comm_ = &newMap->Comm();
02316   }
02317 
02318   // Cleanup (newRowMap is always newMap, so don't delete that)
02319   if(newColMap != newMap)    delete newColMap;
02320   if(newDomainMap != newMap) delete newDomainMap;
02321   if(newRangeMap != newMap)  delete newRangeMap;
02322 
02323   return(0);
02324 }
02325 
02326 // private =====================================================================
02327 int Epetra_CrsGraph::CheckSizes(const Epetra_SrcDistObject& Source) {
02328   try {
02329     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source); // downcast Source from SrcDistObject to CrsGraph
02330     if(!A.GlobalConstantsComputed())
02331       EPETRA_CHK_ERR(-1); // Must have global constants to proceed
02332   }
02333   catch(...) {
02334     return(0); // No error at this point, object could be a RowMatrix
02335   }
02336   return(0);
02337 }
02338 
02339 // private =====================================================================
02340 int Epetra_CrsGraph::CopyAndPermute(const Epetra_SrcDistObject& Source,
02341             int NumSameIDs,
02342             int NumPermuteIDs,
02343             int* PermuteToLIDs,
02344             int* PermuteFromLIDs,
02345                                     const Epetra_OffsetIndex * Indexor)
02346 {
02347   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02348     throw ReportError("Epetra_CrsGraph::CopyAndPermute: Incoming global index type does not match the one for *this",-1);
02349 
02350   try {
02351     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
02352     EPETRA_CHK_ERR(CopyAndPermuteCrsGraph(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
02353             PermuteFromLIDs,Indexor));
02354   }
02355   catch(...) {
02356     try {
02357       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
02358       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
02359                PermuteFromLIDs,Indexor));
02360     }
02361     catch(...) {
02362       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
02363     }
02364   }
02365 
02366   return(0);
02367 }
02368 
02369 // private =====================================================================
02370 template<typename int_type>
02371 int Epetra_CrsGraph::CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
02372                int NumSameIDs,
02373                int NumPermuteIDs,
02374                int* PermuteToLIDs,
02375                int* PermuteFromLIDs,
02376                                              const Epetra_OffsetIndex * Indexor)
02377 {
02378   (void)Indexor;
02379   int i;
02380   int j;
02381   int NumIndices;
02382   int FromRow;
02383   int_type ToRow;
02384   int maxNumIndices = A.MaxNumEntries();
02385   Epetra_IntSerialDenseVector local_indices_vec;
02386 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02387   Epetra_LongLongSerialDenseVector global_indices_vec;
02388 #endif
02389   Epetra_SerialDenseVector Values;
02390 
02391   int* local_indices = 0;
02392   int_type* global_indices = 0;
02393 
02394   if(maxNumIndices > 0) {
02395     local_indices_vec.Size(maxNumIndices);
02396   local_indices = local_indices_vec.Values();
02397 
02398   if(A.RowMatrixRowMap().GlobalIndicesLongLong())
02399   {
02400 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02401     global_indices_vec.Size(maxNumIndices);
02402     global_indices = reinterpret_cast<int_type*>(global_indices_vec.Values());
02403 #else
02404     throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: GlobalIndicesLongLong but no API for long long",-1);
02405 #endif
02406   }
02407   else
02408   {
02409     global_indices = reinterpret_cast<int_type*>(local_indices);
02410   }
02411 
02412     Values.Size(maxNumIndices); // Must extract values even though we discard them
02413   }
02414 
02415   const Epetra_Map& rowMap = A.RowMatrixRowMap();
02416   const Epetra_Map& colMap = A.RowMatrixColMap();
02417 
02418   // Do copy first
02419   for(i = 0; i < NumSameIDs; i++) {
02420     ToRow = (int) rowMap.GID64(i);
02421     EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumIndices, NumIndices, Values.Values(), local_indices));
02422     for(j = 0; j < NumIndices; j++)
02423       global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
02424     // Place into target graph.
02425     int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices);
02426     if(ierr < 0) EPETRA_CHK_ERR(ierr);
02427   }
02428 
02429   // Do local permutation next
02430   for(i = 0; i < NumPermuteIDs; i++) {
02431     FromRow = PermuteFromLIDs[i];
02432     ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02433     EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumIndices, NumIndices, Values.Values(), local_indices));
02434     for(j = 0; j < NumIndices; j++)
02435       global_indices[j] = (int_type) colMap.GID64(local_indices[j]); // convert to GIDs
02436     int ierr = InsertGlobalIndices(ToRow, NumIndices, global_indices); // Place into target graph.
02437     if(ierr < 0) EPETRA_CHK_ERR(ierr);
02438   }
02439 
02440   return(0);
02441 }
02442 
02443 int Epetra_CrsGraph::CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
02444                int NumSameIDs,
02445                int NumPermuteIDs,
02446                int* PermuteToLIDs,
02447                int* PermuteFromLIDs,
02448                                              const Epetra_OffsetIndex * Indexor)
02449 {
02450   if(!A.RowMatrixRowMap().GlobalIndicesTypeMatch(RowMap()))
02451     throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
02452 
02453 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02454   if(A.RowMatrixRowMap().GlobalIndicesInt())
02455     return CopyAndPermuteRowMatrix<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02456   else
02457 #endif
02458 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02459   if(A.RowMatrixRowMap().GlobalIndicesLongLong())
02460     return CopyAndPermuteRowMatrix<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02461   else
02462 #endif
02463   throw ReportError("Epetra_CrsGraph::CopyAndPermuteRowMatrix: Unable to determine global index type of map", -1);
02464 }
02465 
02466 // private =====================================================================
02467 template<typename int_type>
02468 int Epetra_CrsGraph::CopyAndPermuteCrsGraph(const Epetra_CrsGraph& A,
02469               int NumSameIDs,
02470               int NumPermuteIDs,
02471               int* PermuteToLIDs,
02472               int* PermuteFromLIDs,
02473                                             const Epetra_OffsetIndex * Indexor)
02474 {
02475   (void)Indexor;
02476   int i;
02477   int_type Row;
02478   int NumIndices;
02479   int_type* indices = 0;
02480   int_type FromRow, ToRow;
02481   int maxNumIndices = A.MaxNumIndices();
02482   Epetra_IntSerialDenseVector int_IndicesVector;
02483 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02484   Epetra_LongLongSerialDenseVector LL_IndicesVector;
02485 #endif
02486 
02487   if(maxNumIndices > 0 && A.IndicesAreLocal()) {
02488     if(A.RowMap().GlobalIndicesInt())
02489       {
02490         int_IndicesVector.Size(maxNumIndices);
02491         indices = reinterpret_cast<int_type*>(int_IndicesVector.Values());
02492       }
02493     else if(A.RowMap().GlobalIndicesLongLong())
02494       {
02495 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02496   LL_IndicesVector.Size(maxNumIndices);
02497   indices = reinterpret_cast<int_type*>(LL_IndicesVector.Values());
02498 #else
02499     throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02500 #endif
02501       }
02502   }
02503 
02504   // Do copy first
02505   if(NumSameIDs > 0) {
02506     if(A.IndicesAreLocal()) {
02507       for(i = 0; i < NumSameIDs; i++) {
02508         Row = (int_type) GRID64(i);
02509         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumIndices, NumIndices, indices));
02510         // Place into target graph.
02511         int ierr = InsertGlobalIndices(Row, NumIndices, indices);
02512         if(ierr < 0) EPETRA_CHK_ERR(ierr);
02513       }
02514     }
02515     else { // A.IndiceAreGlobal()
02516       for(i = 0; i < NumSameIDs; i++) {
02517         Row = (int_type) GRID64(i);
02518         EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumIndices, indices));
02519         // Place into target graph.
02520         int ierr = InsertGlobalIndices(Row, NumIndices, indices);
02521         if(ierr < 0) EPETRA_CHK_ERR(ierr);
02522       }
02523     }
02524   }
02525 
02526   // Do local permutation next
02527   if(NumPermuteIDs > 0) {
02528     if(A.IndicesAreLocal()) {
02529       for(i = 0; i < NumPermuteIDs; i++) {
02530         FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02531         ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02532         EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
02533         // Place into target graph.
02534         int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
02535         if (ierr < 0) EPETRA_CHK_ERR(ierr);
02536       }
02537     }
02538     else { // A.IndiceAreGlobal()
02539       for(i = 0; i < NumPermuteIDs; i++) {
02540         FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02541         ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02542         EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumIndices, indices));
02543         // Place into target graph.
02544         int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
02545         if (ierr < 0) EPETRA_CHK_ERR(ierr);
02546       }
02547     }
02548   }
02549 
02550   return(0);
02551 }
02552 
02553 int Epetra_CrsGraph::CopyAndPermuteCrsGraph(const Epetra_CrsGraph& A,
02554               int NumSameIDs,
02555               int NumPermuteIDs,
02556               int* PermuteToLIDs,
02557               int* PermuteFromLIDs,
02558                                             const Epetra_OffsetIndex * Indexor)
02559 {
02560   if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
02561     throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Incoming global index type does not match the one for *this",-1);
02562 
02563   if(A.RowMap().GlobalIndicesInt())
02564 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02565     return CopyAndPermuteCrsGraph<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02566 #else
02567     throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesInt but no API for it.",-1);
02568 #endif
02569 
02570   if(A.RowMap().GlobalIndicesLongLong())
02571 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02572     return CopyAndPermuteCrsGraph<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02573 #else
02574     throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02575 #endif
02576 
02577   throw ReportError("Epetra_CrsGraph::CopyAndPermuteCrsGraph: Unable to determine global index type of map", -1);
02578 }
02579 
02580 // private =====================================================================
02581 int Epetra_CrsGraph::PackAndPrepare(const Epetra_SrcDistObject& Source,
02582             int NumExportIDs,
02583             int* ExportLIDs,
02584             int& LenExports,
02585             char*& Exports,
02586             int& SizeOfPacket,
02587                                     int * Sizes,
02588                                     bool& VarSizes,
02589             Epetra_Distributor& Distor)
02590 {
02591   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02592     throw ReportError("Epetra_CrsGraph::PackAndPrepare: Incoming global index type does not match the one for *this",-1);
02593 
02594   int globalMaxNumIndices = 0;
02595   int TotalSendSize = 0;
02596 
02597   VarSizes = true;
02598 
02599   if(Source.Map().GlobalIndicesInt())
02600     SizeOfPacket = (int)sizeof(int);
02601   else if(Source.Map().GlobalIndicesLongLong())
02602     SizeOfPacket = (int)sizeof(long long);
02603   else
02604     throw ReportError("Epetra_CrsGraph::PackAndPrepare: Unable to determine source global index type",-1);
02605 
02606   if(NumExportIDs <= 0) return(0);
02607 
02608   try {
02609     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
02610     globalMaxNumIndices = A.GlobalMaxNumIndices();
02611     for( int i = 0; i < NumExportIDs; ++i )
02612     {
02613       Sizes[i] = (A.NumMyIndices( ExportLIDs[i] ) + 2);
02614       TotalSendSize += Sizes[i];
02615     }
02616   }
02617   catch(...) {
02618     try {
02619       const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
02620       int maxNumIndices = A.MaxNumEntries();
02621       A.Comm().MaxAll(&maxNumIndices, &globalMaxNumIndices, 1);
02622       for( int i = 0; i < NumExportIDs; ++i )
02623       {
02624         int NumEntries;
02625         A.NumMyRowEntries( ExportLIDs[i], NumEntries );
02626         Sizes[i] = (NumEntries + 2);
02627         TotalSendSize += Sizes[i];
02628       }
02629     }
02630     catch(...) {
02631       EPETRA_CHK_ERR(-1); // Bad cast
02632     }
02633   }
02634 
02635   CrsGraphData_->ReAllocateAndCast(Exports, LenExports, TotalSendSize*SizeOfPacket);
02636 
02637   try {
02638     const Epetra_CrsGraph& A = dynamic_cast<const Epetra_CrsGraph&>(Source);
02639     EPETRA_CHK_ERR(PackAndPrepareCrsGraph(A, NumExportIDs, ExportLIDs, LenExports, Exports,
02640             SizeOfPacket, Sizes, VarSizes, Distor));
02641   }
02642   catch(...) {
02643     const Epetra_RowMatrix& A = dynamic_cast<const Epetra_RowMatrix&>(Source);
02644     EPETRA_CHK_ERR(PackAndPrepareRowMatrix(A, NumExportIDs, ExportLIDs, LenExports, Exports,
02645              SizeOfPacket, Sizes, VarSizes, Distor));
02646   }
02647   return(0);
02648 }
02649 
02650 // private =====================================================================
02651 int Epetra_CrsGraph::PackAndPrepareCrsGraph(const Epetra_CrsGraph& A,
02652               int NumExportIDs,
02653               int* ExportLIDs,
02654               int& LenExports,
02655               char*& Exports,
02656               int& SizeOfPacket,
02657                                             int* Sizes,
02658                                             bool& VarSizes,
02659               Epetra_Distributor& Distor)
02660 {
02661   if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
02662     throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Incoming global index type does not match the one for *this",-1);
02663 
02664   (void)LenExports;
02665   (void)SizeOfPacket;
02666   (void)Sizes;
02667   (void)VarSizes;
02668   (void)Distor;
02669   int i;
02670   int NumIndices;
02671 
02672   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02673   // 1st int: GRID of row where GRID is the global row ID for the source graph
02674   // next int:  NumIndices, Number of indices in row.
02675   // next NumIndices: The actual indices for the row.
02676   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
02677   //   sized segments for current communication routines.
02678   int maxNumIndices = A.MaxNumIndices();
02679   //if( maxNumIndices ) indices = new int[maxNumIndices];
02680 
02681   if(A.RowMap().GlobalIndicesInt()) {
02682     int* indices = 0;
02683     int* intptr = (int*) Exports;
02684     int FromRow;
02685     for(i = 0; i < NumExportIDs; i++) {
02686       FromRow = (int) A.GRID64(ExportLIDs[i]);
02687       *intptr = FromRow;
02688       indices = intptr + 2;
02689       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
02690       intptr[1] = NumIndices; // Load second slot of segment
02691       intptr += (NumIndices+2); // Point to next segment
02692     }
02693   }
02694   else if(A.RowMap().GlobalIndicesLongLong()) {
02695 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02696     long long* indices = 0;
02697     long long* LLptr = (long long*) Exports;
02698     long long FromRow;
02699     for(i = 0; i < NumExportIDs; i++) {
02700       FromRow = A.GRID64(ExportLIDs[i]);
02701       *LLptr = FromRow;
02702       indices = LLptr + 2;
02703       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumIndices, NumIndices, indices));
02704       LLptr[1] = NumIndices; // Load second slot of segment
02705       LLptr += (NumIndices+2); // Point to next segment
02706     }
02707 #else
02708     throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02709 #endif
02710   }
02711   else
02712     throw ReportError("Epetra_CrsGraph::PackAndPrepareCrsGraph: Unable to determine source global index type",-1);
02713 
02714   //if( indices ) delete [] indices;
02715 
02716   return(0);
02717 }
02718 
02719 // private =====================================================================
02720 int Epetra_CrsGraph::PackAndPrepareRowMatrix(const Epetra_RowMatrix& A,
02721                int NumExportIDs,
02722                int* ExportLIDs,
02723                int& LenExports,
02724                char*& Exports,
02725                int& SizeOfPacket,
02726                                              int* Sizes,
02727                                              bool& VarSizes,
02728                Epetra_Distributor& Distor)
02729 {
02730   if(!A.Map().GlobalIndicesTypeMatch(RowMap()))
02731     throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Incoming global index type does not match the one for *this",-1);
02732 
02733   (void)LenExports;
02734   (void)SizeOfPacket;
02735   (void)Sizes;
02736   (void)VarSizes;
02737   (void)Distor;
02738   int i;
02739   int j;
02740   int NumIndices;
02741   Epetra_SerialDenseVector Values;
02742 
02743   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02744   // 1st int: GRID of row where GRID is the global row ID for the source graph
02745   // next int:  NumIndices, Number of indices in row.
02746   // next NumIndices: The actual indices for the row.
02747   // Any remaining space (of length GlobalMaxNumIndices - NumIndices ints) will be wasted but we need fixed
02748   //   sized segments for current communication routines.
02749   int maxNumIndices = A.MaxNumEntries();
02750   if(maxNumIndices > 0) {
02751     Values.Size(maxNumIndices);
02752 //    indices = new int[maxNumIndices];
02753   }
02754   const Epetra_Map& rowMap = A.RowMatrixRowMap();
02755   const Epetra_Map& colMap = A.RowMatrixColMap();
02756 
02757   if(rowMap.GlobalIndicesInt() && colMap.GlobalIndicesInt()) {
02758     int* indices = 0;
02759     int FromRow;
02760     int* intptr = (int*) Exports;
02761     for(i = 0; i < NumExportIDs; i++) {
02762       FromRow = (int) rowMap.GID64(ExportLIDs[i]);
02763       *intptr = FromRow;
02764       indices = intptr + 2;
02765       EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), indices));
02766       for(j = 0; j < NumIndices; j++) indices[j] = (int) colMap.GID64(indices[j]); // convert to GIDs
02767       intptr[1] = NumIndices; // Load second slot of segment
02768       intptr += (NumIndices+2); // Point to next segment
02769     }
02770   }
02771   else if(rowMap.GlobalIndicesLongLong() && colMap.GlobalIndicesLongLong()) {
02772     // Bytes of Exports:
02773     // 12345678.12345678....12345678.12345678 ("." means no spaces)
02774     // FromRow  NumIndices  id1 id2  id3 id4  <-- before converting to GIDs
02775     // FromRow  NumIndices  | gid1 | | gid2 | <-- after converting to GIDs
02776 
02777     long long* LL_indices = 0;
02778     long long FromRow;
02779     long long* LLptr = (long long*) Exports;
02780     for(i = 0; i < NumExportIDs; i++) {
02781       FromRow = rowMap.GID64(ExportLIDs[i]);
02782       *LLptr = FromRow;
02783       LL_indices = LLptr + 2;
02784       int* int_indices = reinterpret_cast<int*>(LL_indices);
02785       EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumIndices, NumIndices, Values.Values(), int_indices));
02786 
02787       // convert to GIDs, start from right.
02788       for(j = NumIndices; j > 0;) {
02789        --j;
02790        LL_indices[j] = colMap.GID64(int_indices[j]);
02791       }
02792 
02793       LLptr[1] = NumIndices; // Load second slot of segment
02794       LLptr += (NumIndices+2); // Point to next segment
02795     }
02796   }
02797   else
02798     throw ReportError("Epetra_CrsGraph::PackAndPrepareRowMatrix: Unable to determine source global index type",-1);
02799 
02800 
02801 //  if( indices ) delete [] indices;
02802 
02803   return(0);
02804 }
02805 
02806 // private =====================================================================
02807 int Epetra_CrsGraph::UnpackAndCombine(const Epetra_SrcDistObject& Source,
02808                                       int NumImportIDs,
02809                                       int* ImportLIDs,
02810                                       int LenImports,
02811                                       char* Imports,
02812                                       int& SizeOfPacket,
02813                                       Epetra_Distributor& Distor,
02814                                       Epetra_CombineMode CombineMode,
02815                                       const Epetra_OffsetIndex * Indexor)
02816 {
02817   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02818     throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Incoming global index type does not match the one for *this",-1);
02819 
02820   (void)Source;
02821   (void)LenImports;
02822   (void)SizeOfPacket;
02823   (void)Distor;
02824   (void)CombineMode;
02825   (void)Indexor;
02826   if(NumImportIDs <= 0)
02827     return(0);
02828 
02829   // Unpack it...
02830 
02831   // Each segment of Sends will be filled by a packed row of information for each row as follows:
02832   // 1st int: GRID of row where GRID is the global row ID for the source graph
02833   // next int:  NumIndices, Number of indices in row.
02834   // next NumIndices: The actual indices for the row.
02835 
02836 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02837   if(Source.Map().GlobalIndicesInt()) {
02838     int NumIndices;
02839     int i;
02840     int* indices;
02841     int ToRow;
02842     int* intptr = (int*) Imports;
02843     for(i = 0; i < NumImportIDs; i++) {
02844       ToRow = (int) GRID64(ImportLIDs[i]);
02845       assert((intptr[0])==ToRow); // Sanity check
02846       NumIndices = intptr[1];
02847       indices = intptr + 2;
02848       // Insert indices
02849       int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
02850       if(ierr < 0)
02851         EPETRA_CHK_ERR(ierr);
02852       intptr += (NumIndices+2); // Point to next segment
02853     }
02854   }
02855   else
02856 #endif
02857 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02858   if(Source.Map().GlobalIndicesLongLong()) {
02859     int NumIndices;
02860     int i;
02861     long long* indices;
02862     long long ToRow;
02863     long long* LLptr = (long long*) Imports;
02864     for(i = 0; i < NumImportIDs; i++) {
02865       ToRow = GRID64(ImportLIDs[i]);
02866       assert((LLptr[0])==ToRow); // Sanity check
02867       NumIndices = (int) LLptr[1];
02868       indices = LLptr + 2;
02869       // Insert indices
02870       int ierr = InsertGlobalIndices(ToRow, NumIndices, indices);
02871       if(ierr < 0)
02872         EPETRA_CHK_ERR(ierr);
02873       LLptr += (NumIndices+2); // Point to next segment
02874     }
02875   }
02876   else
02877 #endif
02878     throw ReportError("Epetra_CrsGraph::UnpackAndCombine: Unable to determine source global index type",-1);
02879 
02880   //destroy buffers since this operation is usually only done once
02881   if( LenExports_ ) {
02882     delete [] Exports_;
02883     Exports_ = 0;
02884     LenExports_ = 0;
02885   }
02886   if( LenImports_ ) {
02887     delete [] Imports_;
02888     Imports_ = 0;
02889     LenImports_ = 0;
02890   }
02891 
02892   return(0);
02893 }
02894 
02895 // protected ===================================================================
02896 bool Epetra_CrsGraph::GlobalConstantsComputed() const {
02897   int mineComputed = 0;
02898   int allComputed;
02899   if(CrsGraphData_->GlobalConstantsComputed_)
02900     mineComputed = 1;
02901   RowMap().Comm().MinAll(&mineComputed, &allComputed, 1); // Find out if any PEs changed constants info
02902   // If allComputed comes back zero then at least one PE need global constants recomputed.
02903   return(allComputed==1);
02904 }
02905 
02906 // private =====================================================================
02907 void Epetra_CrsGraph::ComputeIndexState() {
02908   int myIndicesAreLocal = 0;
02909   int myIndicesAreGlobal = 0;
02910   if(CrsGraphData_->IndicesAreLocal_)
02911     myIndicesAreLocal = 1;
02912   if(CrsGraphData_->IndicesAreGlobal_)
02913     myIndicesAreGlobal = 1;
02914   int allIndicesAreLocal;
02915   int allIndicesAreGlobal;
02916   RowMap().Comm().MaxAll(&myIndicesAreLocal, &allIndicesAreLocal, 1); // Find out if any PEs changed Local Index info
02917   RowMap().Comm().MaxAll(&myIndicesAreGlobal, &allIndicesAreGlobal, 1); // Find out if any PEs changed Global Index info
02918   CrsGraphData_->IndicesAreLocal_ = (allIndicesAreLocal==1); // If indices are local on one PE, should be local on all
02919   CrsGraphData_->IndicesAreGlobal_ = (allIndicesAreGlobal==1);  // If indices are global on one PE should be local on all
02920 }
02921 
02922 //==============================================================================
02923 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
02924 int *Epetra_CrsGraph::All_IndicesPlus1() const {
02925   // This functionality is needed because MKL "sparse matrix" "dense matrix"
02926   // functions do not work with column-based dense storage and zero-based
02927   // sparse storage.  So add "1" to indices and save duplicate data.  This means
02928   // we will use one-based indices.  This does not affect sparse-matrix and vector
02929   // operations.
02930 
02931   int* ptr = 0;
02932   if (!StorageOptimized()) {
02933     throw ReportError("Epetra_CrsGraph: int *All_IndicesPlus1() cannot be called when StorageOptimized()==false", -1);
02934   }
02935   else {
02936     Epetra_IntSerialDenseVector& vec = CrsGraphData_->data->All_IndicesPlus1_;
02937 
02938     if(!vec.Length()) {
02939       int* indices = All_Indices();
02940     vec.Size(CrsGraphData_->data->All_Indices_.Length());
02941     ptr = vec.Values();
02942       for(int i = 0; i < CrsGraphData_->NumMyNonzeros_; ++i)
02943       ptr[i] = indices[i] + 1;
02944   }
02945   else {
02946     ptr = vec.Values();
02947   }
02948   }
02949   return ptr;
02950 }
02951 #endif // defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
02952 
02953 //==============================================================================
02954 void Epetra_CrsGraph::Print (std::ostream& os) const {
02955   int MyPID = RowMap().Comm().MyPID();
02956   int NumProc = RowMap().Comm().NumProc();
02957 
02958   for(int iproc = 0; iproc < NumProc; iproc++) {
02959     if(MyPID == iproc) {
02960       if(MyPID == 0) {
02961   os << "\nNumber of Global Block Rows  = " << NumGlobalBlockRows64()      << std::endl;
02962   os <<   "Number of Global Block Cols  = " << NumGlobalBlockCols64()      << std::endl;
02963   os <<   "Number of Global Block Diags = " << NumGlobalBlockDiagonals64() << std::endl;
02964   os <<   "Number of Global Entries     = " << NumGlobalEntries64()        << std::endl;
02965   os << "\nNumber of Global Rows        = " << NumGlobalRows64()           << std::endl;
02966   os <<   "Number of Global Cols        = " << NumGlobalCols64()           << std::endl;
02967   os <<   "Number of Global Diagonals   = " << NumGlobalDiagonals64()      << std::endl;
02968   os <<   "Number of Global Nonzeros    = " << NumGlobalNonzeros64()       << std::endl;
02969   os << "\nGlobal Maximum Block Row Dim = " << GlobalMaxRowDim()         << std::endl;
02970   os <<   "Global Maximum Block Col Dim = " << GlobalMaxColDim()         << std::endl;
02971   os <<   "Global Maximum Num Indices   = " << GlobalMaxNumIndices()     << std::endl;
02972   if(LowerTriangular()) os << " ** Matrix is Lower Triangular **"        << std::endl;
02973   if(UpperTriangular()) os << " ** Matrix is Upper Triangular **"        << std::endl;
02974   if(NoDiagonal())      os << " ** Matrix has no diagonal     **"        << std::endl << std::endl;
02975       }
02976       os << "\nNumber of My Block Rows  = " << NumMyBlockRows()      << std::endl;
02977       os <<   "Number of My Block Cols  = " << NumMyBlockCols()      << std::endl;
02978       os <<   "Number of My Block Diags = " << NumMyBlockDiagonals() << std::endl;
02979       os <<   "Number of My Entries     = " << NumMyEntries()        << std::endl;
02980       os << "\nNumber of My Rows        = " << NumMyRows()           << std::endl;
02981       os <<   "Number of My Cols        = " << NumMyCols()           << std::endl;
02982       os <<   "Number of My Diagonals   = " << NumMyDiagonals()      << std::endl;
02983       os <<   "Number of My Nonzeros    = " << NumMyNonzeros()       << std::endl;
02984       os << "\nMy Maximum Block Row Dim = " << MaxRowDim()           << std::endl;
02985       os <<   "My Maximum Block Col Dim = " << MaxColDim()           << std::endl;
02986       os <<   "My Maximum Num Indices   = " << MaxNumIndices()       << std::endl << std::endl;
02987 
02988       int NumMyBlockRows1 = NumMyBlockRows();
02989       int MaxNumIndices1 = MaxNumIndices();
02990       Epetra_IntSerialDenseVector Indices1_int(MaxNumIndices1);
02991 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02992       Epetra_LongLongSerialDenseVector Indices1_LL(MaxNumIndices1);
02993 #endif
02994 
02995       if(RowMap().GlobalIndicesInt()) {
02996      Indices1_int.Resize(MaxNumIndices1);
02997       }
02998       else if(RowMap().GlobalIndicesLongLong()) {
02999 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
03000      Indices1_LL.Resize(MaxNumIndices1);
03001 #else
03002          throw ReportError("Epetra_CrsGraph::Print: GlobalIndicesLongLong but no long long API",-1);
03003 #endif
03004       }
03005       else
03006          throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
03007 
03008     int NumIndices1;
03009       int i;
03010       int j;
03011 
03012       os.width(14);
03013       os <<  "       Row Index "; os << " ";
03014       for(j = 0; j < MaxNumIndices(); j++) {
03015   os.width(12);
03016   os <<  "Col Index"; os << "      ";
03017       }
03018       os << std::endl;
03019       for(i = 0; i < NumMyBlockRows1; i++) {
03020        if(RowMap().GlobalIndicesInt()) {
03021          int Row = (int) GRID64(i); // Get global row number
03022          ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_int.Values());
03023          os.width(14);
03024          os <<  Row ; os << "    ";
03025          for(j = 0; j < NumIndices1 ; j++) {
03026            os.width(12);
03027            os <<  Indices1_int[j]; os << "    ";
03028          }
03029          os << std::endl;
03030        }
03031        else if(RowMap().GlobalIndicesLongLong()) {
03032 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
03033          long long Row = GRID64(i); // Get global row number
03034          ExtractGlobalRowCopy(Row, MaxNumIndices1, NumIndices1, Indices1_LL.Values());
03035          os.width(14);
03036          os <<  Row ; os << "    ";
03037          for(j = 0; j < NumIndices1 ; j++) {
03038            os.width(12);
03039            os <<  Indices1_LL[j]; os << "    ";
03040          }
03041          os << std::endl;
03042 #else
03043          throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
03044 #endif
03045        }
03046       }
03047       os << std::flush;
03048     }
03049     // Do a few global ops to give I/O a chance to complete
03050     RowMap().Comm().Barrier();
03051     RowMap().Comm().Barrier();
03052     RowMap().Comm().Barrier();
03053   }
03054 }
03055 
03056 //==============================================================================
03057 Epetra_CrsGraph& Epetra_CrsGraph::operator = (const Epetra_CrsGraph& Source) {
03058   if ((this == &Source) || (CrsGraphData_ == Source.CrsGraphData_))
03059     return(*this); // this and Source are same Graph object, or both point to same CrsGraphData object
03060 
03061   CleanupData();
03062   CrsGraphData_ = Source.CrsGraphData_;
03063   CrsGraphData_->IncrementReferenceCount();
03064 
03065   return(*this);
03066 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines