Epetra Package Browser (Single Doxygen Collection) Development
Epetra_CrsGraphData.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_CrsGraphData.h"
00045 #include "Epetra_Import.h"
00046 #include "Epetra_Export.h"
00047 //#include "Epetra_ConfigDefs.h" //DATA_DEBUG
00048 
00049 //=============================================================================
00050 Epetra_CrsGraphData::Epetra_CrsGraphData(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, bool StaticProfile)
00051   // maps
00052   : RowMap_(RowMap),
00053     ColMap_(RowMap),
00054     DomainMap_(RowMap),
00055     RangeMap_(RowMap),
00056     // importer & exporter
00057     Importer_(0),
00058     Exporter_(0),
00059     // booleans
00060     HaveColMap_(false),
00061     Filled_(false),
00062     Allocated_(false),
00063     // for non-static profile, we insert always into sorted lists, so the
00064     // graph will always be sorted. The same holds for the redundancies.
00065     Sorted_(!StaticProfile),
00066     StorageOptimized_(false),
00067     NoRedundancies_(!StaticProfile),
00068     IndicesAreGlobal_(false),
00069     IndicesAreLocal_(false),
00070     IndicesAreContiguous_(false),
00071     LowerTriangular_(true),
00072     UpperTriangular_(true),
00073     NoDiagonal_(true),
00074     GlobalConstantsComputed_(false),
00075     StaticProfile_(StaticProfile),
00076     SortGhostsAssociatedWithEachProcessor_(false),
00077 
00078     // ints
00079     IndexBase_(RowMap.IndexBase()),
00080     NumGlobalEntries_(0),
00081     NumGlobalBlockRows_(RowMap.NumGlobalElements()),
00082     NumGlobalBlockCols_(NumGlobalBlockRows_),
00083     NumGlobalBlockDiagonals_(0),
00084     NumMyEntries_(0),
00085     NumMyBlockRows_(RowMap.NumMyElements()),
00086     NumMyBlockCols_(NumMyBlockRows_),
00087     NumMyBlockDiagonals_(0),
00088     MaxRowDim_(RowMap.MaxElementSize()),
00089     MaxColDim_(MaxRowDim_),
00090     GlobalMaxRowDim_(RowMap.MaxElementSize()),
00091     GlobalMaxColDim_(GlobalMaxRowDim_),
00092     MaxNumNonzeros_(0),
00093     GlobalMaxNumNonzeros_(0),
00094     NumGlobalNonzeros_(0),
00095     NumGlobalRows_(RowMap.NumGlobalPoints()),
00096     NumGlobalCols_(NumGlobalRows_),
00097     NumGlobalDiagonals_(0),
00098     NumMyNonzeros_(0),
00099     NumMyRows_(RowMap.NumMyPoints()),
00100     NumMyCols_(NumMyRows_),
00101     NumMyDiagonals_(0), 
00102     MaxNumIndices_(0),
00103     GlobalMaxNumIndices_(0),
00104     Indices_(new int *[NumMyBlockRows_]),
00105     SortedEntries_(StaticProfile ? 0 : NumMyBlockRows_),
00106     TempColIndices_(NULL),
00107     NumTempColIndices_(0),
00108     NumAllocatedIndicesPerRow_(0),
00109     NumIndicesPerRow_(0),
00110     IndexOffset_(0),
00111     All_Indices_(0),
00112     CV_(CV)
00113 {
00114   //cout << "--CRSGD created(rowmap ctr), addr: " << this << endl; //DATA_DEBUG
00115 }
00116 
00117 //=============================================================================
00118 Epetra_CrsGraphData::Epetra_CrsGraphData(Epetra_DataAccess CV, 
00119            const Epetra_BlockMap& RowMap, 
00120            const Epetra_BlockMap& ColMap, bool StaticProfile)
00121   // maps
00122   : RowMap_(RowMap),
00123     ColMap_(ColMap),
00124     DomainMap_(ColMap),
00125     RangeMap_(RowMap),
00126     // importer & exporter
00127     Importer_(0),
00128     Exporter_(0),
00129     // booleans
00130     HaveColMap_(true),
00131     Filled_(false),
00132     Allocated_(false),
00133     Sorted_(!StaticProfile),
00134     StorageOptimized_(false),
00135     NoRedundancies_(!StaticProfile),
00136     IndicesAreGlobal_(false),
00137     IndicesAreLocal_(false),
00138     IndicesAreContiguous_(false),
00139     LowerTriangular_(true),
00140     UpperTriangular_(true),
00141     NoDiagonal_(true),
00142     GlobalConstantsComputed_(false),
00143     StaticProfile_(StaticProfile),
00144     SortGhostsAssociatedWithEachProcessor_(false),
00145     // ints
00146     IndexBase_(RowMap.IndexBase()),
00147     NumGlobalEntries_(0),
00148     NumGlobalBlockRows_(RowMap.NumGlobalElements()),
00149     NumGlobalBlockCols_(ColMap.NumGlobalElements()),
00150     NumGlobalBlockDiagonals_(0),
00151     NumMyEntries_(0),
00152     NumMyBlockRows_(RowMap.NumMyElements()),
00153     NumMyBlockCols_(ColMap.NumMyElements()),
00154     NumMyBlockDiagonals_(0),
00155     MaxRowDim_(RowMap.MaxElementSize()),
00156     MaxColDim_(ColMap.MaxElementSize()),
00157     GlobalMaxRowDim_(RowMap.MaxElementSize()),
00158     GlobalMaxColDim_(ColMap.MaxElementSize()),
00159     MaxNumNonzeros_(0),
00160     GlobalMaxNumNonzeros_(0),
00161     NumGlobalNonzeros_(0),
00162     NumGlobalRows_(RowMap.NumGlobalPoints()),
00163     NumGlobalCols_(ColMap.NumGlobalPoints()),
00164     NumGlobalDiagonals_(0),
00165     NumMyNonzeros_(0),
00166     NumMyRows_(RowMap.NumMyPoints()),
00167     NumMyCols_(ColMap.NumMyPoints()),
00168     NumMyDiagonals_(0), 
00169     MaxNumIndices_(0),
00170     GlobalMaxNumIndices_(0),
00171     Indices_(new int *[NumMyBlockRows_]),
00172     SortedEntries_(StaticProfile ? 0 : NumMyBlockRows_),
00173     TempColIndices_(NULL),
00174     NumTempColIndices_(0),
00175     NumAllocatedIndicesPerRow_(0),
00176     NumIndicesPerRow_(0),
00177     IndexOffset_(0),
00178     All_Indices_(0),
00179     CV_(CV)
00180 {
00181   //cout << "--CRSGD created(rowmap&colmap ctr), addr: " << this << endl; //DATA_DEBUG
00182 }
00183 
00184 //=============================================================================
00185 Epetra_CrsGraphData::~Epetra_CrsGraphData() {
00186 
00187   if(Indices_ != 0 && !StorageOptimized_) {
00188     for (int i=0; i<NumMyBlockRows_; i++) {
00189       Indices_[i] = 0;
00190     } 
00191     delete[] Indices_;
00192     Indices_ = 0;
00193   }
00194 
00195   if (TempColIndices_ != 0) {
00196     delete [] TempColIndices_;
00197     TempColIndices_ = 0;
00198   }
00199 
00200   if(Importer_ != 0) {
00201     delete Importer_;
00202     Importer_ = 0;
00203   }
00204   if(Exporter_ != 0) {
00205     delete Exporter_;
00206     Importer_ = 0;
00207   }
00208 
00209   NumMyBlockRows_ = 0;  // are these needed?
00210   Filled_ = false;      // they're about to go out of scope, after all
00211   Allocated_ = false;
00212 
00213   //cout << "--CRSGD destroyed, addr: " << this << endl; //DATA_DEBUG
00214 }
00215 
00216 //==========================================================================
00217 int Epetra_CrsGraphData::MakeImportExport() {
00218   // Create Import object for use by matrix classes.    This is only needed if ColMap and DomainMap are different
00219   if (!ColMap_.SameAs(DomainMap_)) {
00220     if (Importer_ != 0) {
00221       delete Importer_;
00222       Importer_ = 0;
00223     }
00224     Importer_ = new Epetra_Import(ColMap_, DomainMap_);
00225   }
00226   
00227   // Now see if we need to define an export map.  This is only needed if RowMap and RangeMap are different
00228   if (!RowMap_.SameAs(RangeMap_)) {
00229     if (Exporter_ != 0) {
00230       delete Exporter_;
00231       Exporter_ = 0;
00232     }
00233     Exporter_ = new Epetra_Export(RowMap_, RangeMap_); // Create Export object. 
00234   }
00235    
00236   return(0);
00237 }
00238 
00239 //==========================================================================
00240 int Epetra_CrsGraphData::ReAllocateAndCast(char*& UserPtr, int& Length, const int IntPacketSizeTimesNumTrans) {
00241   if(IntPacketSizeTimesNumTrans > Length) {
00242     if(Length > 0) 
00243       delete[] UserPtr;
00244     Length = IntPacketSizeTimesNumTrans;
00245     int* newPtr = new int[Length];
00246     UserPtr = reinterpret_cast<char*> (newPtr);
00247   }
00248   return(0);
00249 }
00250 
00251 //==========================================================================
00252 void Epetra_CrsGraphData::Print(ostream& os, int level) const {
00253   bool four_bit = (level >= 4);      // 4-bit = BlockMaps
00254   bool two_bit = ((level % 4) >= 2); // 2-bit = Indices
00255   bool one_bit = ((level % 2) == 1); // 1-bit = Everything else
00256 
00257   os << "\n***** CrsGraphData (output level " << level << ") *****" << endl;
00258 
00259   if(four_bit) {
00260     os << "RowMap_:\n" << RowMap_ << endl;
00261     os << "ColMap_:\n" << ColMap_ << endl;
00262     os << "DomainMap_:\n" << DomainMap_ << endl;
00263     os << "RangeMap_:\n" << RangeMap_ << endl;
00264   }
00265   
00266   if(one_bit) {
00267     os.width(26); os << "HaveColMap_: "              << HaveColMap_;
00268     os.width(25); os << "Filled_: "                  << Filled_;
00269     os.width(25); os << "Allocated_: "               << Allocated_;
00270     os.width(25); os << "Sorted_: "                  << Sorted_ << endl;
00271     os.width(26); os << "StorageOptimized_: "        << StorageOptimized_;
00272     os.width(25); os << "SortGhostsAssociatedWithEachProcessor_: " << SortGhostsAssociatedWithEachProcessor_;
00273     os.width(25); os << "NoRedundancies_: "          << NoRedundancies_;
00274     os.width(25); os << "IndicesAreGlobal_: "        << IndicesAreGlobal_;
00275     os.width(25); os << "IndicesAreLocal_: "         << IndicesAreLocal_ << endl;
00276     os.width(26); os << "IndicesAreContiguous_: "    << IndicesAreContiguous_;
00277     os.width(25); os << "LowerTriangular_: "         << LowerTriangular_;
00278     os.width(25); os << "UpperTriangular_: "         << UpperTriangular_;
00279     os.width(25); os << "NoDiagonal_: "              << NoDiagonal_ << endl;
00280     os.width(25); os << "GlobalConstantsComputed_: " << GlobalConstantsComputed_ << endl;
00281     os.width(25); os << "StaticProfile_: " << StaticProfile_ << endl << endl;
00282     
00283     os.width(10); os << "NGBR_: " << NumGlobalBlockRows_;
00284     os.width(10); os << "NGBC_: " << NumGlobalBlockCols_;
00285     os.width(10); os << "NGBD_: " << NumGlobalBlockDiagonals_;
00286     os.width(10); os << "NGE_: "  << NumGlobalEntries_;
00287     os.width(10); os << "NGR_: "  << NumGlobalRows_;
00288     os.width(10); os << "NGC_: "  << NumGlobalCols_;
00289     os.width(10); os << "NGD_: "  << NumGlobalDiagonals_;
00290     os.width(10); os << "NGN_: "  << NumGlobalNonzeros_;
00291     os.width(10); os << "IB_: "   << IndexBase_ << endl;
00292     os.width(10); os << "GMRD_: " << GlobalMaxRowDim_;
00293     os.width(11); os << "GMCD_: " << GlobalMaxColDim_;
00294     os.width(11); os << "GMNI_: " << GlobalMaxNumIndices_;
00295     os.width(11); os << "NMBR_: " << NumMyBlockRows_;
00296     os.width(10); os << "NMBC_: " << NumMyBlockCols_;
00297     os.width(10); os << "NMBD_: " << NumMyBlockDiagonals_;
00298     os.width(10); os << "NME_: "  << NumMyEntries_;
00299     os.width(10); os << "NMR_: "  << NumMyRows_;
00300     os.width(10); os << "CV_: " << CV_ << endl;
00301     os.width(10); os << "NMC_: "  << NumMyCols_;
00302     os.width(10); os << "NMD_: "  << NumMyDiagonals_;
00303     os.width(10); os << "NMN_: "  << NumMyNonzeros_;
00304     os.width(10); os << "MRD_: "  << MaxRowDim_;
00305     os.width(11); os << "MCD_: "  << MaxColDim_;
00306     os.width(11); os << "MNI_: "  << MaxNumIndices_;
00307     os.width(11); os << "MNN_: "  << MaxNumNonzeros_;
00308     os.width(11); os << "GMNN_: " << GlobalMaxNumNonzeros_;
00309     os.width(11); os << "RC: " << ReferenceCount() << endl << endl;
00310     
00311     os << "NIPR_: " << NumIndicesPerRow_ << endl;
00312     os << "NAIPR_: " << NumAllocatedIndicesPerRow_ << endl;
00313     os << "IndexOffset_: " << IndexOffset_ << endl;
00314     os << "All_Indices_: " << All_Indices_ << endl;
00315   }
00316     
00317   if(two_bit) {
00318     os << "Indices_: " << Indices_ << endl;
00319     if(Indices_ != 0) {
00320       for(int i = 0; i < NumMyBlockRows_; i++) {
00321   os << "Indices_[" << i << "]: (" << Indices_[i] << ") ";
00322   if(Indices_[i] != 0) {
00323     for(int j = 0; j < NumAllocatedIndicesPerRow_[i]; j++)
00324       os << Indices_[i][j] << " ";
00325   }
00326   os << endl;
00327       }
00328     }
00329   }
00330   
00331   os << "***** End CrsGraphData *****" << endl;
00332 }
00333 
00334 
00335 //==========================================================================
00336 void
00337 Epetra_CrsGraphData::EntriesInOneRow::AddEntry (const int Col)
00338 {
00339   // first check the last element (or if line is still empty)
00340   if ( (entries_.size()==0) || ( entries_.back() < Col) )
00341     {
00342       entries_.push_back(Col);
00343       return;
00344     }
00345  
00346   // do a binary search to find the place where to insert:
00347   std::vector<int>::iterator it = std::lower_bound(entries_.begin(),
00348                 entries_.end(),
00349                 Col);
00350  
00351   // If this entry is a duplicate, exit immediately
00352   if (*it == Col)
00353     return;
00354  
00355   // Insert at the right place in the vector. Vector grows automatically to
00356   // fit elements. Always doubles its size.
00357   entries_.insert(it, Col);
00358 }
00359 
00360 
00361 //==========================================================================
00362 void
00363 Epetra_CrsGraphData::EntriesInOneRow::AddEntries (const int  numCols,
00364              const int *Indices)
00365 {
00366   if (numCols == 0)
00367     return;
00368 
00369   // Check whether the indices are sorted. Can do more efficient then.
00370   bool indicesAreSorted = true;
00371   for (int i=1; i<numCols; ++i)
00372     if (Indices[i] <= Indices[i-1]) {
00373       indicesAreSorted = false;
00374       break;
00375     }
00376 
00377   if (indicesAreSorted && numCols > 3) {
00378     const int * curInput = &Indices[0];
00379     int col = *curInput;
00380     const int * endInput = &Indices[numCols];
00381 
00382     // easy case: list of entries is empty or all entries are smaller than
00383     // the ones to be inserted
00384     if (entries_.size() == 0 || entries_.back() < col)
00385     {
00386       entries_.insert(entries_.end(), &Indices[0], &Indices[numCols]);
00387       return;
00388     }
00389 
00390     // find a possible insertion point for the first entry. check whether
00391     // the first entry is a duplicate before actually doing something.
00392     std::vector<int>::iterator it = 
00393       std::lower_bound(entries_.begin(), entries_.end(), col);
00394     while (*it == col) {
00395       ++curInput;
00396       if (curInput == endInput)
00397         break;
00398       col = *curInput;
00399 
00400       // check the very next entry in the current array
00401       ++it;
00402       if (it == entries_.end())
00403         break;
00404       if (*it > col)
00405         break;
00406       if (*it == col)
00407         continue;
00408 
00409       // ok, it wasn't the very next one, do a binary search to find the
00410       // insert point
00411       it = std::lower_bound(it, entries_.end(), col);
00412       if (it == entries_.end())
00413         break;
00414     }
00415 
00416     // all input entries were duplicates.
00417     if (curInput == endInput)
00418       return;
00419 
00420     // Resize vector by just inserting the list at the correct point. Note
00421     // that the list will not yet be sorted, but rather have the insert
00422     // indices in the middle and the old indices from the list on the
00423     // end. Next we will have to merge the two lists.
00424     const int pos1 = it - entries_.begin();
00425     entries_.insert (it, curInput, endInput);
00426     it = entries_.begin() + pos1;
00427 
00428     // Now merge the two lists...
00429     std::vector<int>::iterator it2 = it + (endInput - curInput);
00430 
00431     // As long as there are indices both in the end of the entries list and
00432     // in the input list, always continue with the smaller index.
00433     while (curInput != endInput && it2 != entries_.end())
00434     {
00435       if (*curInput < *it2)
00436         *it++ = *curInput++;
00437       else if (*curInput == *it2)
00438       {
00439         *it++ = *it2++;
00440         ++curInput;
00441       }
00442       else
00443         *it++ = *it2++;
00444     }
00445     // in case there are indices left in the input list or in the end of
00446     // entries, we have to use the entries that are left. Only one of the
00447     // two while loops will actually be entered.
00448     while (curInput != endInput)
00449       *it++ = *curInput++;
00450 
00451     while (it2 != entries_.end())
00452       *it++ = *it2++;
00453 
00454     // resize and return
00455     const int new_size = it - entries_.begin();
00456     entries_.resize (new_size);
00457     return;
00458   }
00459 
00460   // for unsorted or just a few indices, go to the one-by-one entry
00461   // function.
00462   for (int i=0; i<numCols; ++i)
00463     AddEntry(Indices[i]);
00464 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines