|
Tpetra Matrix/Vector Services
Version of the Day
|
|
Tpetra Matrix/Vector Services
Version of the Day
|
00001 // @HEADER 00002 // *********************************************************************** 00003 // 00004 // Tpetra: Templated Linear Algebra Services Package 00005 // Copyright (2008) Sandia Corporation 00006 // 00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, 00008 // the U.S. Government retains certain rights in this software. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00038 // 00039 // ************************************************************************ 00040 // @HEADER 00041 00042 #ifndef TPETRA_CRSGRAPH_DEF_HPP 00043 #define TPETRA_CRSGRAPH_DEF_HPP 00044 00045 #include <Tpetra_Distributor.hpp> 00046 #include <Teuchos_Assert.hpp> 00047 #include <Teuchos_NullIteratorTraits.hpp> 00048 #include <Teuchos_as.hpp> 00049 #include <algorithm> 00050 #include <string> 00051 #include <utility> 00052 #include <Teuchos_SerialDenseMatrix.hpp> 00053 00054 #ifdef DOXYGEN_USE_ONLY 00055 #include "Tpetra_CrsGraph_decl.hpp" 00056 #endif 00057 00058 namespace Tpetra { 00059 00062 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00063 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00064 CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00065 size_t maxNumEntriesPerRow, 00066 ProfileType pftype, 00067 const RCP<ParameterList>& params) 00068 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00069 , rowMap_(rowMap) 00070 , nodeNumEntries_(0) 00071 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00072 , pftype_(pftype) 00073 , numAllocForAllRows_(maxNumEntriesPerRow) 00074 , indicesAreAllocated_(false) 00075 , indicesAreLocal_(false) 00076 , indicesAreGlobal_(false) 00077 , fillComplete_(false) 00078 , indicesAreSorted_ (true) 00079 , noRedundancies_ (true) 00080 , haveLocalConstants_ (false) 00081 , haveGlobalConstants_ (false) 00082 , sortGhostsAssociatedWithEachProcessor_(true) 00083 { 00084 typedef Teuchos::OrdinalTraits<size_t> OTST; 00085 staticAssertions(); 00086 TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(), 00087 std::invalid_argument, "The allocation hint must be a valid size_t value, " 00088 "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::" 00089 "invalid()."); 00090 resumeFill(params); 00091 checkInternalState(); 00092 } 00093 00096 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00097 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00098 CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00099 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 00100 size_t maxNumEntriesPerRow, 00101 ProfileType pftype, 00102 const RCP<ParameterList>& params) 00103 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00104 , rowMap_(rowMap) 00105 , colMap_(colMap) 00106 , nodeNumEntries_(0) 00107 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00108 , pftype_(pftype) 00109 , numAllocForAllRows_(maxNumEntriesPerRow) 00110 , indicesAreAllocated_(false) 00111 , indicesAreLocal_(false) 00112 , indicesAreGlobal_(false) 00113 , fillComplete_(false) 00114 , indicesAreSorted_(true) 00115 , noRedundancies_(true) 00116 , haveLocalConstants_ (false) 00117 , haveGlobalConstants_ (false) 00118 , sortGhostsAssociatedWithEachProcessor_(true) 00119 { 00120 typedef Teuchos::OrdinalTraits<size_t> OTST; 00121 staticAssertions(); 00122 TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(), 00123 std::invalid_argument, "The allocation hint must be a valid size_t value, " 00124 "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::" 00125 "invalid()."); 00126 resumeFill(params); 00127 checkInternalState(); 00128 } 00129 00132 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00133 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00134 CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00135 const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, 00136 ProfileType pftype, 00137 const RCP<ParameterList>& params) 00138 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00139 , rowMap_(rowMap) 00140 , nodeNumEntries_(0) 00141 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00142 , pftype_(pftype) 00143 , numAllocPerRow_(NumEntriesPerRowToAlloc) 00144 , numAllocForAllRows_(0) 00145 , indicesAreAllocated_(false) 00146 , indicesAreLocal_(false) 00147 , indicesAreGlobal_(false) 00148 , fillComplete_(false) 00149 , indicesAreSorted_(true) 00150 , noRedundancies_(true) 00151 , haveLocalConstants_ (false) 00152 , haveGlobalConstants_ (false) 00153 , sortGhostsAssociatedWithEachProcessor_(true) 00154 { 00155 typedef Teuchos::OrdinalTraits<size_t> OTST; 00156 const char tfecfFuncName[] = "CrsGraph(rowMap,NumEntriesPerRowToAlloc)"; 00157 staticAssertions(); 00158 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument, 00159 ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node."); 00160 for (size_t r=0; r < getNodeNumRows(); ++r) { 00161 const size_t curRowCount = NumEntriesPerRowToAlloc[r]; 00162 TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(), 00163 std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies " 00164 "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::" 00165 "invalid())."); 00166 } 00167 resumeFill(params); 00168 checkInternalState(); 00169 } 00170 00171 00174 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00175 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00176 CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00177 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 00178 const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, 00179 ProfileType pftype, 00180 const RCP<ParameterList>& params) 00181 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00182 , rowMap_(rowMap) 00183 , colMap_(colMap) 00184 , nodeNumEntries_(0) 00185 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00186 , pftype_(pftype) 00187 , numAllocPerRow_(NumEntriesPerRowToAlloc) 00188 , numAllocForAllRows_(0) 00189 , indicesAreAllocated_(false) 00190 , indicesAreLocal_(false) 00191 , indicesAreGlobal_(false) 00192 , fillComplete_(false) 00193 , indicesAreSorted_(true) 00194 , noRedundancies_(true) 00195 , haveLocalConstants_ (false) 00196 , haveGlobalConstants_ (false) 00197 , sortGhostsAssociatedWithEachProcessor_(true) 00198 { 00199 typedef Teuchos::OrdinalTraits<size_t> OTST; 00200 const char tfecfFuncName[] = "CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc)"; 00201 staticAssertions(); 00202 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument, 00203 ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node."); 00204 for (size_t r=0; r < getNodeNumRows(); ++r) { 00205 const size_t curRowCount = NumEntriesPerRowToAlloc[r]; 00206 TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(), 00207 std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies " 00208 "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::" 00209 "invalid())."); 00210 } 00211 resumeFill(params); 00212 checkInternalState(); 00213 } 00214 00215 00218 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00219 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00220 CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00221 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 00222 const ArrayRCP<size_t> & rowPointers, 00223 const ArrayRCP<LocalOrdinal> & columnIndices, 00224 const RCP<ParameterList>& params) 00225 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00226 , rowMap_(rowMap) 00227 , colMap_(colMap) 00228 , nodeNumEntries_(0) 00229 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00230 , pftype_(StaticProfile) 00231 , numAllocForAllRows_(0) 00232 , indicesAreAllocated_(true) 00233 , indicesAreLocal_(true) 00234 , indicesAreGlobal_(false) 00235 , fillComplete_(false) 00236 , indicesAreSorted_(true) 00237 , noRedundancies_(true) 00238 , haveLocalConstants_ (false) 00239 , haveGlobalConstants_ (false) 00240 , sortGhostsAssociatedWithEachProcessor_(true) 00241 { 00242 staticAssertions(); 00243 globalNumEntries_ = globalNumDiags_ = globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid(); 00244 setAllIndices(rowPointers,columnIndices); 00245 checkInternalState(); 00246 } 00247 00250 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00251 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::~CrsGraph() 00252 {} 00253 00254 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00255 RCP<const ParameterList> 00256 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00257 getValidParameters () const 00258 { 00259 RCP<ParameterList> params = parameterList ("Tpetra::CrsGraph"); 00260 00261 // Make a sublist for the Import. 00262 RCP<ParameterList> importSublist = parameterList ("Import"); 00263 00264 // FIXME (mfh 02 Apr 2012) We should really have the Import and 00265 // Export objects fill in these lists. However, we don't want to 00266 // create an Import or Export unless we need them. For now, we 00267 // know that the Import and Export just pass the list directly to 00268 // their Distributor, so we can create a Distributor here 00269 // (Distributor's constructor is a lightweight operation) and have 00270 // it fill in the list. 00271 00272 // Fill in Distributor default parameters by creating a 00273 // Distributor and asking it to do the work. 00274 Distributor distributor (rowMap_->getComm(), importSublist); 00275 params->set ("Import", *importSublist, "How the Import performs communication."); 00276 00277 // Make a sublist for the Export. For now, it's a clone of the 00278 // Import sublist. It's not a shallow copy, though, since we 00279 // might like the Import to do communication differently than the 00280 // Export. 00281 params->set ("Export", *importSublist, "How the Export performs communication."); 00282 00283 return params; 00284 } 00285 00286 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00287 void 00288 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00289 setParameterList (const RCP<ParameterList>& params) 00290 { 00291 RCP<const ParameterList> validParams = getValidParameters (); 00292 params->validateParametersAndSetDefaults (*validParams); 00293 this->setMyParamList (params); 00294 } 00295 00296 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00297 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumRows() const 00298 { 00299 return rowMap_->getGlobalNumElements(); 00300 } 00301 00302 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00303 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumCols() const 00304 { 00305 const char tfecfFuncName[] = "getGlobalNumCols()"; 00306 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00307 isFillComplete() == false, std::runtime_error, 00308 ": requires domain map, which requires that fillComplete() has been " 00309 "called."); 00310 return getDomainMap()->getGlobalNumElements(); 00311 } 00312 00313 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00314 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumRows() const 00315 { 00316 return rowMap_->getNodeNumElements(); 00317 } 00318 00319 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00320 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumCols() const 00321 { 00322 const char tfecfFuncName[] = "getNodeNumCols()"; 00323 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00324 hasColMap() == false, std::runtime_error, 00325 ": requires column map. This requires either that a custom column Map " 00326 "was given to the constructor, or that fillComplete() has been called."); 00327 return colMap_->getNodeNumElements(); 00328 } 00329 00330 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00331 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumDiags() const 00332 { 00333 return nodeNumDiags_; 00334 } 00335 00336 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00337 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumDiags() const 00338 { 00339 return globalNumDiags_; 00340 } 00341 00342 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00343 RCP<Node> 00344 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNode() const 00345 { 00346 return rowMap_->getNode(); 00347 } 00348 00349 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00350 RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > 00351 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRowMap() const 00352 { 00353 return rowMap_; 00354 } 00355 00356 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00357 RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > 00358 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getColMap() const 00359 { 00360 return colMap_; 00361 } 00362 00363 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00364 RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > 00365 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const 00366 { 00367 return domainMap_; 00368 } 00369 00370 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00371 RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > 00372 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const 00373 { 00374 return rangeMap_; 00375 } 00376 00377 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00378 RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > 00379 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getImporter() const 00380 { 00381 return importer_; 00382 } 00383 00384 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00385 RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > 00386 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getExporter() const 00387 { 00388 return exporter_; 00389 } 00390 00391 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00392 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::hasColMap() const 00393 { 00394 return (colMap_ != null); 00395 } 00396 00397 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00398 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isStorageOptimized() const 00399 { 00400 bool isOpt = indicesAreAllocated_ && (numRowEntries_ == null) && (getNodeNumRows() > 0); 00401 #ifdef HAVE_TPETRA_DEBUG 00402 const char tfecfFuncName[] = "isStorageOptimized()"; 00403 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (isOpt == true) && (getProfileType() == DynamicProfile), std::logic_error, 00404 ": Violated stated post-conditions. Please contact Tpetra team."); 00405 #endif 00406 return isOpt; 00407 } 00408 00409 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00410 ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getProfileType() const 00411 { 00412 return pftype_; 00413 } 00414 00415 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00416 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumEntries() const 00417 { 00418 return globalNumEntries_; 00419 } 00420 00421 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00422 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumEntries() const 00423 { 00424 return nodeNumEntries_; 00425 } 00426 00427 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00428 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalMaxNumRowEntries() const 00429 { 00430 return globalMaxNumRowEntries_; 00431 } 00432 00433 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00434 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeMaxNumRowEntries() const 00435 { 00436 return nodeMaxNumRowEntries_; 00437 } 00438 00439 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00440 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isFillComplete() const 00441 { 00442 return fillComplete_; 00443 } 00444 00445 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00446 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isFillActive() const 00447 { 00448 return !fillComplete_; 00449 } 00450 00451 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00452 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isUpperTriangular() const 00453 { 00454 return upperTriangular_; 00455 } 00456 00457 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00458 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isLowerTriangular() const 00459 { 00460 return lowerTriangular_; 00461 } 00462 00463 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00464 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isLocallyIndexed() const 00465 { 00466 return indicesAreLocal_; 00467 } 00468 00469 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00470 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isGloballyIndexed() const 00471 { 00472 return indicesAreGlobal_; 00473 } 00474 00475 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00476 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeAllocationSize() const 00477 { 00478 return nodeNumAllocated_; 00479 } 00480 00481 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00482 RCP<const Comm<int> > 00483 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getComm() const 00484 { 00485 return rowMap_->getComm(); 00486 } 00487 00488 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00489 GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getIndexBase() const 00490 { 00491 return rowMap_->getIndexBase(); 00492 } 00493 00494 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00495 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::indicesAreAllocated() const 00496 { 00497 return indicesAreAllocated_; 00498 } 00499 00500 00503 // // 00504 // Internal utility methods // 00505 // // 00508 00509 00512 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00513 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isSorted() const 00514 { 00515 return indicesAreSorted_; 00516 } 00517 00518 00521 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00522 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isMerged() const 00523 { 00524 return noRedundancies_; 00525 } 00526 00529 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00530 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::setLocallyModified () 00531 { 00532 // FIXME (mfh 07 May 2013) How do we know that the change 00533 // introduced a redundancy, or even that it invalidated the sorted 00534 // order of indices? CrsGraph has always made this conservative 00535 // guess. It could be a bit costly to check at insertion time, 00536 // though. 00537 indicesAreSorted_ = false; 00538 noRedundancies_ = false; 00539 00540 // We've modified the graph, so we'll have to recompute local 00541 // constants like the number of diagonal entries on this process. 00542 haveLocalConstants_ = false; 00543 } 00544 00547 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00548 void 00549 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00550 allocateIndices (ELocalGlobal lg) 00551 { 00552 // This is a protected function, only callable by us. If it was 00553 // called incorrectly, it is our fault. That's why the tests 00554 // below throw std::logic_error instead of std::invalid_argument. 00555 const char tfecfFuncName[] = "allocateIndices: "; 00556 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00557 isLocallyIndexed () && lg == GlobalIndices, std::logic_error, 00558 "The graph is locally indexed, but Tpetra code is calling this method " 00559 "with lg=GlobalIndices. Please report this bug to the Tpetra developers."); 00560 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00561 isGloballyIndexed () && lg == LocalIndices, std::logic_error, 00562 "The graph is globally indexed, but Tpetra code is calling this method " 00563 "with lg=LocalIndices. Please report this bug to the Tpetra developers."); 00564 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00565 indicesAreAllocated (), std::logic_error, "The graph's indices are " 00566 "already allocated, but Tpetra code is calling allocateIndices() again. " 00567 "Please report this bug to the Tpetra developers."); 00568 00569 const size_t numRows = getNodeNumRows(); 00570 indicesAreLocal_ = (lg == LocalIndices); 00571 indicesAreGlobal_ = (lg == GlobalIndices); 00572 nodeNumAllocated_ = 0; 00573 if (numAllocPerRow_ == null && getNodeNumRows() > 0) { 00574 // this wastes memory, temporarily, but it simplifies the code 00575 // and interfaces to follow 00576 ArrayRCP<size_t> tmpnumallocperrow = arcp<size_t>(numRows); 00577 std::fill(tmpnumallocperrow.begin(), tmpnumallocperrow.end(), numAllocForAllRows_); 00578 numAllocPerRow_ = tmpnumallocperrow; 00579 } 00580 // 00581 if (getProfileType() == StaticProfile) { 00582 // 00583 // STATIC ALLOCATION PROFILE 00584 // 00585 // Have the local sparse kernels object allocate row offsets for 00586 // us, with first-touch allocation if applicable. This is not 00587 // as important for global indices, because we never use global 00588 // indices in sparse kernels, but we might as well use the code 00589 // that we have for both the local and global indices cases. 00590 // Furthermore, first-touch allocation ensures that we don't 00591 // take up too much memory in any one NUMA affinity region. 00592 rowPtrs_ = sparse_ops_type::allocRowPtrs( getRowMap()->getNode(), numAllocPerRow_() ); 00593 if (lg == LocalIndices) { 00594 lclInds1D_ = sparse_ops_type::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), rowPtrs_() ); 00595 } 00596 else { 00597 gblInds1D_ = sparse_ops_type::template allocStorage<GlobalOrdinal>( getRowMap()->getNode(), rowPtrs_() ); 00598 } 00599 nodeNumAllocated_ = rowPtrs_[numRows]; 00600 } 00601 else { 00602 // 00603 // DYNAMIC ALLOCATION PROFILE 00604 // 00605 typename ArrayRCP<const size_t>::iterator numalloc = numAllocPerRow_.begin(); 00606 size_t howmany = numAllocForAllRows_; 00607 if (lg == LocalIndices) { 00608 lclInds2D_ = arcp< Array<LocalOrdinal> >(numRows); 00609 nodeNumAllocated_ = 0; 00610 for (size_t i=0; i < numRows; ++i) { 00611 if (numAllocPerRow_ != null) howmany = *numalloc++; 00612 nodeNumAllocated_ += howmany; 00613 if (howmany > 0) lclInds2D_[i].resize(howmany); 00614 } 00615 } 00616 else { // allocate global indices 00617 gblInds2D_ = arcp< Array<GlobalOrdinal> >(numRows); 00618 nodeNumAllocated_ = 0; 00619 for (size_t i=0; i < numRows; ++i) { 00620 if (numAllocPerRow_ != null) howmany = *numalloc++; 00621 nodeNumAllocated_ += howmany; 00622 if (howmany > 0) gblInds2D_[i].resize(howmany); 00623 } 00624 } 00625 } 00626 if (numRows > 0) { 00627 numRowEntries_ = arcp<size_t>(numRows); 00628 std::fill( numRowEntries_.begin(), numRowEntries_.end(), 0 ); 00629 } 00630 // done with these 00631 numAllocForAllRows_ = 0; 00632 numAllocPerRow_ = null; 00633 indicesAreAllocated_ = true; 00634 checkInternalState(); 00635 } 00636 00637 00640 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00641 template <class T> 00642 ArrayRCP<T> CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::allocateValues1D() const 00643 { 00644 const char tfecfFuncName[] = "allocateValues()"; 00645 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00646 indicesAreAllocated() == false, 00647 std::runtime_error, ": graph indices must be allocated before values." 00648 ); 00649 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00650 getProfileType() != StaticProfile, 00651 std::runtime_error, ": graph indices must be allocated in a static profile." 00652 ); 00653 ArrayRCP<T> values1D; 00654 values1D = sparse_ops_type::template allocStorage<T>( getRowMap()->getNode(), rowPtrs_() ); 00655 return values1D; 00656 } 00657 00658 00661 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00662 template <class T> 00663 ArrayRCP<Array<T> > 00664 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::allocateValues2D() const 00665 { 00666 const char tfecfFuncName[] = "allocateValues2D()"; 00667 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00668 indicesAreAllocated() == false, 00669 std::runtime_error, ": graph indices must be allocated before values." 00670 ); 00671 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00672 getProfileType() != DynamicProfile, 00673 std::runtime_error, ": graph indices must be allocated in a dynamic profile." 00674 ); 00675 ArrayRCP<Array<T> > values2D; 00676 values2D = arcp<Array<T> > (getNodeNumRows ()); 00677 if (lclInds2D_ != null) { 00678 const size_t numRows = lclInds2D_.size (); 00679 for (size_t r = 0; r < numRows; ++r) { 00680 values2D[r].resize (lclInds2D_[r].size ()); 00681 } 00682 } 00683 else if (gblInds2D_ != null) { 00684 const size_t numRows = gblInds2D_.size (); 00685 for (size_t r = 0; r < numRows; ++r) { 00686 values2D[r].resize (gblInds2D_[r].size ()); 00687 } 00688 } 00689 return values2D; 00690 } 00691 00692 00693 // ///////////////////////////////////////////////////////////////////////////// 00694 // ///////////////////////////////////////////////////////////////////////////// 00695 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 00696 // template <ELocalGlobal lg, class T> 00697 // RowInfo 00698 // CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00699 // updateAllocAndValues (RowInfo rowinfo, 00700 // size_t newAllocSize, 00701 // Array<T>& rowVals) 00702 // { 00703 // #ifdef HAVE_TPETRA_DEBUG 00704 // TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement(rowinfo.localRow) ); 00705 // TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize ); 00706 // TEUCHOS_TEST_FOR_EXCEPT( (lg == LocalIndices && ! isLocallyIndexed()) || 00707 // (lg == GlobalIndices && ! isGloballyIndexed()) ); 00708 // TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 ); 00709 // TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated() ); 00710 // #endif 00711 // // ArrayRCP::resize automatically copies over values on reallocation. 00712 // if (lg == LocalIndices) { 00713 // lclInds2D_[rowinfo.localRow].resize (newAllocSize); 00714 // } 00715 // else { // lg == GlobalIndices 00716 // gblInds2D_[rowinfo.localRow].resize (newAllocSize); 00717 // } 00718 // rowVals.resize (newAllocSize); 00719 // nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize); 00720 // rowinfo.allocSize = newAllocSize; 00721 // return rowinfo; 00722 // } 00723 00724 00727 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00728 ArrayView<const LocalOrdinal> 00729 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00730 getLocalView (RowInfo rowinfo) const 00731 { 00732 ArrayView<const LocalOrdinal> view; 00733 if (rowinfo.allocSize > 0) { 00734 if (lclInds1D_ != null) { 00735 view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize); 00736 } 00737 else if (! lclInds2D_[rowinfo.localRow].empty ()) { 00738 view = lclInds2D_[rowinfo.localRow] (); 00739 } 00740 } 00741 return view; 00742 } 00743 00744 00747 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00748 ArrayView<LocalOrdinal> 00749 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00750 getLocalViewNonConst (RowInfo rowinfo) 00751 { 00752 ArrayView<LocalOrdinal> view; 00753 if (rowinfo.allocSize > 0) { 00754 if (lclInds1D_ != null) { 00755 view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize); 00756 } 00757 else if (! lclInds2D_[rowinfo.localRow].empty()) { 00758 view = lclInds2D_[rowinfo.localRow] (); 00759 } 00760 } 00761 return view; 00762 } 00763 00764 00767 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00768 ArrayView<const GlobalOrdinal> 00769 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00770 getGlobalView (RowInfo rowinfo) const 00771 { 00772 ArrayView<const GlobalOrdinal> view; 00773 if (rowinfo.allocSize > 0) { 00774 if (gblInds1D_ != null) { 00775 view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize); 00776 } 00777 else if (! gblInds2D_[rowinfo.localRow].empty()) { 00778 view = gblInds2D_[rowinfo.localRow] (); 00779 } 00780 } 00781 return view; 00782 } 00783 00784 00787 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00788 ArrayView<GlobalOrdinal> 00789 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00790 getGlobalViewNonConst (RowInfo rowinfo) 00791 { 00792 ArrayView<GlobalOrdinal> view; 00793 if (rowinfo.allocSize > 0) { 00794 if (gblInds1D_ != null) { 00795 view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize); 00796 } 00797 else if (!gblInds2D_[rowinfo.localRow].empty()) { 00798 view = gblInds2D_[rowinfo.localRow] (); 00799 } 00800 } 00801 return view; 00802 } 00803 00804 00807 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00808 RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRowInfo(size_t myRow) const 00809 { 00810 #ifdef HAVE_TPETRA_DEBUG 00811 const char tfecfFuncName[] = "getRowInfo"; 00812 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00813 rowMap_->isNodeLocalElement (myRow) == false, 00814 std::logic_error, 00815 ": The given (local) row index myRow = " << myRow 00816 << " does not belong to the graph's row Map. " 00817 "This probably indicates a bug in Tpetra::CrsGraph or Tpetra::CrsMatrix. " 00818 "Please report this to the Tpetra developers."); 00819 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00820 hasRowInfo() == false, std::logic_error, 00821 ": Late catch! Graph does not have row info anymore. " 00822 "Error should have been caught earlier. Please contact Tpetra team."); 00823 #endif // HAVE_TPETRA_DEBUG 00824 const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid (); 00825 RowInfo ret; 00826 ret.localRow = myRow; 00827 if (nodeNumAllocated_ != 0 && nodeNumAllocated_ != STINV) { 00828 // graph data structures have the info that we need 00829 // 00830 // if static graph, offsets tell us the allocation size 00831 if (getProfileType() == StaticProfile) { 00832 ret.offset1D = rowPtrs_[myRow]; 00833 ret.allocSize = rowPtrs_[myRow+1] - rowPtrs_[myRow]; 00834 if (numRowEntries_ == null) { 00835 ret.numEntries = ret.allocSize; 00836 } 00837 else { 00838 ret.numEntries = numRowEntries_[myRow]; 00839 } 00840 } 00841 else { 00842 ret.offset1D = STINV; 00843 if (isLocallyIndexed ()) { 00844 ret.allocSize = lclInds2D_[myRow].size(); 00845 } 00846 else { 00847 ret.allocSize = gblInds2D_[myRow].size(); 00848 } 00849 ret.numEntries = numRowEntries_[myRow]; 00850 } 00851 } 00852 else if (nodeNumAllocated_ == 0) { 00853 // have performed allocation, but the graph has no allocation or entries 00854 ret.allocSize = 0; 00855 ret.numEntries = 0; 00856 ret.offset1D = STINV; 00857 } 00858 else if (indicesAreAllocated () == false) { 00859 // haven't performed allocation yet; probably won't hit this code 00860 if (numAllocPerRow_ == null) { 00861 ret.allocSize = numAllocForAllRows_; 00862 } 00863 else { 00864 ret.allocSize = numAllocPerRow_[myRow]; 00865 } 00866 ret.numEntries = 0; 00867 ret.offset1D = STINV; 00868 } 00869 else { 00870 // don't know how we ended up here... 00871 TEUCHOS_TEST_FOR_EXCEPT(true); 00872 } 00873 return ret; 00874 } 00875 00876 00879 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00880 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::staticAssertions() const 00881 { 00882 // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal): 00883 // This is so that we can store local indices in the memory 00884 // formerly occupied by global indices. 00885 // 00886 // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal) and 00887 // max(size_t) >= max(LocalOrdinal) 00888 // This is so that we can represent any LocalOrdinal as a 00889 // size_t, and any LocalOrdinal as a GlobalOrdinal 00890 Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size1; (void)cta_size1; 00891 Teuchos::CompileTimeAssert<sizeof(global_size_t) < sizeof(size_t) > cta_size2; (void)cta_size2; 00892 // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check 00893 std::string msg = typeName(*this) + ": Object cannot be allocated with stated template arguments: size assumptions are not valid."; 00894 TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(), std::runtime_error, msg); 00895 TEUCHOS_TEST_FOR_EXCEPTION( (global_size_t)OrdinalTraits<LocalOrdinal>::max() > (global_size_t)OrdinalTraits<GlobalOrdinal>::max(), std::runtime_error, msg); 00896 TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(), std::runtime_error, msg); 00897 TEUCHOS_TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(), std::runtime_error, msg); 00898 } 00899 00900 00903 template <class LocalOrdinal, class GlobalOrdinal, class Node> 00904 template <ELocalGlobal lg> 00905 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::filterIndices(const SLocalGlobalNCViews &inds) const 00906 { 00907 const map_type& cmap = *colMap_; 00908 Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; 00909 (void)cta_lg; 00910 size_t numFiltered = 0; 00911 #ifdef HAVE_TPETRA_DEBUG 00912 size_t numFiltered_debug = 0; 00913 #endif 00914 if (lg == GlobalIndices) { 00915 ArrayView<GlobalOrdinal> ginds = inds.ginds; 00916 typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(); 00917 typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin(); 00918 while (cptr != ginds.end()) { 00919 if (cmap.isNodeGlobalElement(*cptr)) { 00920 *fend++ = *cptr; 00921 #ifdef HAVE_TPETRA_DEBUG 00922 ++numFiltered_debug; 00923 #endif 00924 } 00925 ++cptr; 00926 } 00927 numFiltered = fend - ginds.begin(); 00928 } 00929 else if (lg == LocalIndices) { 00930 ArrayView<LocalOrdinal> linds = inds.linds; 00931 typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(); 00932 typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin(); 00933 while (cptr != linds.end()) { 00934 if (cmap.isNodeLocalElement(*cptr)) { 00935 *fend++ = *cptr; 00936 #ifdef HAVE_TPETRA_DEBUG 00937 ++numFiltered_debug; 00938 #endif 00939 } 00940 ++cptr; 00941 } 00942 numFiltered = fend - linds.begin(); 00943 } 00944 #ifdef HAVE_TPETRA_DEBUG 00945 TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug ); 00946 #endif 00947 return numFiltered; 00948 } 00949 00950 00951 // ///////////////////////////////////////////////////////////////////////////// 00952 // ///////////////////////////////////////////////////////////////////////////// 00953 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 00954 // template <class T> 00955 // size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00956 // filterGlobalIndicesAndValues (const ArrayView<GlobalOrdinal>& ginds, 00957 // const ArrayView<T>& vals) const 00958 // { 00959 // const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_; 00960 // size_t numFiltered = 0; 00961 // typename ArrayView<T>::iterator fvalsend = vals.begin(); 00962 // typename ArrayView<T>::iterator valscptr = vals.begin(); 00963 // #ifdef HAVE_TPETRA_DEBUG 00964 // size_t numFiltered_debug = 0; 00965 // #endif 00966 // typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(); 00967 // typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin(); 00968 // while (cptr != ginds.end()) { 00969 // if (cmap.isNodeGlobalElement (*cptr)) { 00970 // *fend++ = *cptr; 00971 // *fvalsend++ = *valscptr; 00972 // #ifdef HAVE_TPETRA_DEBUG 00973 // ++numFiltered_debug; 00974 // #endif 00975 // } 00976 // ++cptr; 00977 // ++valscptr; 00978 // } 00979 // numFiltered = fend - ginds.begin(); 00980 // #ifdef HAVE_TPETRA_DEBUG 00981 // TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug ); 00982 // TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() ); 00983 // const size_t numFilteredActual = 00984 // Teuchos::as<size_t> (fvalsend - vals.begin ()); 00985 // TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual ); 00986 // #endif 00987 // return numFiltered; 00988 // } 00989 00990 00991 // ///////////////////////////////////////////////////////////////////////////// 00992 // ///////////////////////////////////////////////////////////////////////////// 00993 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 00994 // template <class T> 00995 // size_t 00996 // CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 00997 // filterLocalIndicesAndValues (const ArrayView<LocalOrdinal>& linds, 00998 // const ArrayView<T>& vals) const 00999 // { 01000 // const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_; 01001 // size_t numFiltered = 0; 01002 // typename ArrayView<T>::iterator fvalsend = vals.begin(); 01003 // typename ArrayView<T>::iterator valscptr = vals.begin(); 01004 // #ifdef HAVE_TPETRA_DEBUG 01005 // size_t numFiltered_debug = 0; 01006 // #endif 01007 // typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(); 01008 // typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin(); 01009 // while (cptr != linds.end()) { 01010 // if (cmap.isNodeLocalElement (*cptr)) { 01011 // *fend++ = *cptr; 01012 // *fvalsend++ = *valscptr; 01013 // #ifdef HAVE_TPETRA_DEBUG 01014 // ++numFiltered_debug; 01015 // #endif 01016 // } 01017 // ++cptr; 01018 // ++valscptr; 01019 // } 01020 // numFiltered = fend - linds.begin(); 01021 // #ifdef HAVE_TPETRA_DEBUG 01022 // TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug ); 01023 // TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() ); 01024 // const size_t numFilteredActual = 01025 // Teuchos::as<size_t> (fvalsend - vals.begin ()); 01026 // TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual ); 01027 // #endif 01028 // return numFiltered; 01029 // } 01030 01031 01034 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01035 size_t 01036 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01037 insertIndices (const RowInfo& rowinfo, 01038 const SLocalGlobalViews &newInds, 01039 const ELocalGlobal lg, 01040 const ELocalGlobal I) 01041 { 01042 #ifdef HAVE_TPETRA_DEBUG 01043 TEUCHOS_TEST_FOR_EXCEPTION( 01044 lg != GlobalIndices && lg != LocalIndices, std::invalid_argument, 01045 "Tpetra::CrsGraph::insertIndices: lg must be either GlobalIndices or " 01046 "LocalIndices."); 01047 #endif // HAVE_TPETRA_DEBUG 01048 size_t numNewInds = 0; 01049 if (lg == GlobalIndices) { // input indices are global 01050 ArrayView<const GlobalOrdinal> new_ginds = newInds.ginds; 01051 numNewInds = new_ginds.size(); 01052 if (I == GlobalIndices) { // store global indices 01053 ArrayView<GlobalOrdinal> gind_view = getGlobalViewNonConst(rowinfo); 01054 std::copy(new_ginds.begin(), new_ginds.end(), gind_view.begin()+rowinfo.numEntries); 01055 } 01056 else if (I == LocalIndices) { // store local indices 01057 ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo); 01058 typename ArrayView<const GlobalOrdinal>::iterator in = new_ginds.begin(); 01059 const typename ArrayView<const GlobalOrdinal>::iterator stop = new_ginds.end(); 01060 typename ArrayView<LocalOrdinal>::iterator out = lind_view.begin()+rowinfo.numEntries; 01061 while (in != stop) { 01062 *out++ = colMap_->getLocalElement (*in++); 01063 } 01064 } 01065 } 01066 else if (lg == LocalIndices) { // input indices are local 01067 ArrayView<const LocalOrdinal> new_linds = newInds.linds; 01068 numNewInds = new_linds.size(); 01069 if (I == LocalIndices) { // store local indices 01070 ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo); 01071 std::copy(new_linds.begin(), new_linds.end(), lind_view.begin()+rowinfo.numEntries); 01072 } 01073 else if (I == GlobalIndices) { 01074 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Tpetra::CrsGraph::" 01075 "insertIndices: the case where the input indices are local and the " 01076 "indices to write are global (lg=LocalIndices, I=GlobalIndices) is " 01077 "not implemented, because it does not make sense." << std::endl << 01078 "If you have correct local column indices, that means the graph has " 01079 "a column Map. In that case, you should be storing local indices."); 01080 } 01081 } 01082 numRowEntries_[rowinfo.localRow] += numNewInds; 01083 nodeNumEntries_ += numNewInds; 01084 setLocallyModified (); 01085 return numNewInds; 01086 } 01087 01090 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01091 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01092 insertGlobalIndicesImpl (const LocalOrdinal myRow, 01093 const ArrayView<const GlobalOrdinal> &indices) 01094 { 01095 const char* tfecfFuncName ("insertGlobalIndicesImpl"); 01096 01097 RowInfo rowInfo = getRowInfo(myRow); 01098 const size_t numNewInds = indices.size(); 01099 const size_t newNumEntries = rowInfo.numEntries + numNewInds; 01100 if (newNumEntries > rowInfo.allocSize) { 01101 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01102 getProfileType() == StaticProfile, std::runtime_error, 01103 ": new indices exceed statically allocated graph structure."); 01104 01105 // update allocation, doubling size to reduce number of reallocations 01106 size_t newAllocSize = 2*rowInfo.allocSize; 01107 if (newAllocSize < newNumEntries) 01108 newAllocSize = newNumEntries; 01109 gblInds2D_[myRow].resize(newAllocSize); 01110 nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize); 01111 } 01112 01113 // Copy new indices at end of global index array 01114 if (gblInds1D_ != null) 01115 std::copy(indices.begin(), indices.end(), 01116 gblInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries); 01117 else 01118 std::copy(indices.begin(), indices.end(), 01119 gblInds2D_[myRow].begin()+rowInfo.numEntries); 01120 numRowEntries_[myRow] += numNewInds; 01121 nodeNumEntries_ += numNewInds; 01122 setLocallyModified (); 01123 01124 #ifdef HAVE_TPETRA_DEBUG 01125 { 01126 const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow); 01127 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01128 chkNewNumEntries != newNumEntries, std::logic_error, 01129 ": Internal logic error. Please contact Tpetra team."); 01130 } 01131 #endif 01132 } 01133 01136 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01137 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01138 insertLocalIndicesImpl (const LocalOrdinal myRow, 01139 const ArrayView<const LocalOrdinal> &indices) 01140 { 01141 const char* tfecfFuncName ("insertLocallIndicesImpl"); 01142 01143 RowInfo rowInfo = getRowInfo(myRow); 01144 const size_t numNewInds = indices.size(); 01145 const size_t newNumEntries = rowInfo.numEntries + numNewInds; 01146 if (newNumEntries > rowInfo.allocSize) { 01147 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01148 getProfileType() == StaticProfile, std::runtime_error, 01149 ": new indices exceed statically allocated graph structure."); 01150 01151 // update allocation, doubling size to reduce number of reallocations 01152 size_t newAllocSize = 2*rowInfo.allocSize; 01153 if (newAllocSize < newNumEntries) 01154 newAllocSize = newNumEntries; 01155 lclInds2D_[myRow].resize(newAllocSize); 01156 nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize); 01157 } 01158 01159 // Insert new indices at end of lclInds array 01160 if (lclInds1D_ != null) 01161 std::copy(indices.begin(), indices.end(), 01162 lclInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries); 01163 else 01164 std::copy(indices.begin(), indices.end(), 01165 lclInds2D_[myRow].begin()+rowInfo.numEntries); 01166 numRowEntries_[myRow] += numNewInds; 01167 nodeNumEntries_ += numNewInds; 01168 setLocallyModified (); 01169 #ifdef HAVE_TPETRA_DEBUG 01170 { 01171 const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow); 01172 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01173 chkNewNumEntries != newNumEntries, std::logic_error, 01174 ": Internal logic error. Please contact Tpetra team."); 01175 } 01176 #endif 01177 } 01178 01179 // ///////////////////////////////////////////////////////////////////////////// 01180 // ///////////////////////////////////////////////////////////////////////////// 01181 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 01182 // template <class Scalar> 01183 // void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01184 // insertIndicesAndValues (const RowInfo& rowInfo, 01185 // const SLocalGlobalViews& newInds, 01186 // const ArrayView<Scalar>& oldRowVals, 01187 // const ArrayView<const Scalar>& newRowVals, 01188 // const ELocalGlobal lg, 01189 // const ELocalGlobal I) 01190 // { 01191 // const size_t numNewInds = insertIndices (rowInfo, newInds, lg, I); 01192 // typename ArrayView<const Scalar>::const_iterator newRowValsBegin = 01193 // newRowVals.begin (); 01194 // std::copy (newRowValsBegin, newRowValsBegin + numNewInds, 01195 // oldRowVals.begin () + rowInfo.numEntries); 01196 // } 01197 01198 // ///////////////////////////////////////////////////////////////////////////// 01199 // ///////////////////////////////////////////////////////////////////////////// 01200 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 01201 // template <class Scalar, class BinaryFunction> 01202 // void 01203 // CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01204 // transformLocalValues (RowInfo rowInfo, 01205 // const Teuchos::ArrayView<Scalar>& rowVals, 01206 // const Teuchos::ArrayView<const LocalOrdinal>& inds, 01207 // const Teuchos::ArrayView<const Scalar>& newVals, 01208 // BinaryFunction f) const 01209 // { 01210 // const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid(); 01211 // const size_t numElts = Teuchos::as<size_t> (inds.size ()); 01212 // size_t hint = 0; // Guess for the current index k into rowVals 01213 01214 // // Get a view of the column indices in the row. This amortizes 01215 // // the cost of getting the view over all the entries of inds. 01216 // ArrayView<const LocalOrdinal> colInds = getLocalView (rowInfo); 01217 01218 // for (size_t j = 0; j < numElts; ++j) { 01219 // const size_t k = findLocalIndex (rowInfo, inds[j], colInds, hint); 01220 // if (k != STINV) { 01221 // rowVals[k] = f( rowVals[k], newVals[j] ); 01222 // hint = k+1; 01223 // } 01224 // } 01225 // } 01226 01227 01228 // ///////////////////////////////////////////////////////////////////////////// 01229 // ///////////////////////////////////////////////////////////////////////////// 01230 // template <class LocalOrdinal, class GlobalOrdinal, class Node> 01231 // template <class Scalar, class BinaryFunction> 01232 // void 01233 // CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01234 // transformGlobalValues (RowInfo rowInfo, 01235 // const Teuchos::ArrayView<Scalar>& rowVals, 01236 // const Teuchos::ArrayView<const GlobalOrdinal>& inds, 01237 // const Teuchos::ArrayView<const Scalar>& newVals, 01238 // BinaryFunction f) const 01239 // { 01240 // const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid(); 01241 // const size_t numElts = Teuchos::as<size_t> (inds.size ()); 01242 // size_t hint = 0; // hint is a guess as to wheter the index is 01243 01244 // for (size_t j = 0; j < numElts; ++j) { 01245 // const size_t k = findGlobalIndex (rowInfo, inds[j], hint); 01246 // if (k != STINV) { 01247 // rowVals[k] = f( rowVals[k], newVals[j] ); 01248 // hint = k+1; 01249 // } 01250 // } 01251 // } 01252 01253 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01254 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01255 sortRowIndices (RowInfo rowinfo) 01256 { 01257 if (rowinfo.numEntries > 0) { 01258 Teuchos::ArrayView<LocalOrdinal> inds_view = 01259 this->getLocalViewNonConst (rowinfo); 01260 std::sort (inds_view.begin (), inds_view.begin () + rowinfo.numEntries); 01261 } 01262 } 01263 01264 // in the future, perhaps this could use std::sort with a boost::zip_iterator 01265 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01266 template <class Scalar> 01267 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01268 sortRowIndicesAndValues (const RowInfo rowinfo, 01269 const Teuchos::ArrayView<Scalar>& values) 01270 { 01271 if (rowinfo.numEntries > 0) { 01272 Teuchos::ArrayView<LocalOrdinal> indsView = 01273 this->getLocalViewNonConst (rowinfo); 01274 sort2 (indsView.begin (), indsView.begin () + rowinfo.numEntries, 01275 values.begin ()); 01276 } 01277 } 01278 01279 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01280 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01281 mergeRowIndices (RowInfo rowinfo) 01282 { 01283 using Teuchos::ArrayView; 01284 const char tfecfFuncName[] = "mergeRowIndices: "; 01285 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01286 isStorageOptimized (), std::logic_error, "The graph is already storage " 01287 "optimized, so we shouldn't be merging any indices. " 01288 "Please report this bug to the Tpetra developers."); 01289 01290 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo); 01291 typename ArrayView<LocalOrdinal>::iterator beg, end, newend; 01292 beg = inds_view.begin(); 01293 end = inds_view.begin() + rowinfo.numEntries; 01294 newend = std::unique(beg,end); 01295 const size_t mergedEntries = newend - beg; 01296 #ifdef HAVE_TPETRA_DEBUG 01297 // merge should not have eliminated any entries; if so, the 01298 // assignment below will destroy the packed structure 01299 TEUCHOS_TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries ); 01300 #endif 01301 numRowEntries_[rowinfo.localRow] = mergedEntries; 01302 nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries); 01303 } 01304 01305 // in the future, this could use std::unique with a boost::zip_iterator 01306 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01307 template<class Scalar> 01308 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01309 mergeRowIndicesAndValues (RowInfo rowinfo, const ArrayView<Scalar>& rowValues) 01310 { 01311 const char tfecfFuncName[] = "mergeRowIndicesAndValues"; 01312 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01313 isStorageOptimized(), std::logic_error, ": It is invalid to call this " 01314 "method if the graph's storage has already been optimized." << std::endl 01315 << "Please report this bug to the Tpetra developers."); 01316 01317 typedef typename ArrayView<Scalar>::iterator Iter; 01318 Iter rowValueIter = rowValues.begin (); 01319 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst (rowinfo); 01320 typename ArrayView<LocalOrdinal>::iterator beg, end, newend; 01321 01322 // beg,end define a half-exclusive interval over which to iterate. 01323 beg = inds_view.begin(); 01324 end = inds_view.begin() + rowinfo.numEntries; 01325 newend = beg; 01326 if (beg != end) { 01327 typename ArrayView<LocalOrdinal>::iterator cur = beg + 1; 01328 Iter vcur = rowValueIter + 1; 01329 Iter vend = rowValueIter; 01330 cur = beg+1; 01331 while (cur != end) { 01332 if (*cur != *newend) { 01333 // new entry; save it 01334 ++newend; 01335 ++vend; 01336 (*newend) = (*cur); 01337 (*vend) = (*vcur); 01338 } 01339 else { 01340 // old entry; merge it 01341 //(*vend) = f (*vend, *vcur); 01342 (*vend) += *vcur; 01343 } 01344 ++cur; 01345 ++vcur; 01346 } 01347 ++newend; // one past the last entry, per typical [beg,end) semantics 01348 } 01349 const size_t mergedEntries = newend - beg; 01350 #ifdef HAVE_TPETRA_DEBUG 01351 // merge should not have eliminated any entries; if so, the 01352 // assignment below will destroy the packed structure 01353 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01354 isStorageOptimized() && mergedEntries != rowinfo.numEntries, 01355 std::logic_error, 01356 ": Merge was incorrect; it eliminated entries from the graph. " 01357 << std::endl << "Please report this bug to the Tpetra developers."); 01358 #endif // HAVE_TPETRA_DEBUG 01359 numRowEntries_[rowinfo.localRow] = mergedEntries; 01360 nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries); 01361 } 01362 01363 01366 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01367 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01368 setDomainRangeMaps (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& domainMap, 01369 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rangeMap) 01370 { 01371 // simple pointer comparison for equality 01372 if (domainMap_ != domainMap) { 01373 domainMap_ = domainMap; 01374 importer_ = null; 01375 } 01376 if (rangeMap_ != rangeMap) { 01377 rangeMap_ = rangeMap; 01378 exporter_ = null; 01379 } 01380 } 01381 01382 01385 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01386 size_t 01387 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01388 findLocalIndex (RowInfo rowinfo, LocalOrdinal ind, size_t hint) const 01389 { 01390 ArrayView<const LocalOrdinal> colInds = getLocalView (rowinfo); 01391 return this->findLocalIndex (rowinfo, ind, colInds, hint); 01392 } 01393 01394 01397 template <class LocalOrdinal, 01398 class GlobalOrdinal, 01399 class Node> 01400 size_t 01401 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01402 findLocalIndex (RowInfo rowinfo, 01403 LocalOrdinal ind, 01404 ArrayView<const LocalOrdinal> colInds, 01405 size_t hint) const 01406 { 01407 typedef typename ArrayView<const LocalOrdinal>::iterator IT; 01408 01409 // If the hint was correct, then the hint is the offset to return. 01410 if (hint < rowinfo.numEntries && colInds[hint] == ind) { 01411 return hint; 01412 } 01413 01414 // The hint was wrong, so we must search for the given column 01415 // index in the column indices for the given row. How we do the 01416 // search depends on whether the graph's column indices are 01417 // sorted. 01418 IT beg = colInds.begin (); 01419 IT end = beg + rowinfo.numEntries; 01420 IT ptr = beg + rowinfo.numEntries; // "null" 01421 bool found = true; 01422 01423 if (isSorted ()) { 01424 // binary search 01425 std::pair<IT,IT> p = std::equal_range (beg, end, ind); 01426 if (p.first == p.second) { 01427 found = false; 01428 } 01429 else { 01430 ptr = p.first; 01431 } 01432 } 01433 else { 01434 // direct search 01435 ptr = std::find (beg, end, ind); 01436 if (ptr == end) { 01437 found = false; 01438 } 01439 } 01440 01441 if (found) { 01442 return Teuchos::as<size_t> (ptr - beg); 01443 } 01444 else { 01445 return Teuchos::OrdinalTraits<size_t>::invalid (); 01446 } 01447 } 01448 01449 01452 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01453 size_t 01454 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01455 findGlobalIndex (RowInfo rowinfo, GlobalOrdinal ind, size_t hint) const 01456 { 01457 using Teuchos::ArrayView; 01458 typedef typename ArrayView<const GlobalOrdinal>::iterator IT; 01459 01460 // Don't let an invalid global column index through. 01461 if (ind == Teuchos::OrdinalTraits<GlobalOrdinal>::invalid ()) { 01462 return Teuchos::OrdinalTraits<size_t>::invalid (); 01463 } 01464 01465 ArrayView<const GlobalOrdinal> indices = getGlobalView (rowinfo); 01466 01467 // We don't actually require that the hint be a valid index. 01468 // If it is not in range, we just ignore it. 01469 if (hint < rowinfo.numEntries && indices[hint] == ind) { 01470 return hint; 01471 } 01472 01473 IT beg = indices.begin (); 01474 IT end = indices.begin () + rowinfo.numEntries; // not indices.end() 01475 if (isSorted ()) { // use binary search 01476 const std::pair<IT,IT> p = std::equal_range (beg, end, ind); 01477 if (p.first == p.second) { // range of matching entries is empty 01478 return Teuchos::OrdinalTraits<size_t>::invalid (); 01479 } else { 01480 return p.first - beg; 01481 } 01482 } 01483 else { // not sorted; must use linear search 01484 const IT loc = std::find (beg, end, ind); 01485 if (loc == end) { 01486 return Teuchos::OrdinalTraits<size_t>::invalid (); 01487 } else { 01488 return loc - beg; 01489 } 01490 } 01491 } 01492 01493 01496 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01497 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::clearGlobalConstants() 01498 { 01499 globalNumEntries_ = OrdinalTraits<global_size_t>::invalid(); 01500 globalNumDiags_ = OrdinalTraits<global_size_t>::invalid(); 01501 globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid(); 01502 haveGlobalConstants_ = false; 01503 } 01504 01505 01508 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01509 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01510 checkInternalState () const 01511 { 01512 #ifdef HAVE_TPETRA_DEBUG 01513 const global_size_t GSTI = OrdinalTraits<global_size_t>::invalid(); 01514 const size_t STI = OrdinalTraits<size_t>::invalid(); 01515 std::string err = typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team."; 01516 // check the internal state of this data structure 01517 // this is called by numerous state-changing methods, in a debug build, to ensure that the object 01518 // always remains in a valid state 01519 // the graph should have been allocated with a row map 01520 TEUCHOS_TEST_FOR_EXCEPTION( rowMap_ == null, std::logic_error, err ); 01521 // am either complete or active 01522 TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() == isFillComplete(), std::logic_error, err ); 01523 // if active, i have no local graph 01524 TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() && lclGraph_ != null, std::logic_error, err ); 01525 // if the graph has been fill completed, then all maps should be present 01526 TEUCHOS_TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err ); 01527 // if storage has been optimized, then indices should have been allocated (even if trivially so) 01528 TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err ); 01529 // if storage has been optimized, then number of allocated is now the number of entries 01530 TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && getNodeAllocationSize() != getNodeNumEntries(), std::logic_error, err ); 01531 // if graph doesn't have the global constants, then they should all be marked as invalid 01532 TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != GSTI || globalNumDiags_ != GSTI || globalMaxNumRowEntries_ != GSTI ), std::logic_error, err ); 01533 // if the graph has global cosntants, then they should be valid. 01534 TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == GSTI || globalNumDiags_ == GSTI || globalMaxNumRowEntries_ == GSTI ), std::logic_error, err ); 01535 TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ), 01536 std::logic_error, err ); 01537 // if indices are allocated, then the allocation specifications should have been released 01538 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null), std::logic_error, err ); 01539 // if indices are not allocated, then information dictating allocation quantities should be present 01540 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != STI || nodeNumEntries_ != 0), std::logic_error, err ); 01541 // if storage is optimized, then profile should be static 01542 TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile, std::logic_error, err ); 01543 // if rowPtrs_ exists, it is required to have N+1 rows, and rowPtrs_[N] == gblInds1D_.size()/lclInds1D_.size() 01544 TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)gblInds1D_.size()), std::logic_error, err ); 01545 TEUCHOS_TEST_FOR_EXCEPTION( isLocallyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)lclInds1D_.size()), std::logic_error, err ); 01546 // if profile is dynamic and we have allocated, then 2D allocations should be present 01547 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && lclInds2D_ == null && gblInds2D_ == null, 01548 std::logic_error, err ); 01549 // if profile is dynamic, then numentries and 2D indices are needed and should be present 01550 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (numRowEntries_ == null || (lclInds2D_ == null && gblInds2D_ == null)), 01551 std::logic_error, err ); 01552 // if profile is dynamic, then 1D allocations should not be present 01553 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (lclInds1D_ != null || gblInds1D_ != null), std::logic_error, err ); 01554 // if profile is dynamic, then row offsets should not be present 01555 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && rowPtrs_ != null, std::logic_error, err ); 01556 // if profile is static and we have allocated non-trivially, then 1D allocations should be present 01557 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeAllocationSize() > 0 && lclInds1D_ == null && gblInds1D_ == null, 01558 std::logic_error, err ); 01559 // if profile is static, then 2D allocations should not be present 01560 TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (lclInds2D_ != null || gblInds2D_ != null), std::logic_error, err ); 01561 // if indices are not allocated, then none of the buffers should be. 01562 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (rowPtrs_ != null || numRowEntries_ != null || 01563 lclInds1D_ != null || lclInds2D_ != null || 01564 gblInds1D_ != null || gblInds2D_ != null), std::logic_error, err ); 01565 // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_) 01566 TEUCHOS_TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false, std::logic_error, err ); 01567 // indices may be local or global, but not both 01568 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true, std::logic_error, err ); 01569 // if indices are local, then global allocations should not be present 01570 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (gblInds1D_ != null || gblInds2D_ != null), std::logic_error, err ); 01571 // if indices are global, then local allocations should not be present 01572 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (lclInds1D_ != null || lclInds2D_ != null), std::logic_error, err ); 01573 // if indices are local, then local allocations should be present 01574 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && getNodeAllocationSize() > 0 && lclInds1D_ == null && getNodeNumRows() > 0 && lclInds2D_ == null, 01575 std::logic_error, err ); 01576 // if indices are global, then global allocations should be present 01577 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && getNodeAllocationSize() > 0 && gblInds1D_ == null && getNodeNumRows() > 0 && gblInds2D_ == null, 01578 std::logic_error, err ); 01579 // if indices are allocated, then we should have recorded how many were allocated 01580 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true && nodeNumAllocated_ == STI, std::logic_error, err ); 01581 // if indices are not allocated, then the allocation size should be marked invalid 01582 TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != STI, std::logic_error, err ); 01583 // check the actual allocations 01584 if (indicesAreAllocated()) { 01585 size_t actualNumAllocated = 0; 01586 if (pftype_ == DynamicProfile) { 01587 if (isGloballyIndexed() && gblInds2D_ != null) { 01588 for (size_t r = 0; r < getNodeNumRows(); ++r) { 01589 actualNumAllocated += gblInds2D_[r].size(); 01590 } 01591 } 01592 else if (isLocallyIndexed() && lclInds2D_ != null) { 01593 for (size_t r = 0; r < getNodeNumRows(); ++r) { 01594 actualNumAllocated += lclInds2D_[r].size(); 01595 } 01596 } 01597 TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err ); 01598 } 01599 else if (rowPtrs_ != null) { // pftype_ == StaticProfile 01600 actualNumAllocated = rowPtrs_[getNodeNumRows()]; 01601 TEUCHOS_TEST_FOR_EXCEPTION( isLocallyIndexed() == true && (size_t)lclInds1D_.size() != actualNumAllocated, std::logic_error, err ); 01602 TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)gblInds1D_.size() != actualNumAllocated, std::logic_error, err ); 01603 TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err ); 01604 } 01605 } 01606 #endif 01607 } 01608 01609 01612 // // 01613 // User-visible class methods // 01614 // // 01617 01618 01621 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01622 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const 01623 { 01624 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01625 size_t ret = OrdinalTraits<size_t>::invalid(); 01626 if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid()) 01627 { 01628 RowInfo rowinfo = getRowInfo(lrow); 01629 ret = rowinfo.numEntries; 01630 } 01631 return ret; 01632 } 01633 01634 01637 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01638 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumEntriesInLocalRow(LocalOrdinal localRow) const 01639 { 01640 size_t ret = OrdinalTraits<size_t>::invalid(); 01641 if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) { 01642 RowInfo rowinfo = getRowInfo(localRow); 01643 ret = rowinfo.numEntries; 01644 } 01645 return ret; 01646 } 01647 01648 01651 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01652 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const 01653 { 01654 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01655 size_t ret = OrdinalTraits<size_t>::invalid(); 01656 if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid()) 01657 { 01658 RowInfo rowinfo = getRowInfo(lrow); 01659 ret = rowinfo.allocSize; 01660 } 01661 return ret; 01662 } 01663 01664 01667 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01668 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const 01669 { 01670 size_t ret = OrdinalTraits<size_t>::invalid(); 01671 if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) { 01672 RowInfo rowinfo = getRowInfo(localRow); 01673 ret = rowinfo.allocSize; 01674 } 01675 return ret; 01676 } 01677 01678 01681 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01682 ArrayRCP<const size_t> CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeRowPtrs() const 01683 { 01684 return rowPtrs_; 01685 } 01686 01687 01690 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01691 ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodePackedIndices() const 01692 { 01693 return lclInds1D_; 01694 } 01695 01696 01699 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01700 void 01701 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01702 getLocalRowCopy (LocalOrdinal localRow, 01703 const ArrayView<LocalOrdinal>& indices, 01704 size_t& numEntries) const 01705 { 01706 using Teuchos::ArrayView; 01707 typedef LocalOrdinal LO; 01708 typedef GlobalOrdinal GO; 01709 01710 TEUCHOS_TEST_FOR_EXCEPTION( 01711 isGloballyIndexed () && ! hasColMap (), std::runtime_error, 01712 "Tpetra::CrsGraph::getLocalRowCopy: The graph is globally indexed and " 01713 "does not have a column Map yet. That means we don't have local indices " 01714 "for columns yet, so it doesn't make sense to call this method. If the " 01715 "graph doesn't have a column Map yet, you should call fillComplete on " 01716 "it first."); 01717 TEUCHOS_TEST_FOR_EXCEPTION( 01718 ! hasRowInfo(), std::runtime_error, 01719 "Tpetra::CrsGraph::getLocalRowCopy: graph row information was deleted " 01720 "at fillComplete()."); 01721 01722 01723 if (! getRowMap ()->isNodeLocalElement (localRow)) { 01724 numEntries = 0; 01725 return; 01726 } 01727 01728 const RowInfo rowinfo = getRowInfo (localRow); 01729 const size_t theNumEntries = rowinfo.numEntries; 01730 01731 TEUCHOS_TEST_FOR_EXCEPTION( 01732 static_cast<size_t> (indices.size ()) < theNumEntries, 01733 std::runtime_error, 01734 "Tpetra::CrsGraph::getLocalRowCopy: The given row " << localRow << " has " 01735 << theNumEntries << " entries, but indices.size() = " << indices.size () 01736 << ", which does not suffice to store the row's indices."); 01737 01738 numEntries = theNumEntries; 01739 01740 if (isLocallyIndexed ()) { 01741 ArrayView<const LO> lview = getLocalView (rowinfo); 01742 std::copy (lview.begin (), lview.begin () + numEntries, indices.begin ()); 01743 } 01744 else if (isGloballyIndexed ()) { 01745 ArrayView<const GO> gview = getGlobalView (rowinfo); 01746 const map_type& colMap = * (this->getColMap ()); 01747 for (size_t j = 0; j < numEntries; ++j) { 01748 indices[j] = colMap.getLocalElement (gview[j]); 01749 } 01750 } 01751 else { 01752 // If the graph on the calling process is neither locally nor 01753 // globally indexed, that means it owns no column indices. 01754 // 01755 // FIXME (mfh 21 Oct 2013) It's not entirely clear to me whether 01756 // we can reach this branch, given the checks above. However, 01757 // if that is the case, it should still be correct to call this 01758 // function if the calling process owns no column indices. 01759 numEntries = 0; 01760 } 01761 } 01762 01763 01766 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01767 void 01768 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01769 getGlobalRowCopy (GlobalOrdinal globalRow, 01770 const ArrayView<GlobalOrdinal>& indices, 01771 size_t& NumIndices) const 01772 { 01773 // One of the following is true: 01774 // 01775 // 1. The calling process owns no column indices. 01776 // 2. The graph stores global column indices. 01777 // 3. The graph has a column Map that we can use to convert any 01778 // stored local indices into global indices for output. 01779 01780 const char tfecfFuncName[] = "getGlobalRowCopy"; 01781 const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow); 01782 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01783 localRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), 01784 std::runtime_error, ": globalRow (== " << globalRow << ") is not owned " 01785 "by (i.e., is not in the row Map on) the calling process with rank " 01786 << getComm ()->getRank () << "."); 01787 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01788 ! hasRowInfo (), std::runtime_error, 01789 ": graph row information was deleted at fillComplete()."); 01790 01791 const RowInfo rowinfo = getRowInfo (static_cast<size_t> (localRow)); 01792 NumIndices = rowinfo.numEntries; 01793 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01794 static_cast<size_t> (indices.size ()) < NumIndices, std::runtime_error, 01795 ": specified storage (size == " << indices.size () << ") does not suffice " 01796 "to hold all entries in this row (NumIndices == " << NumIndices << ")."); 01797 01798 if (isLocallyIndexed ()) { 01799 ArrayView<const LocalOrdinal> lview = getLocalView (rowinfo); 01800 for (size_t j = 0; j < NumIndices; ++j) { 01801 indices[j] = colMap_->getGlobalElement (lview[j]); 01802 } 01803 } 01804 else if (isGloballyIndexed ()) { 01805 ArrayView<const GlobalOrdinal> gview = getGlobalView (rowinfo); 01806 std::copy (gview.begin (), gview.begin () + NumIndices, indices.begin ()); 01807 } 01808 } 01809 01810 01813 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01814 void 01815 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01816 getLocalRowView (LocalOrdinal localRow, ArrayView<const LocalOrdinal> &indices) const 01817 { 01818 const char tfecfFuncName[] = "getLocalRowView"; 01819 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01820 isGloballyIndexed(), std::runtime_error, 01821 ": The graph is globally indexed, so we cannot return a view with local " 01822 "column indices. Use getLocalRowCopy() instead."); 01823 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01824 ! hasRowInfo (), std::runtime_error, 01825 ": graph row information was deleted at fillComplete()."); 01826 01827 indices = null; 01828 if (rowMap_->isNodeLocalElement (localRow)) { 01829 const RowInfo rowinfo = getRowInfo(localRow); 01830 if (rowinfo.numEntries > 0) { 01831 indices = getLocalView (rowinfo); 01832 indices = indices (0, rowinfo.numEntries); 01833 } 01834 #ifdef HAVE_TPETRA_DEBUG 01835 const size_t numEntriesInView = static_cast<size_t> (indices.size ()); 01836 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01837 numEntriesInView != getNumEntriesInLocalRow (localRow), 01838 std::logic_error, ": The number of entries " << numEntriesInView 01839 << " in the view to return must be the same as getNumEntriesInLocalRow" 01840 "(localRow=" << localRow << ") = " << getNumEntriesInLocalRow (localRow) 01841 << ", but this is not the case. " 01842 "Please report this bug to the Tpetra developers."); 01843 #endif // HAVE_TPETRA_DEBUG 01844 } 01845 } 01846 01847 01850 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01851 void 01852 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01853 getGlobalRowView (GlobalOrdinal globalRow, 01854 ArrayView<const GlobalOrdinal>& indices) const 01855 { 01856 using Teuchos::as; 01857 const char tfecfFuncName[] = "getGlobalRowView"; 01858 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01859 isLocallyIndexed (), std::runtime_error, 01860 ": The graph is locally indexed, so we cannot return a view with global " 01861 "column indices. Use getGlobalRowCopy() instead."); 01862 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01863 ! hasRowInfo (), std::runtime_error, 01864 ": graph row information was deleted at fillComplete()."); 01865 01866 // We don't have to call isNodeGlobalElement() to test whether 01867 // globalRow belongs to the calling process. That method requires 01868 // a global to local lookup anyway, and getLocalElement() returns 01869 // invalid() anyway if the global index wasn't found. 01870 const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow); 01871 indices = null; 01872 if (localRow != Teuchos::OrdinalTraits<LocalOrdinal>::invalid ()) { 01873 const RowInfo rowInfo = getRowInfo (as<size_t> (localRow)); 01874 if (rowInfo.numEntries > 0) { 01875 indices = (getGlobalView (rowInfo)) (0, rowInfo.numEntries); 01876 } 01877 } 01878 #ifdef HAVE_TPETRA_DEBUG 01879 using std::endl; 01880 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01881 as<size_t> (indices.size ()) != getNumEntriesInGlobalRow (globalRow), 01882 std::logic_error, 01883 ": Violated stated post-conditions:" 01884 << " indices.size () = " << indices.size () << endl 01885 << " as<size_t> (indices.size ()) = " << as<size_t> (indices.size ()) 01886 << endl << " getNumEntriesInGlobalRow (globalRow = " << globalRow 01887 << ") = " << getNumEntriesInGlobalRow (globalRow) << endl 01888 << "Please report this bug to the Tpetra developers."); 01889 #endif // HAVE_TPETRA_DEBUG 01890 } 01891 01892 01895 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01896 void 01897 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01898 insertLocalIndices (const LocalOrdinal localRow, 01899 const ArrayView<const LocalOrdinal> &indices) 01900 { 01901 using Teuchos::ArrayView; 01902 const char tfecfFuncName[] = "insertLocalIndices"; 01903 01904 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01905 isFillActive() == false, std::runtime_error, 01906 ": requires that fill is active."); 01907 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01908 isGloballyIndexed() == true, std::runtime_error, 01909 ": graph indices are global; use insertGlobalIndices()."); 01910 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01911 hasColMap() == false, std::runtime_error, 01912 ": cannot insert local indices without a column map."); 01913 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01914 rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, 01915 ": row does not belong to this node."); 01916 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01917 hasRowInfo() == false, std::runtime_error, 01918 ": graph row information was deleted at fillComplete()."); 01919 if (! indicesAreAllocated ()) { 01920 allocateIndices (LocalIndices); 01921 } 01922 01923 #ifdef HAVE_TPETRA_DEBUG 01924 // In a debug build, if the graph has a column Map, test whether 01925 // any of the given column indices are not in the column Map. 01926 // Keep track of the invalid column indices so we can tell the 01927 // user about them. 01928 if (hasColMap ()) { 01929 using Teuchos::Array; 01930 using Teuchos::toString; 01931 using std::endl; 01932 typedef LocalOrdinal LO; 01933 typedef typename ArrayView<const LO>::size_type size_type; 01934 01935 const map_type& colMap = * (getColMap ()); 01936 Array<LO> badColInds; 01937 bool allInColMap = true; 01938 for (size_type k = 0; k < indices.size (); ++k) { 01939 if (! colMap.isNodeLocalElement (indices[k])) { 01940 allInColMap = false; 01941 badColInds.push_back (indices[k]); 01942 } 01943 } 01944 if (! allInColMap) { 01945 std::ostringstream os; 01946 os << "Tpetra::CrsMatrix::insertLocalIndices: You attempted to insert " 01947 "entries in owned row " << localRow << ", at the following column " 01948 "indices: " << toString (indices) << "." << endl; 01949 os << "Of those, the following indices are not in the column Map on " 01950 "this process: " << toString (badColInds) << "." << endl << "Since " 01951 "the graph has a column Map already, it is invalid to insert entries " 01952 "at those locations."; 01953 TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ()); 01954 } 01955 } 01956 #endif // HAVE_TPETRA_DEBUG 01957 01958 insertLocalIndicesImpl (localRow, indices); 01959 01960 #ifdef HAVE_TPETRA_DEBUG 01961 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01962 indicesAreAllocated() == false || isLocallyIndexed() == false, 01963 std::logic_error, 01964 ": Violated stated post-conditions. Please contact Tpetra team."); 01965 #endif 01966 } 01967 01968 01971 template <class LocalOrdinal, class GlobalOrdinal, class Node> 01972 void 01973 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 01974 insertLocalIndicesFiltered (const LocalOrdinal localRow, 01975 const ArrayView<const LocalOrdinal> &indices) 01976 { 01977 typedef LocalOrdinal LO; 01978 const char tfecfFuncName[] = "insertLocalIndicesFiltered"; 01979 01980 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01981 isFillActive() == false, std::runtime_error, 01982 ": requires that fill is active."); 01983 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01984 isGloballyIndexed() == true, std::runtime_error, 01985 ": graph indices are global; use insertGlobalIndices()."); 01986 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01987 hasColMap() == false, std::runtime_error, 01988 ": cannot insert local indices without a column map."); 01989 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01990 rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, 01991 ": row does not belong to this node."); 01992 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 01993 hasRowInfo() == false, std::runtime_error, 01994 ": graph row information was deleted at fillComplete()."); 01995 if (indicesAreAllocated() == false) { 01996 allocateIndices(LocalIndices); 01997 } 01998 01999 // If we have a column map, use it to filter the entries. 02000 if (hasColMap ()) { 02001 Array<LO> filtered_indices(indices); 02002 SLocalGlobalViews inds_view; 02003 SLocalGlobalNCViews inds_ncview; 02004 inds_ncview.linds = filtered_indices(); 02005 const size_t numFilteredEntries = 02006 filterIndices<LocalIndices>(inds_ncview); 02007 inds_view.linds = filtered_indices (0, numFilteredEntries); 02008 insertLocalIndicesImpl(localRow, inds_view.linds); 02009 } 02010 else { 02011 insertLocalIndicesImpl(localRow, indices); 02012 } 02013 #ifdef HAVE_TPETRA_DEBUG 02014 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02015 indicesAreAllocated() == false || isLocallyIndexed() == false, 02016 std::logic_error, 02017 ": Violated stated post-conditions. Please contact Tpetra team."); 02018 #endif 02019 } 02020 02021 02024 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02025 void 02026 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02027 insertGlobalIndices (const GlobalOrdinal grow, 02028 const ArrayView<const GlobalOrdinal> &indices) 02029 { 02030 typedef LocalOrdinal LO; 02031 typedef GlobalOrdinal GO; 02032 typedef typename ArrayView<const GO>::size_type size_type; 02033 const char tfecfFuncName[] = "insertGlobalIndices"; 02034 02035 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02036 isLocallyIndexed() == true, std::runtime_error, 02037 ": graph indices are local; use insertLocalIndices()."); 02038 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02039 hasRowInfo() == false, std::runtime_error, 02040 ": graph row information was deleted at fillComplete()."); 02041 // This can't really be satisfied for now, because if we are 02042 // fillComplete(), then we are local. In the future, this may 02043 // change. However, the rule that modification require active 02044 // fill will not change. 02045 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02046 isFillActive() == false, std::runtime_error, 02047 ": You are not allowed to call this method if fill is not active. " 02048 "If fillComplete has been called, you must first call resumeFill " 02049 "before you may insert indices."); 02050 if (! indicesAreAllocated ()) { 02051 allocateIndices (GlobalIndices); 02052 } 02053 const LO myRow = rowMap_->getLocalElement (grow); 02054 if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) { 02055 #ifdef HAVE_TPETRA_DEBUG 02056 if (hasColMap ()) { 02057 using std::endl; 02058 const map_type& colMap = * (getColMap ()); 02059 // In a debug build, keep track of the nonowned ("bad") column 02060 // indices, so that we can display them in the exception 02061 // message. In a release build, just ditch the loop early if 02062 // we encounter a nonowned column index. 02063 Array<GO> badColInds; 02064 bool allInColMap = true; 02065 for (size_type k = 0; k < indices.size (); ++k) { 02066 if (! colMap.isNodeGlobalElement (indices[k])) { 02067 allInColMap = false; 02068 badColInds.push_back (indices[k]); 02069 } 02070 } 02071 if (! allInColMap) { 02072 std::ostringstream os; 02073 os << "Tpetra::CrsGraph::insertGlobalIndices: You attempted to insert " 02074 "entries in owned row " << grow << ", at the following column " 02075 "indices: " << toString (indices) << "." << endl; 02076 os << "Of those, the following indices are not in the column Map on " 02077 "this process: " << toString (badColInds) << "." << endl << "Since " 02078 "the matrix has a column Map already, it is invalid to insert " 02079 "entries at those locations."; 02080 TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ()); 02081 } 02082 } 02083 #endif // HAVE_TPETRA_DEBUG 02084 insertGlobalIndicesImpl (myRow, indices); 02085 } 02086 else { // a nonlocal row 02087 const size_type numIndices = indices.size (); 02088 // This creates the Array if it doesn't exist yet. std::map's 02089 // operator[] does a lookup each time, so it's better to pull 02090 // nonlocals_[grow] out of the loop. 02091 std::vector<GO>& nonlocalRow = nonlocals_[grow]; 02092 for (size_type k = 0; k < numIndices; ++k) { 02093 nonlocalRow.push_back (indices[k]); 02094 } 02095 } 02096 #ifdef HAVE_TPETRA_DEBUG 02097 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02098 indicesAreAllocated() == false || isGloballyIndexed() == false, 02099 std::logic_error, 02100 ": Violated stated post-conditions. Please contact Tpetra team."); 02101 #endif 02102 } 02103 02104 02107 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02108 void 02109 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02110 insertGlobalIndicesFiltered (const GlobalOrdinal grow, 02111 const ArrayView<const GlobalOrdinal> &indices) 02112 { 02113 typedef LocalOrdinal LO; 02114 typedef GlobalOrdinal GO; 02115 const char tfecfFuncName[] = "insertGlobalIndicesFiltered"; 02116 02117 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02118 isLocallyIndexed() == true, std::runtime_error, 02119 ": graph indices are local; use insertLocalIndices()."); 02120 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02121 hasRowInfo() == false, std::runtime_error, 02122 ": graph row information was deleted at fillComplete()."); 02123 // This can't really be satisfied for now, because if we are 02124 // fillComplete(), then we are local. In the future, this may 02125 // change. However, the rule that modification require active 02126 // fill will not change. 02127 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02128 isFillActive() == false, std::runtime_error, 02129 ": You are not allowed to call this method if fill is not active. " 02130 "If fillComplete has been called, you must first call resumeFill " 02131 "before you may insert indices."); 02132 if (indicesAreAllocated() == false) { 02133 allocateIndices (GlobalIndices); 02134 } 02135 const LO myRow = rowMap_->getLocalElement (grow); 02136 if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) { 02137 // If we have a column map, use it to filter the entries. 02138 if (hasColMap ()) { 02139 Array<GO> filtered_indices(indices); 02140 SLocalGlobalViews inds_view; 02141 SLocalGlobalNCViews inds_ncview; 02142 inds_ncview.ginds = filtered_indices(); 02143 const size_t numFilteredEntries = 02144 filterIndices<GlobalIndices> (inds_ncview); 02145 inds_view.ginds = filtered_indices (0, numFilteredEntries); 02146 insertGlobalIndicesImpl(myRow, inds_view.ginds); 02147 } 02148 else { 02149 insertGlobalIndicesImpl(myRow, indices); 02150 } 02151 } 02152 else { // a nonlocal row 02153 typedef typename ArrayView<const GO>::size_type size_type; 02154 const size_type numIndices = indices.size (); 02155 for (size_type k = 0; k < numIndices; ++k) { 02156 nonlocals_[grow].push_back (indices[k]); 02157 } 02158 } 02159 #ifdef HAVE_TPETRA_DEBUG 02160 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02161 indicesAreAllocated() == false || isGloballyIndexed() == false, 02162 std::logic_error, 02163 ": Violated stated post-conditions. Please contact Tpetra team."); 02164 #endif 02165 } 02166 02169 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02170 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::removeLocalIndices(LocalOrdinal lrow) 02171 { 02172 const char tfecfFuncName[] = "removeLocalIndices()"; 02173 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active."); 02174 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isStorageOptimized() == true, std::runtime_error, ": cannot remove indices after optimizeStorage() has been called."); 02175 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true, std::runtime_error, ": graph indices are global."); 02176 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error, ": row does not belong to this node."); 02177 if (indicesAreAllocated() == false) { 02178 allocateIndices(LocalIndices); 02179 } 02180 // 02181 clearGlobalConstants(); 02182 // 02183 if (numRowEntries_ != null) { 02184 nodeNumEntries_ -= numRowEntries_[lrow]; 02185 numRowEntries_[lrow] = 0; 02186 } 02187 #ifdef HAVE_TPETRA_DEBUG 02188 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getNumEntriesInLocalRow(lrow) != 0 || indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error, 02189 ": Violated stated post-conditions. Please contact Tpetra team."); 02190 #endif 02191 } 02192 02195 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02196 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::setAllIndices(const ArrayRCP<size_t> & rowPointers,const ArrayRCP<LocalOrdinal> & columnIndices) 02197 { 02198 const char tfecfFuncName[] = "setAllIndices()"; 02199 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( hasColMap() == false, std::runtime_error, ": requires a ColMap."); 02200 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)rowPointers.size() != getNodeNumRows()+1, std::runtime_error, ": requires rowPointers.size() == getNodeNumRows()+1"); 02201 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds1D_ != Teuchos::null || gblInds1D_ != Teuchos::null, std::runtime_error, ": cannot have 1D data structures allocated."); 02202 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds2D_ != Teuchos::null || gblInds2D_ != Teuchos::null, std::runtime_error, ": cannot have 2D data structures allocated."); 02203 02204 indicesAreAllocated_ = true; 02205 indicesAreLocal_ = true; 02206 pftype_ = StaticProfile; // if the profile wasn't static before, it sure is now. 02207 lclInds1D_ = columnIndices; 02208 rowPtrs_ = rowPointers; 02209 nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()]; 02210 numAllocForAllRows_ = 0; 02211 numAllocPerRow_ = null; 02212 checkInternalState(); 02213 } 02214 02217 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02218 void 02219 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02220 getNumEntriesPerLocalRowUpperBound (Teuchos::ArrayRCP<const size_t>& boundPerLocalRow, 02221 size_t& boundForAllLocalRows, 02222 bool& boundSameForAllLocalRows) const 02223 { 02224 // The three output arguments. We assign them to the actual 02225 // output arguments at the end, in order to implement 02226 // transactional semantics. 02227 Teuchos::ArrayRCP<const size_t> numEntriesPerRow; 02228 size_t numEntriesForAll = 0; 02229 bool allRowsSame = true; 02230 02231 const ptrdiff_t numRows = static_cast<ptrdiff_t> (this->getNodeNumRows ()); 02232 02233 if (! this->indicesAreAllocated ()) { 02234 if (! this->numAllocPerRow_.is_null ()) { 02235 numEntriesPerRow = this->numAllocPerRow_; 02236 allRowsSame = false; // conservatively; we don't check the array 02237 } 02238 else { 02239 numEntriesForAll = this->numAllocForAllRows_; 02240 allRowsSame = true; 02241 } 02242 } 02243 else if (! this->numRowEntries_.is_null ()) { 02244 numEntriesPerRow = numRowEntries_; 02245 allRowsSame = false; // conservatively; we don't check the array 02246 } 02247 else if (this->nodeNumAllocated_ == 0) { 02248 numEntriesForAll = 0; 02249 allRowsSame = true; 02250 } 02251 else { 02252 // left with the case that we have optimized storage. in this 02253 // case, we have to construct a list of row sizes. 02254 TEUCHOS_TEST_FOR_EXCEPTION( 02255 this->getProfileType () != StaticProfile, std::logic_error, 02256 "Tpetra::CrsGraph::getNumEntriesPerRowUpperBound: " 02257 "The graph is not StaticProfile, but storage appears to be optimized. " 02258 "Please report this bug to the Tpetra developers."); 02259 TEUCHOS_TEST_FOR_EXCEPTION( 02260 numRows != 0 && this->rowPtrs_.size () == 0, std::logic_error, 02261 "Tpetra::CrsGraph::getNumEntriesPerRowUpperBound: " 02262 "The graph has " << numRows << " (> 0) row" << (numRows != 1 ? "s" : "") 02263 << " on the calling process, but the rowPtrs_ array has zero entries. " 02264 "Please report this bug to the Tpetra developers."); 02265 02266 Teuchos::ArrayRCP<size_t> numEnt; 02267 if (numRows != 0) { 02268 numEnt = Teuchos::arcp<size_t> (numRows); 02269 } 02270 02271 // We have to iterate through the row offsets anyway, so we 02272 // might as well check whether all rows' bounds are the same. 02273 bool allRowsReallySame = false; 02274 for (ptrdiff_t i = 0; i < numRows; ++i) { 02275 numEnt[i] = this->rowPtrs_[i+1] - this->rowPtrs_[i]; 02276 if (i != 0 && numEnt[i] != numEnt[i-1]) { 02277 allRowsReallySame = false; 02278 } 02279 } 02280 if (allRowsReallySame) { 02281 if (numRows == 0) { 02282 numEntriesForAll = 0; 02283 } else { 02284 numEntriesForAll = numEnt[1] - numEnt[0]; 02285 } 02286 allRowsSame = true; 02287 } 02288 else { 02289 numEntriesPerRow = numEnt; // Teuchos::arcp_const_cast<const size_t> (numEnt); 02290 allRowsSame = false; // conservatively; we don't check the array 02291 } 02292 } 02293 02294 TEUCHOS_TEST_FOR_EXCEPTION( 02295 numEntriesForAll != 0 && numEntriesPerRow.size () != 0, std::logic_error, 02296 "Tpetra::CrsGraph::getNumEntriesPerLocalRowUpperBound: " 02297 "numEntriesForAll and numEntriesPerRow are not consistent. The former " 02298 "is " << numEntriesForAll << ", but the latter has nonzero size " 02299 << numEntriesPerRow.size () << ". " 02300 "Please report this bug to the Tpetra developers."); 02301 02302 boundPerLocalRow = numEntriesPerRow; 02303 boundForAllLocalRows = numEntriesForAll; 02304 boundSameForAllLocalRows = allRowsSame; 02305 } 02306 02307 // TODO: in the future, globalAssemble() should use import/export functionality 02310 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02311 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::globalAssemble() 02312 { 02313 using Teuchos::Array; 02314 using Teuchos::as; 02315 using Teuchos::Comm; 02316 using Teuchos::gatherAll; 02317 using Teuchos::ireceive; 02318 using Teuchos::isend; 02319 using Teuchos::outArg; 02320 using Teuchos::REDUCE_MAX; 02321 using Teuchos::reduceAll; 02322 using Teuchos::toString; 02323 using Teuchos::waitAll; 02324 using std::endl; 02325 using std::make_pair; 02326 using std::pair; 02327 typedef GlobalOrdinal GO; 02328 typedef typename std::map<GO, std::vector<GO> >::const_iterator NLITER; 02329 typedef typename Array<GO>::size_type size_type; 02330 02331 const char tfecfFuncName[] = "globalAssemble"; // for exception macro 02332 RCP<const Comm<int> > comm = getComm(); 02333 02334 const int numImages = comm->getSize(); 02335 const int myImageID = comm->getRank(); 02336 #ifdef HAVE_TPETRA_DEBUG 02337 Teuchos::barrier (*comm); 02338 #endif // HAVE_TPETRA_DEBUG 02339 02340 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive(), std::runtime_error, 02341 ": requires that fill is active."); 02342 // Determine if any nodes have global entries to share 02343 { 02344 size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals; 02345 reduceAll (*comm, REDUCE_MAX, MyNonlocals, outArg (MaxGlobalNonlocals)); 02346 if (MaxGlobalNonlocals == 0) { 02347 return; // no entries to share 02348 } 02349 } 02350 02351 // compute a list of NLRs from nonlocals_ and use it to compute: 02352 // IdsAndRows: a vector of (id,row) pairs 02353 // NLR2Id: a map from NLR to the Id that owns it 02354 // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i 02355 // sendIDs: a list of all images I send to 02356 // recvIDs: a list of all images I receive from (constructed later) 02357 Array<pair<int, GO> > IdsAndRows; 02358 std::map<GO, int> NLR2Id; 02359 Teuchos::SerialDenseMatrix<int, char> globalNeighbors; 02360 Array<int> sendIDs, recvIDs; 02361 { 02362 // nonlocals_ contains the entries we are holding for all non-local rows 02363 // we want a list of the rows for which we have data 02364 Array<GO> NLRs; 02365 std::set<GO> setOfRows; 02366 for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) { 02367 setOfRows.insert(iter->first); 02368 } 02369 // copy the elements in the set into an Array 02370 NLRs.resize(setOfRows.size()); 02371 std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin()); 02372 02373 // get a list of ImageIDs for the non-local rows (NLRs) 02374 Array<int> NLRIds(NLRs.size()); 02375 { 02376 LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds()); 02377 int lclerror = ( stat == IDNotPresent ? 1 : 0 ); 02378 int gblerror; 02379 reduceAll<int, int> (*comm, REDUCE_MAX, lclerror, outArg (gblerror)); 02380 if (gblerror != 0) { 02381 const int myRank = comm->getRank (); 02382 std::ostringstream os; 02383 os << "On one or more processes in the communicator, " 02384 << "there were insertions into rows of the graph that do not " 02385 << "exist in the row Map on any process in the communicator." 02386 << endl << "This process " << myRank << " is " 02387 << (lclerror == 0 ? "not " : "") << "one of those offending " 02388 << "processes." << endl; 02389 if (lclerror != 0) { 02390 // If NLRIds[k] is -1, then NLRs[k] is a row index not in 02391 // the row Map. Collect this list of invalid row indices 02392 // for display in the exception message. 02393 Array<GO> invalidNonlocalRows; 02394 for (size_type k = 0; k < NLRs.size (); ++k) { 02395 if (NLRIds[k] == -1) { 02396 invalidNonlocalRows.push_back (NLRs[k]); 02397 } 02398 } 02399 const size_type numInvalid = invalidNonlocalRows.size (); 02400 os << "On this process, " << numInvalid << " nonlocal row" 02401 << (numInvalid != 1 ? "s " : " ") << " were inserted that are " 02402 << "not in the row Map on any process." << endl; 02403 // Don't print _too_ many nonlocal rows. 02404 if (numInvalid <= 100) { 02405 os << "Offending row indices: " 02406 << toString (invalidNonlocalRows ()) << endl; 02407 } 02408 } 02409 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02410 gblerror != 0, std::runtime_error, 02411 ": nonlocal entries correspond to invalid rows." 02412 << endl << os.str ()); 02413 } 02414 } 02415 02416 // build up a list of neighbors, as well as a map between NLRs and Ids 02417 // localNeighbors[i] != 0 iff I have data to send to image i 02418 // put NLRs,Ids into an array of pairs 02419 IdsAndRows.reserve(NLRs.size()); 02420 Array<char> localNeighbors(numImages, 0); 02421 typename Array<GO>::const_iterator nlr; 02422 typename Array<int>::const_iterator id; 02423 for (nlr = NLRs.begin(), id = NLRIds.begin(); 02424 nlr != NLRs.end(); ++nlr, ++id) { 02425 NLR2Id[*nlr] = *id; 02426 localNeighbors[*id] = 1; 02427 // IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr)); 02428 IdsAndRows.push_back(make_pair(*id,*nlr)); 02429 } 02430 for (int j=0; j<numImages; ++j) { 02431 if (localNeighbors[j]) { 02432 sendIDs.push_back(j); 02433 } 02434 } 02435 // sort IdsAndRows, by Ids first, then rows 02436 std::sort(IdsAndRows.begin(),IdsAndRows.end()); 02437 // gather from other nodes to form the full graph 02438 globalNeighbors.shapeUninitialized(numImages,numImages); 02439 gatherAll (*getComm(), numImages, localNeighbors.getRawPtr(), 02440 numImages * numImages, globalNeighbors.values()); 02441 // globalNeighbors at this point contains (on all images) the 02442 // connectivity between the images. 02443 // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j 02444 } 02445 02447 // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH 02448 // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID 02450 02451 // loop over all columns to know from which images I can expect to receive something 02452 for (int j=0; j<numImages; ++j) { 02453 if (globalNeighbors(myImageID,j)) { 02454 recvIDs.push_back(j); 02455 } 02456 } 02457 const size_t numRecvs = recvIDs.size(); 02458 02459 // we know how many we're sending to already 02460 // form a contiguous list of all data to be sent 02461 // track the number of entries for each ID 02462 Array<pair<GO, GO> > IJSendBuffer; 02463 Array<size_t> sendSizes(sendIDs.size(), 0); 02464 size_t numSends = 0; 02465 for (typename Array<pair<int, GO> >::const_iterator IdAndRow = IdsAndRows.begin(); 02466 IdAndRow != IdsAndRows.end(); ++IdAndRow) { 02467 int id = IdAndRow->first; 02468 GO row = IdAndRow->second; 02469 // have we advanced to a new send? 02470 if (sendIDs[numSends] != id) { 02471 numSends++; 02472 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(sendIDs[numSends] != id, 02473 std::logic_error, ": internal logic error. Contact Tpetra team."); 02474 } 02475 // copy data for row into contiguous storage 02476 for (typename std::vector<GO>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j) 02477 { 02478 IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) ); 02479 sendSizes[numSends]++; 02480 } 02481 } 02482 if (IdsAndRows.size() > 0) { 02483 numSends++; // one last increment, to make it a count instead of an index 02484 } 02485 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, ": internal logic error. Contact Tpetra team."); 02486 02487 // don't need this data anymore 02488 nonlocals_.clear(); 02489 02491 // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS 02493 02494 // Array of pending nonblocking communication requests. It's OK 02495 // to mix nonblocking send and receive requests in the same 02496 // waitAll() call. 02497 Array<RCP<Teuchos::CommRequest<int> > > requests; 02498 02499 // perform non-blocking sends: send sizes to our recipients 02500 for (size_t s = 0; s < numSends ; ++s) { 02501 // We're using a nonowning RCP because all communication 02502 // will be local to this method and the scope of our data 02503 requests.push_back (isend<int, size_t> (*comm, 02504 rcp (&sendSizes[s], false), 02505 sendIDs[s])); 02506 } 02507 // perform non-blocking receives: receive sizes from our senders 02508 Array<size_t> recvSizes (numRecvs); 02509 for (size_t r = 0; r < numRecvs; ++r) { 02510 // We're using a nonowning RCP because all communication 02511 // will be local to this method and the scope of our data 02512 requests.push_back (ireceive (*comm, rcp (&recvSizes[r], false), recvIDs[r])); 02513 } 02514 // Wait on all the nonblocking sends and receives. 02515 if (! requests.empty()) { 02516 waitAll (*comm, requests()); 02517 } 02518 #ifdef HAVE_TPETRA_DEBUG 02519 Teuchos::barrier (*comm); 02520 #endif // HAVE_TPETRA_DEBUG 02521 02522 // This doesn't necessarily deallocate the array. 02523 requests.resize (0); 02524 02526 // NOW SEND/RECEIVE ALL ROW DATA 02528 // from the size info, build the ArrayViews into IJSendBuffer 02529 Array<ArrayView<pair<GO,GO> > > sendBuffers(numSends,null); 02530 { 02531 size_t cur = 0; 02532 for (size_t s=0; s<numSends; ++s) { 02533 sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]); 02534 cur += sendSizes[s]; 02535 } 02536 } 02537 // perform non-blocking sends 02538 for (size_t s=0; s < numSends ; ++s) 02539 { 02540 // We're using a nonowning RCP because all communication 02541 // will be local to this method and the scope of our data 02542 ArrayRCP<pair<GO,GO> > tmpSendBuf = 02543 arcp (sendBuffers[s].getRawPtr(), 0, sendBuffers[s].size(), false); 02544 requests.push_back (isend<int, pair<GO,GO> > (*comm, tmpSendBuf, sendIDs[s])); 02545 } 02546 // calculate amount of storage needed for receives 02547 // setup pointers for the receives as well 02548 size_t totalRecvSize = std::accumulate (recvSizes.begin(), recvSizes.end(), 0); 02549 Array<pair<GO,GO> > IJRecvBuffer (totalRecvSize); 02550 // from the size info, build the ArrayViews into IJRecvBuffer 02551 Array<ArrayView<pair<GO,GO> > > recvBuffers (numRecvs, null); 02552 { 02553 size_t cur = 0; 02554 for (size_t r=0; r<numRecvs; ++r) { 02555 recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]); 02556 cur += recvSizes[r]; 02557 } 02558 } 02559 // perform non-blocking recvs 02560 for (size_t r = 0; r < numRecvs; ++r) { 02561 // We're using a nonowning RCP because all communication 02562 // will be local to this method and the scope of our data 02563 ArrayRCP<pair<GO,GO> > tmpRecvBuf = 02564 arcp (recvBuffers[r].getRawPtr(), 0, recvBuffers[r].size(), false); 02565 requests.push_back (ireceive (*comm, tmpRecvBuf, recvIDs[r])); 02566 } 02567 // perform waits 02568 if (! requests.empty()) { 02569 waitAll (*comm, requests()); 02570 } 02571 #ifdef HAVE_TPETRA_DEBUG 02572 Teuchos::barrier (*comm); 02573 #endif // HAVE_TPETRA_DEBUG 02574 02576 // NOW PROCESS THE RECEIVED ROW DATA 02578 // TODO: instead of adding one entry at a time, add one row at a time. 02579 // this requires resorting; they arrived sorted by sending node, 02580 // so that entries could be non-contiguous if we received 02581 // multiple entries for a particular row from different processors. 02582 // it also requires restoring the data, which may make it not worth the trouble. 02583 for (typename Array<pair<GO,GO> >::const_iterator ij = IJRecvBuffer.begin(); 02584 ij != IJRecvBuffer.end(); ++ij) 02585 { 02586 insertGlobalIndicesFiltered (ij->first, tuple<GO> (ij->second)); 02587 } 02588 checkInternalState(); 02589 } 02590 02591 02594 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02595 void 02596 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02597 resumeFill (const RCP<ParameterList> ¶ms) 02598 { 02599 const char tfecfFuncName[] = "resumeFill"; 02600 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! hasRowInfo(), std::runtime_error, 02601 ": Sorry, you cannot resume fill of the CrsGraph, since the graph's row " 02602 "information was deleted in fillComplete()."); 02603 02604 #ifdef HAVE_TPETRA_DEBUG 02605 Teuchos::barrier( *rowMap_->getComm() ); 02606 #endif // HAVE_TPETRA_DEBUG 02607 clearGlobalConstants(); 02608 lclGraph_ = null; 02609 if (params != null) this->setParameterList (params); 02610 lowerTriangular_ = false; 02611 upperTriangular_ = false; 02612 // either still sorted/merged or initially sorted/merged 02613 indicesAreSorted_ = true; 02614 noRedundancies_ = true; 02615 fillComplete_ = false; 02616 #ifdef HAVE_TPETRA_DEBUG 02617 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02618 ! isFillActive() || isFillComplete(), std::logic_error, 02619 "::resumeFill(): At end of method, either fill is not active or fill is " 02620 "complete. This violates stated post-conditions. Please report this bug " 02621 "to the Tpetra developers."); 02622 #endif // HAVE_TPETRA_DEBUG 02623 } 02624 02625 02626 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02627 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02628 fillComplete (const Teuchos::RCP<Teuchos::ParameterList>& params) 02629 { 02630 // If the graph already has domain and range Maps, don't clobber 02631 // them. If it doesn't, use the current row Map for both the 02632 // domain and range Maps. 02633 // 02634 // NOTE (mfh 28 Sep 2014): If the graph was constructed without a 02635 // column Map, and column indices are inserted which are not in 02636 // the row Map on any process, this will cause troubles. However, 02637 // that is not a common case for most applications that we 02638 // encounter, and checking for it might require more 02639 // communication. 02640 Teuchos::RCP<const map_type> domMap = this->getDomainMap (); 02641 if (domMap.is_null ()) { 02642 domMap = this->getRowMap (); 02643 } 02644 Teuchos::RCP<const map_type> ranMap = this->getRangeMap (); 02645 if (ranMap.is_null ()) { 02646 ranMap = this->getRowMap (); 02647 } 02648 this->fillComplete (domMap, ranMap, params); 02649 } 02650 02651 02652 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02653 void 02654 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02655 fillComplete (const Teuchos::RCP<const map_type>& domainMap, 02656 const Teuchos::RCP<const map_type>& rangeMap, 02657 const Teuchos::RCP<Teuchos::ParameterList>& params) 02658 { 02659 const char tfecfFuncName[] = "fillComplete"; 02660 02661 #ifdef HAVE_TPETRA_DEBUG 02662 rowMap_->getComm ()->barrier (); 02663 #endif // HAVE_TPETRA_DEBUG 02664 02665 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(), 02666 std::runtime_error, ": Graph fill state must be active (isFillActive() " 02667 "must be true) before calling fillComplete()."); 02668 02669 const int numProcs = getComm ()->getSize (); 02670 02671 // 02672 // Read and set parameters 02673 // 02674 02675 // Does the caller want to sort remote GIDs (within those owned by 02676 // the same process) in makeColMap()? 02677 if (! params.is_null ()) { 02678 if (params->isParameter ("sort column map ghost gids")) { 02679 sortGhostsAssociatedWithEachProcessor_ = 02680 params->get<bool> ("sort column map ghost gids", 02681 sortGhostsAssociatedWithEachProcessor_); 02682 } 02683 else if (params->isParameter ("Sort column Map ghost GIDs")) { 02684 sortGhostsAssociatedWithEachProcessor_ = 02685 params->get<bool> ("Sort column Map ghost GIDs", 02686 sortGhostsAssociatedWithEachProcessor_); 02687 } 02688 } 02689 02690 // If true, the caller promises that no process did nonlocal 02691 // changes since the last call to fillComplete. 02692 bool assertNoNonlocalInserts = false; 02693 if (! params.is_null ()) { 02694 assertNoNonlocalInserts = 02695 params->get<bool> ("No Nonlocal Changes", assertNoNonlocalInserts); 02696 } 02697 02698 // 02699 // Allocate indices, if they haven't already been allocated 02700 // 02701 if (! indicesAreAllocated ()) { 02702 if (hasColMap ()) { 02703 // We have a column Map, so use local indices. 02704 allocateIndices (LocalIndices); 02705 } else { 02706 // We don't have a column Map, so use global indices. 02707 allocateIndices (GlobalIndices); 02708 } 02709 } 02710 02711 // 02712 // Do global assembly, if requested and if the communicator 02713 // contains more than one process. 02714 // 02715 const bool mayNeedGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1; 02716 if (mayNeedGlobalAssemble) { 02717 // This first checks if we need to do global assembly. 02718 // The check costs a single all-reduce. 02719 globalAssemble (); 02720 } 02721 else { 02722 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02723 numProcs > 1 && nonlocals_.size() > 0, std::runtime_error, 02724 ":" << std::endl << "The graph's communicator contains only one " 02725 "process, but there are nonlocal entries. " << std::endl << 02726 "This probably means that invalid entries were added to the graph."); 02727 } 02728 02729 // Set domain and range Map. This may clear the Import / Export 02730 // objects if the new Maps differ from any old ones. 02731 setDomainRangeMaps (domainMap, rangeMap); 02732 02733 // If the graph does not already have a column Map (either from 02734 // the user constructor calling the version of the constructor 02735 // that takes a column Map, or from a previous fillComplete call), 02736 // then create it. 02737 if (! hasColMap ()) { 02738 makeColMap (); 02739 } 02740 02741 // Make indices local, if they aren't already. 02742 // The method doesn't do any work if the indices are already local. 02743 makeIndicesLocal (); 02744 02745 if (! isSorted ()) { 02746 // If this process has no indices, then CrsGraph considers it 02747 // already trivially sorted. Thus, this method need not be 02748 // called on all processes in the row Map's communicator. 02749 sortAllIndices (); 02750 } 02751 if (! isMerged ()) { 02752 mergeAllIndices (); 02753 } 02754 makeImportExport(); // Make Import and Export objects, if necessary 02755 computeGlobalConstants (); 02756 fillLocalGraph (params); 02757 fillComplete_ = true; 02758 02759 #ifdef HAVE_TPETRA_DEBUG 02760 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02761 isFillActive() == true || isFillComplete() == false, std::logic_error, 02762 ": Violated stated post-conditions. Please contact Tpetra team."); 02763 #endif 02764 02765 checkInternalState (); 02766 } 02767 02768 02771 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02772 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::expertStaticFillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap, 02773 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & rangeMap, 02774 const RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > &importer, 02775 const RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > &exporter, 02776 const RCP<ParameterList> ¶ms) 02777 { 02778 #ifdef HAVE_TPETRA_DEBUG 02779 rowMap_->getComm ()->barrier (); 02780 #endif // HAVE_TPETRA_DEBUG 02781 const char tfecfFuncName[] = "expertStaticFillComplete()"; 02782 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillComplete() == true || hasColMap() == false, 02783 std::runtime_error, ": fillComplete cannot have already been called and a ColMap is required."); 02784 02785 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( getNodeNumRows() > 0 && rowPtrs_==Teuchos::null, 02786 std::runtime_error, ": a matrix will getNodeNumRows()>0 requires rowptr to be set."); 02787 02788 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( domainMap == Teuchos::null || rangeMap == Teuchos::null, 02789 std::runtime_error, ": requires a non-null domainMap & rangeMap."); 02790 02791 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( pftype_ !=StaticProfile, 02792 std::runtime_error, ": requires StaticProfile."); 02793 02794 // Note: We don't need to do the following things which are normally done in fillComplete: 02795 // allocateIndices, globalAssemble, makeColMap, makeIndicesLocal, sortAllIndices, mergeAllIndices 02796 02797 // Note: Need to do this so computeGlobalConstants & fillLocalGraph work 02798 nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()]; 02799 02800 // Constants from allocateIndices 02801 numAllocForAllRows_ = 0; 02802 numAllocPerRow_ = null; 02803 indicesAreAllocated_ = true; 02804 02805 // Constants from makeIndicesLocal 02806 indicesAreLocal_ = true; 02807 indicesAreGlobal_ = false; 02808 02809 // set domain/range map: may clear the import/export objects 02810 setDomainRangeMaps(domainMap,rangeMap); 02811 02812 // Presume the user sorted and merged the arrays first 02813 indicesAreSorted_ = true; 02814 noRedundancies_ = true; 02815 02816 // makeImportExport won't create a new importer/exporter if I set one here first. 02817 importer_=Teuchos::null; 02818 exporter_=Teuchos::null; 02819 if(importer != Teuchos::null) { 02820 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!importer->getSourceMap()->isSameAs(*getDomainMap()) || !importer->getTargetMap()->isSameAs(*getColMap()), 02821 std::invalid_argument,": importer does not match matrix maps."); 02822 importer_ = importer; 02823 02824 } 02825 if(exporter != Teuchos::null) { 02826 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!exporter->getSourceMap()->isSameAs(*getRowMap()) || !exporter->getTargetMap()->isSameAs(*getRangeMap()), 02827 std::invalid_argument,": exporter does not match matrix maps."); 02828 exporter_ = exporter; 02829 } 02830 makeImportExport(); 02831 02832 // Compute the constants 02833 computeGlobalConstants(); 02834 02835 // Since we have a StaticProfile, fillLocalGraph will do the right thing... 02836 fillLocalGraph(params); 02837 fillComplete_ = true; 02838 02839 #ifdef HAVE_TPETRA_DEBUG 02840 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team."); 02841 #endif 02842 checkInternalState(); 02843 } 02844 02845 02846 02849 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02850 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::fillLocalGraph(const RCP<ParameterList> ¶ms) 02851 { 02852 using Teuchos::ArrayRCP; 02853 using Teuchos::null; 02854 using Teuchos::ParameterList; 02855 using Teuchos::parameterList; 02856 using Teuchos::RCP; 02857 using Teuchos::sublist; 02858 02859 const size_t numRows = getNodeNumRows (); 02860 ArrayRCP<size_t> ptrs; 02861 ArrayRCP<LocalOrdinal> inds; 02862 bool requestOptimizedStorage = true; 02863 if (params != null && params->get ("Optimize Storage", true) == false) { 02864 requestOptimizedStorage = false; 02865 } 02866 02867 RCP<Node> node = getRowMap ()->getNode (); 02868 if (getProfileType () == DynamicProfile) { 02869 // Storage is currently in 2-D ("array of arrays") format. 02870 // Pack into 1-D packed storage. 02871 ptrs = sparse_ops_type::allocRowPtrs (node, numRowEntries_ ()); 02872 inds = sparse_ops_type::template allocStorage<LocalOrdinal> (node, ptrs ()); 02873 for (size_t row = 0; row < numRows; ++row) { 02874 const size_t numentrs = numRowEntries_[row]; 02875 std::copy (lclInds2D_[row].begin (), 02876 lclInds2D_[row].begin () + numentrs, 02877 inds + ptrs[row]); 02878 } 02879 } 02880 else if (getProfileType () == StaticProfile) { 02881 // We already have 1-D storage, but it might not yet be packed. 02882 if (nodeNumEntries_ != nodeNumAllocated_) { // need to pack storage 02883 ptrs = sparse_ops_type::allocRowPtrs (node, numRowEntries_ ()); 02884 inds = sparse_ops_type::template allocStorage<LocalOrdinal> (node, ptrs ()); 02885 for (size_t row = 0; row < numRows; ++row) { 02886 const size_t numentrs = numRowEntries_[row]; 02887 std::copy (lclInds1D_ + rowPtrs_[row], 02888 lclInds1D_ + rowPtrs_[row] + numentrs, 02889 inds + ptrs[row]); 02890 } 02891 } 02892 else { // don't need to pack storage 02893 inds = lclInds1D_; 02894 ptrs = rowPtrs_; 02895 } 02896 } 02897 // can we ditch the old allocations for the packed one? 02898 if (requestOptimizedStorage) { 02899 lclInds2D_ = null; 02900 numRowEntries_ = null; 02901 // keep the new stuff 02902 lclInds1D_ = inds; 02903 rowPtrs_ = ptrs; 02904 nodeNumAllocated_ = nodeNumEntries_; 02905 pftype_ = StaticProfile; 02906 } 02907 // build the local graph, hand over the indices 02908 RCP<ParameterList> lclparams = params.is_null () ? 02909 parameterList () : 02910 sublist (params, "Local Graph"); 02911 lclGraph_ = 02912 rcp (new local_graph_type (getRowMap ()->getNodeNumElements (), 02913 getColMap ()->getNodeNumElements (), 02914 getRowMap ()->getNode (), lclparams)); 02915 lclGraph_->setStructure (ptrs, inds); 02916 ptrs = null; 02917 inds = null; 02918 // finalize local graph 02919 const Teuchos::EDiag diag = getNodeNumDiags () < getNodeNumRows () ? 02920 Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG; 02921 Teuchos::EUplo uplo = Teuchos::UNDEF_TRI; 02922 if (isUpperTriangular ()) { 02923 uplo = Teuchos::UPPER_TRI; 02924 } 02925 else if (isLowerTriangular ()) { 02926 uplo = Teuchos::LOWER_TRI; 02927 } 02928 sparse_ops_type::finalizeGraph (uplo, diag, *lclGraph_, params); 02929 } 02930 02931 02932 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02933 void 02934 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02935 replaceColMap (const Teuchos::RCP<const map_type>& newColMap) 02936 { 02937 // NOTE: This safety check matches the code, but not the documentation of Crsgraph 02938 // 02939 // FIXME (mfh 18 Aug 2014) This will break if the calling process 02940 // has no entries, because in that case, it is neither locally nor 02941 // globally indexed. 02942 const char tfecfFuncName[] = "replaceColMap: "; 02943 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02944 isLocallyIndexed () || isGloballyIndexed (), std::runtime_error, 02945 "Requires matching maps and non-static graph."); 02946 colMap_ = newColMap; 02947 } 02948 02949 02950 template <class LocalOrdinal, class GlobalOrdinal, class Node> 02951 void 02952 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 02953 reindexColumns (const Teuchos::RCP<const map_type>& newColMap, 02954 const Teuchos::RCP<const import_type>& newImport, 02955 const bool sortIndicesInEachRow) 02956 { 02957 using Teuchos::REDUCE_MIN; 02958 using Teuchos::reduceAll; 02959 using Teuchos::RCP; 02960 typedef GlobalOrdinal GO; 02961 typedef LocalOrdinal LO; 02962 const char tfecfFuncName[] = "reindexColumns: "; 02963 02964 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 02965 isFillComplete (), std::runtime_error, "The graph is fill complete " 02966 "(isFillComplete() returns true). You must call resumeFill() before " 02967 "you may call this method."); 02968 02969 // mfh 19 Aug 2014: This method does NOT redistribute data; it 02970 // doesn't claim to do the work of an Import or Export. This 02971 // means that for all processes, the calling process MUST own all 02972 // column indices, in both the old column Map (if it exists) and 02973 // the new column Map. We check this via an all-reduce. 02974 // 02975 // Some processes may be globally indexed, others may be locally 02976 // indexed, and others (that have no graph entries) may be 02977 // neither. This method will NOT change the graph's current 02978 // state. If it's locally indexed, it will stay that way, and 02979 // vice versa. It would easy to add an option to convert indices 02980 // from global to local, so as to save a global-to-local 02981 // conversion pass. However, we don't do this here. The intended 02982 // typical use case is that the graph already has a column Map and 02983 // is locally indexed, and this is the case for which we optimize. 02984 02985 const size_t lclNumRows = getNodeNumRows (); 02986 02987 // Attempt to convert indices to the new column Map's version of 02988 // local. This will fail if on the calling process, the graph has 02989 // indices that are not on that process in the new column Map. 02990 // After the local conversion attempt, we will do an all-reduce to 02991 // see if any processes failed. 02992 02993 // If this is false, then either the graph contains a column index 02994 // which is invalid in the CURRENT column Map, or the graph is 02995 // locally indexed but currently has no column Map. In either 02996 // case, there is no way to convert the current local indices into 02997 // global indices, so that we can convert them into the new column 02998 // Map's local indices. It's possible for this to be true on some 02999 // processes but not others, due to replaceColMap. 03000 bool allCurColIndsValid = true; 03001 // On the calling process, are all valid current column indices 03002 // also in the new column Map on the calling process? In other 03003 // words, does local reindexing suffice, or should the user have 03004 // done an Import or Export instead? 03005 bool localSuffices = true; 03006 03007 // Final arrays for the local indices. We will allocate exactly 03008 // one of these ONLY if the graph is locally indexed on the 03009 // calling process, and ONLY if the graph has one or more entries 03010 // (is not empty) on the calling process. In that case, we 03011 // allocate the first (1-D storage) if the graph has a static 03012 // profile, else we allocate the second (2-D storage). 03013 Teuchos::ArrayRCP<LO> newLclInds1D; 03014 Teuchos::ArrayRCP<Teuchos::Array<LO> > newLclInds2D; 03015 03016 // If indices aren't allocated, that means the calling process 03017 // owns no entries in the graph. Thus, there is nothing to 03018 // convert, and it trivially succeeds locally. 03019 if (indicesAreAllocated ()) { 03020 if (isLocallyIndexed ()) { 03021 if (hasColMap ()) { // locally indexed, and currently has a column Map 03022 const map_type& oldColMap = * (getColMap ()); 03023 if (pftype_ == StaticProfile) { 03024 // Allocate storage for the new local indices. 03025 RCP<Node> node = getRowMap ()->getNode (); 03026 newLclInds1D = 03027 sparse_ops_type::template allocStorage<LO> (node, rowPtrs_ ()); 03028 03029 // Attempt to convert the new indices locally. 03030 for (size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) { 03031 const RowInfo rowInfo = getRowInfo (lclRow); 03032 const size_t beg = rowInfo.offset1D; 03033 const size_t end = beg + rowInfo.numEntries; 03034 for (size_t k = beg; k < end; ++k) { 03035 const LO oldLclCol = lclInds1D_[k]; 03036 if (oldLclCol == Teuchos::OrdinalTraits<LO>::invalid ()) { 03037 allCurColIndsValid = false; 03038 break; // Stop at the first invalid index 03039 } 03040 const GO gblCol = oldColMap.getGlobalElement (oldLclCol); 03041 03042 // The above conversion MUST succeed. Otherwise, the 03043 // current local index is invalid, which means that 03044 // the graph was constructed incorrectly. 03045 if (gblCol == Teuchos::OrdinalTraits<GO>::invalid ()) { 03046 allCurColIndsValid = false; 03047 break; // Stop at the first invalid index 03048 } 03049 else { 03050 const LO newLclCol = newColMap->getLocalElement (gblCol); 03051 if (newLclCol == Teuchos::OrdinalTraits<LO>::invalid ()) { 03052 localSuffices = false; 03053 break; // Stop at the first invalid index 03054 } 03055 newLclInds1D[k] = newLclCol; 03056 } 03057 } // for each entry in the current row 03058 } // for each locally owned row 03059 } 03060 else { // pftype_ == DynamicProfile 03061 // Allocate storage for the new local indices. We only 03062 // allocate the outer array here; we will allocate the 03063 // inner arrays below. 03064 newLclInds2D = Teuchos::arcp<Teuchos::Array<LO> > (lclNumRows); 03065 03066 // Attempt to convert the new indices locally. 03067 for (size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) { 03068 const RowInfo rowInfo = getRowInfo (lclRow); 03069 newLclInds2D.resize (rowInfo.allocSize); 03070 03071 Teuchos::ArrayView<const LO> oldLclRowView = getLocalView (rowInfo); 03072 Teuchos::ArrayView<LO> newLclRowView = (newLclInds2D[lclRow]) (); 03073 03074 for (size_t k = 0; k < rowInfo.numEntries; ++k) { 03075 const LO oldLclCol = oldLclRowView[k]; 03076 if (oldLclCol == Teuchos::OrdinalTraits<LO>::invalid ()) { 03077 allCurColIndsValid = false; 03078 break; // Stop at the first invalid index 03079 } 03080 const GO gblCol = oldColMap.getGlobalElement (oldLclCol); 03081 03082 // The above conversion MUST succeed. Otherwise, the 03083 // local index is invalid and the graph is wrong. 03084 if (gblCol == Teuchos::OrdinalTraits<GO>::invalid ()) { 03085 allCurColIndsValid = false; 03086 break; // Stop at the first invalid index 03087 } 03088 else { 03089 const LO newLclCol = newColMap->getLocalElement (gblCol); 03090 if (newLclCol == Teuchos::OrdinalTraits<LO>::invalid ()) { 03091 localSuffices = false; 03092 break; // Stop at the first invalid index. 03093 } 03094 newLclRowView[k] = newLclCol; 03095 } 03096 } // for each entry in the current row 03097 } // for each locally owned row 03098 } // pftype_ 03099 } 03100 else { // locally indexed, but no column Map 03101 // This case is only possible if replaceColMap() was called 03102 // with a null argument on the calling process. It's 03103 // possible, but it means that this method can't possibly 03104 // succeed, since we have no way of knowing how to convert 03105 // the current local indices to global indices. 03106 allCurColIndsValid = false; 03107 } 03108 } 03109 else { // globally indexed 03110 // If the graph is globally indexed, we don't need to save 03111 // local indices, but we _do_ need to know whether the current 03112 // global indices are valid in the new column Map. We may 03113 // need to do a getRemoteIndexList call to find this out. 03114 // 03115 // In this case, it doesn't matter whether the graph currently 03116 // has a column Map. We don't need the old column Map to 03117 // convert from global indices to the _new_ column Map's local 03118 // indices. Furthermore, we can use the same code, whether 03119 // the graph is static or dynamic profile. 03120 03121 // Test whether the current global indices are in the new 03122 // column Map on the calling process. 03123 for (size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) { 03124 const RowInfo rowInfo = getRowInfo (lclRow); 03125 Teuchos::ArrayView<const GO> oldGblRowView = getGlobalView (rowInfo); 03126 for (size_t k = 0; k < rowInfo.numEntries; ++k) { 03127 const GO gblCol = oldGblRowView[k]; 03128 if (! newColMap->isNodeGlobalElement (gblCol)) { 03129 localSuffices = false; 03130 break; // Stop at the first invalid index 03131 } 03132 } // for each entry in the current row 03133 } // for each locally owned row 03134 } // locally or globally indexed 03135 } // whether indices are allocated 03136 03137 // Do an all-reduce to check both possible error conditions. 03138 int lclSuccess[2]; 03139 lclSuccess[0] = allCurColIndsValid ? 1 : 0; 03140 lclSuccess[1] = localSuffices ? 1 : 0; 03141 int gblSuccess[2]; 03142 gblSuccess[0] = 0; 03143 gblSuccess[1] = 0; 03144 RCP<const Teuchos::Comm<int> > comm = 03145 getRowMap ().is_null () ? Teuchos::null : getRowMap ()->getComm (); 03146 if (! comm.is_null ()) { 03147 reduceAll<int, int> (*comm, REDUCE_MIN, 2, lclSuccess, gblSuccess); 03148 } 03149 03150 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03151 gblSuccess[0] == 0, std::runtime_error, "It is not possible to continue." 03152 " The most likely reason is that the graph is locally indexed, but the " 03153 "column Map is missing (null) on some processes, due to a previous call " 03154 "to replaceColMap()."); 03155 03156 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03157 gblSuccess[1] == 0, std::runtime_error, "On some process, the graph " 03158 "contains column indices that are in the old column Map, but not in the " 03159 "new column Map (on that process). This method does NOT redistribute " 03160 "data; it does not claim to do the work of an Import or Export operation." 03161 " This means that for all processess, the calling process MUST own all " 03162 "column indices, in both the old column Map and the new column Map. In " 03163 "this case, you will need to do an Import or Export operation to " 03164 "redistribute data."); 03165 03166 // Commit the results. 03167 if (isLocallyIndexed ()) { 03168 if (pftype_ == StaticProfile) { 03169 lclInds1D_ = newLclInds1D; 03170 } else { // dynamic profile 03171 lclInds2D_ = newLclInds2D; 03172 } 03173 // We've reindexed, so we don't know if the indices are sorted. 03174 // 03175 // FIXME (mfh 17 Sep 2014) It could make sense to check this, 03176 // since we're already going through all the indices above. We 03177 // could also sort each row in place; that way, we would only 03178 // have to make one pass over the rows. 03179 indicesAreSorted_ = false; 03180 if (sortIndicesInEachRow) { 03181 // NOTE (mfh 17 Sep 2014) The graph must be locally indexed in 03182 // order to call this method. 03183 // 03184 // FIXME (mfh 17 Sep 2014) This violates the strong exception 03185 // guarantee. It would be better to sort the new index arrays 03186 // before committing them. 03187 sortAllIndices (); 03188 } 03189 } 03190 colMap_ = newColMap; 03191 03192 if (newImport.is_null ()) { 03193 // FIXME (mfh 19 Aug 2014) Should use the above all-reduce to 03194 // check whether the input Import is null on any process. 03195 // 03196 // If the domain Map hasn't been set yet, we can't compute a new 03197 // Import object. Leave it what it is; it should be null, but 03198 // it doesn't matter. If the domain Map _has_ been set, then 03199 // compute a new Import object if necessary. 03200 if (! domainMap_.is_null ()) { 03201 if (! domainMap_->isSameAs (* newColMap)) { 03202 importer_ = Teuchos::rcp (new import_type (domainMap_, newColMap)); 03203 } else { 03204 importer_ = Teuchos::null; // don't need an Import 03205 } 03206 } 03207 } else { 03208 // The caller gave us an Import object. Assume that it's valid. 03209 importer_ = newImport; 03210 } 03211 } 03212 03213 03216 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03217 void 03218 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 03219 replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap, 03220 const Teuchos::RCP<const import_type>& newImporter) 03221 { 03222 const char prefix[] = "Tpetra::CrsGraph::replaceDomainMapAndImporter: "; 03223 TEUCHOS_TEST_FOR_EXCEPTION( 03224 colMap_.is_null (), std::invalid_argument, prefix << "You may not call " 03225 "this method unless the graph already has a column Map."); 03226 TEUCHOS_TEST_FOR_EXCEPTION( 03227 newDomainMap.is_null (), std::invalid_argument, 03228 prefix << "The new domain Map must be nonnull."); 03229 03230 #ifdef HAVE_TPETRA_DEBUG 03231 if (newImporter.is_null ()) { 03232 // It's not a good idea to put expensive operations in a macro 03233 // clause, even if they are side effect - free, because macros 03234 // don't promise that they won't evaluate their arguments more 03235 // than once. It's polite for them to do so, but not required. 03236 const bool colSameAsDom = colMap_->isSameAs (*newDomainMap); 03237 TEUCHOS_TEST_FOR_EXCEPTION( 03238 colSameAsDom, std::invalid_argument, "If the new Import is null, " 03239 "then the new domain Map must be the same as the current column Map."); 03240 } 03241 else { 03242 const bool colSameAsTgt = 03243 colMap_->isSameAs (* (newImporter->getTargetMap ())); 03244 const bool newDomSameAsSrc = 03245 newDomainMap->isSameAs (* (newImporter->getSourceMap ())); 03246 TEUCHOS_TEST_FOR_EXCEPTION( 03247 ! colSameAsTgt || ! newDomSameAsSrc, std::invalid_argument, "If the " 03248 "new Import is nonnull, then the current column Map must be the same " 03249 "as the new Import's target Map, and the new domain Map must be the " 03250 "same as the new Import's source Map."); 03251 } 03252 #endif // HAVE_TPETRA_DEBUG 03253 03254 domainMap_ = newDomainMap; 03255 importer_ = Teuchos::rcp_const_cast<import_type> (newImporter); 03256 } 03257 03258 03261 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03262 const Teuchos::RCP< 03263 const typename CrsGraph< 03264 LocalOrdinal, GlobalOrdinal, Node>::local_graph_type> 03265 CrsGraph<LocalOrdinal, GlobalOrdinal, Node>:: 03266 getLocalGraph () const 03267 { 03268 return lclGraph_; 03269 } 03270 03271 03274 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03275 const Teuchos::RCP< 03276 typename CrsGraph< 03277 LocalOrdinal, GlobalOrdinal, Node>::local_graph_type> 03278 CrsGraph<LocalOrdinal, GlobalOrdinal, Node>:: 03279 getLocalGraphNonConst () 03280 { 03281 return lclGraph_; 03282 } 03283 03284 03287 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03288 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::computeGlobalConstants() 03289 { 03290 using Teuchos::as; 03291 using Teuchos::outArg; 03292 using Teuchos::reduceAll; 03293 typedef LocalOrdinal LO; 03294 typedef GlobalOrdinal GO; 03295 typedef global_size_t GST; 03296 03297 #ifdef HAVE_TPETRA_DEBUG 03298 TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap(), std::logic_error, "Tpetra::" 03299 "CrsGraph::computeGlobalConstants: At this point, the graph should have " 03300 "a column Map, but it does not. Please report this bug to the Tpetra " 03301 "developers."); 03302 #endif // HAVE_TPETRA_DEBUG 03303 03304 // If necessary, (re)compute the local constants: nodeNumDiags_, 03305 // lowerTriangular_, upperTriangular_, and nodeMaxNumRowEntries_. 03306 if (! haveLocalConstants_) { 03307 // We have actually already computed nodeNumEntries_. 03308 // nodeNumEntries_ gets updated by methods that insert or remove 03309 // indices (including setAllIndices and 03310 // expertStaticFillComplete). Before fillComplete, its count 03311 // may include duplicate column indices in the same row. 03312 // However, mergeRowIndices and mergeRowIndicesAndValues both 03313 // subtract off merged indices in each row from the total count. 03314 // Thus, nodeNumEntries_ _should_ be accurate at this point, 03315 // meaning that we don't have to re-count it here. 03316 03317 // Reset local properties 03318 upperTriangular_ = true; 03319 lowerTriangular_ = true; 03320 nodeMaxNumRowEntries_ = 0; 03321 nodeNumDiags_ = 0; 03322 03323 // At this point, we know that we have both a row Map and a column Map. 03324 const Map<LO,GO,Node>& rowMap = *rowMap_; 03325 const Map<LO,GO,Node>& colMap = *colMap_; 03326 03327 // Go through all the entries of the graph. Count the number of 03328 // diagonal elements we encounter, and figure out whether the 03329 // graph is lower or upper triangular. Diagonal elements are 03330 // determined using global indices, with respect to the whole 03331 // graph. However, lower or upper triangularity is a local 03332 // property, and is determined using local indices. 03333 // 03334 // At this point, indices have already been sorted in each row. 03335 // That makes finding out whether the graph is lower / upper 03336 // triangular easier. 03337 if (indicesAreAllocated () && nodeNumAllocated_ > 0) { 03338 const size_t numLocalRows = getNodeNumRows (); 03339 for (size_t localRow = 0; localRow < numLocalRows; ++localRow) { 03340 const GO globalRow = rowMap.getGlobalElement (localRow); 03341 // Find the local (column) index for the diagonal element. 03342 const LO rlcid = colMap.getLocalElement (globalRow); 03343 RowInfo rowInfo = getRowInfo (localRow); 03344 ArrayView<const LO> rview = getLocalView (rowInfo); 03345 typename ArrayView<const LO>::iterator beg, end, cur; 03346 beg = rview.begin(); 03347 end = beg + rowInfo.numEntries; 03348 if (beg != end) { 03349 for (cur = beg; cur != end; ++cur) { 03350 // is this the diagonal? 03351 if (rlcid == *cur) ++nodeNumDiags_; 03352 } 03353 // Local column indices are sorted in each row. That means 03354 // the smallest column index in this row (on this process) 03355 // is *beg, and the largest column index in this row (on 03356 // this process) is *(end - 1). We know that end - 1 is 03357 // valid because beg != end. 03358 const size_t smallestCol = as<size_t> (*beg); 03359 const size_t largestCol = as<size_t> (*(end - 1)); 03360 03361 if (smallestCol < localRow) { 03362 upperTriangular_ = false; 03363 } 03364 if (localRow < largestCol) { 03365 lowerTriangular_ = false; 03366 } 03367 } 03368 // Update the max number of entries over all rows. 03369 nodeMaxNumRowEntries_ = std::max (nodeMaxNumRowEntries_, rowInfo.numEntries); 03370 } 03371 } 03372 haveLocalConstants_ = true; 03373 } // if my process doesn't have local constants 03374 03375 // Compute global constants from local constants. Processes that 03376 // already have local constants still participate in the 03377 // all-reduces, using their previously computed values. 03378 if (haveGlobalConstants_ == false) { 03379 // Promote all the nodeNum* and nodeMaxNum* quantities from 03380 // size_t to global_size_t, when doing the all-reduces for 03381 // globalNum* / globalMaxNum* results. 03382 // 03383 // FIXME (mfh 07 May 2013) Unfortunately, we either have to do 03384 // this in two all-reduces (one for the sum and the other for 03385 // the max), or use a custom MPI_Op that combines the sum and 03386 // the max. The latter might even be slower than two 03387 // all-reduces on modern network hardware. It would also be a 03388 // good idea to use nonblocking all-reduces (MPI 3), so that we 03389 // don't have to wait around for the first one to finish before 03390 // starting the second one. 03391 GST lcl[2], gbl[2]; 03392 lcl[0] = as<GST> (nodeNumEntries_); 03393 lcl[1] = as<GST> (nodeNumDiags_); 03394 reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_SUM, 03395 2, lcl, gbl); 03396 globalNumEntries_ = gbl[0]; 03397 globalNumDiags_ = gbl[1]; 03398 reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_MAX, 03399 as<GST> (nodeMaxNumRowEntries_), 03400 outArg (globalMaxNumRowEntries_)); 03401 haveGlobalConstants_ = true; 03402 } 03403 } 03404 03405 03408 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03409 void 03410 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeIndicesLocal () 03411 { 03412 typedef LocalOrdinal LO; 03413 typedef GlobalOrdinal GO; 03414 #ifdef HAVE_TPETRA_DEBUG 03415 const char tfecfFuncName[] = "makeIndicesLocal"; 03416 #endif // HAVE_TPETRA_DEBUG 03417 03418 // Transform indices to local index space 03419 const size_t nlrs = getNodeNumRows (); 03420 if (isGloballyIndexed () && nlrs > 0) { 03421 // allocate data for local indices 03422 if (getProfileType() == StaticProfile) { 03423 // If GO and LO are the same size, we can reuse the existing 03424 // array of 1-D index storage to convert column indices from 03425 // GO to LO. Otherwise, we'll just allocate a new buffer. 03426 if (nodeNumAllocated_ && sizeof (LO) == sizeof (GO)) { 03427 lclInds1D_ = arcp_reinterpret_cast<LO> (gblInds1D_).persistingView (0, nodeNumAllocated_); 03428 } 03429 else { 03430 lclInds1D_ = sparse_ops_type::template allocStorage<LO> (getRowMap ()->getNode (), rowPtrs_ ()); 03431 } 03432 for (size_t r = 0; r < nlrs; ++r) { 03433 const size_t offset = rowPtrs_[r]; 03434 const size_t numentry = numRowEntries_[r]; 03435 for (size_t j = 0; j < numentry; ++j) { 03436 const GO gid = gblInds1D_[offset + j]; 03437 const LO lid = colMap_->getLocalElement (gid); 03438 lclInds1D_[offset + j] = lid; 03439 #ifdef HAVE_TPETRA_DEBUG 03440 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03441 lclInds1D_[offset + j] == Teuchos::OrdinalTraits<LO>::invalid(), 03442 std::logic_error, 03443 ": In local row r=" << r << ", global column " << gid << " is " 03444 "not in the column Map. This should never happen. Please " 03445 "report this bug to the Tpetra developers."); 03446 #endif // HAVE_TPETRA_DEBUG 03447 } 03448 } 03449 // We've converted column indices from global to local, so we 03450 // can deallocate the global column indices (which we know are 03451 // in 1-D storage, because the graph has static profile). 03452 gblInds1D_ = null; 03453 } 03454 else { // the graph has dynamic profile (2-D index storage) 03455 lclInds2D_ = arcp<Array<LO> > (nlrs); 03456 for (size_t r = 0; r < getNodeNumRows (); ++r) { 03457 if (! gblInds2D_[r].empty ()) { 03458 const GO* const ginds = gblInds2D_[r].getRawPtr (); 03459 const size_t rna = gblInds2D_[r].size (); 03460 const size_t numentry = numRowEntries_[r]; 03461 lclInds2D_[r].resize (rna); 03462 LO* const linds = lclInds2D_[r].getRawPtr (); 03463 for (size_t j = 0; j < numentry; ++j) { 03464 GO gid = ginds[j]; 03465 LO lid = colMap_->getLocalElement (gid); 03466 linds[j] = lid; 03467 #ifdef HAVE_TPETRA_DEBUG 03468 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03469 linds[j] == Teuchos::OrdinalTraits<LO>::invalid(), 03470 std::logic_error, 03471 ": Global column ginds[j=" << j << "]=" << ginds[j] 03472 << " of local row r=" << r << " is not in the column Map. " 03473 "This should never happen. Please report this bug to the " 03474 "Tpetra developers."); 03475 #endif // HAVE_TPETRA_DEBUG 03476 } 03477 } 03478 } 03479 gblInds2D_ = null; 03480 } 03481 } // globallyIndexed() && nlrs > 0 03482 indicesAreLocal_ = true; 03483 indicesAreGlobal_ = false; 03484 checkInternalState(); 03485 } 03486 03487 03490 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03491 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::sortAllIndices() 03492 { 03493 TEUCHOS_TEST_FOR_EXCEPT(isGloballyIndexed()==true); // this should be called only after makeIndicesLocal() 03494 if (isSorted() == false) { 03495 for (size_t row=0; row < getNodeNumRows(); ++row) { 03496 RowInfo rowInfo = getRowInfo(row); 03497 sortRowIndices(rowInfo); 03498 } 03499 } 03500 // we just sorted every row 03501 indicesAreSorted_ = true; 03502 } 03503 03504 03507 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03508 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeColMap () 03509 { 03510 using Teuchos::Array; 03511 using Teuchos::ArrayView; 03512 using Teuchos::rcp; 03513 using Teuchos::REDUCE_MAX; 03514 using Teuchos::reduceAll; 03515 using std::endl; 03516 typedef LocalOrdinal LO; 03517 typedef GlobalOrdinal GO; 03518 const char tfecfFuncName[] = "makeColMap"; 03519 03520 if (hasColMap ()) { // The graph already has a column Map. 03521 // FIXME (mfh 26 Feb 2013): This currently prevents structure 03522 // changes that affect the column Map. 03523 return; 03524 } 03525 03526 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03527 isLocallyIndexed (), std::runtime_error, 03528 ": The graph is locally indexed. Calling makeColMap() to make the " 03529 "column Map requires that the graph be globally indexed."); 03530 03531 // After the calling process is done going through all of the rows 03532 // it owns, myColumns will contain the list of indices owned by 03533 // this process in the column Map. 03534 Array<GO> myColumns; 03535 03536 // If we reach this point, we don't have a column Map yet, so the 03537 // graph can't be locally indexed. Thus, isGloballyIndexed() == 03538 // false means that the graph is empty on this process, so 03539 // myColumns will be left empty. 03540 if (isGloballyIndexed ()) { 03541 // Go through all the rows, finding the populated column indices. 03542 // 03543 // Our final list of indices for the column Map constructor will 03544 // have the following properties (all of which are with respect 03545 // to the calling process): 03546 // 03547 // 1. Indices in the domain Map go first. 03548 // 2. Indices not in the domain Map follow, ordered first 03549 // contiguously by their owning process rank (in the domain 03550 // Map), then in increasing order within that. 03551 // 3. No duplicate indices. 03552 // 03553 // This imitates the ordering used by Aztec(OO) and Epetra. 03554 // Storing indices owned by the same process (in the domain Map) 03555 // contiguously permits the use of contiguous send and receive 03556 // buffers. 03557 // 03558 // We begin by partitioning the column indices into "local" GIDs 03559 // (owned by the domain Map) and "remote" GIDs (not owned by the 03560 // domain Map). We use the same order for local GIDs as the 03561 // domain Map, so we track them in place in their array. We use 03562 // an std::set (RemoteGIDSet) to keep track of remote GIDs, so 03563 // that we don't have to merge duplicates later. 03564 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid (); 03565 size_t numLocalColGIDs = 0, numRemoteColGIDs = 0; 03566 03567 // GIDisLocal[lid] == 0 if and only if local index lid in the 03568 // domain Map is remote (not local). 03569 Array<char> GIDisLocal (domainMap_->getNodeNumElements (), 0); 03570 std::set<GO> RemoteGIDSet; 03571 // This preserves the not-sorted Epetra order of GIDs. 03572 std::vector<GO> RemoteGIDUnorderedVector; 03573 const size_t myNumRows = getNodeNumRows (); 03574 for (size_t r = 0; r < myNumRows; ++r) { 03575 RowInfo rowinfo = getRowInfo (r); 03576 if (rowinfo.numEntries > 0) { 03577 // NOTE (mfh 02 Sep 2014) getGlobalView() returns a view of 03578 // all the space in the row, not just the occupied entries. 03579 // (This matters for the case of unpacked 1-D storage. We 03580 // might not have packed it yet.) That's why we need to 03581 // take a subview. 03582 ArrayView<const GO> rowGids = getGlobalView (rowinfo); 03583 rowGids = rowGids (0, rowinfo.numEntries); 03584 03585 for (size_t k = 0; k < rowinfo.numEntries; ++k) { 03586 const GO gid = rowGids[k]; 03587 const LO lid = domainMap_->getLocalElement (gid); 03588 if (lid != LINV) { // gid is locally owned 03589 const char alreadyFound = GIDisLocal[lid]; 03590 if (alreadyFound == 0) { 03591 GIDisLocal[lid] = static_cast<char> (1); 03592 ++numLocalColGIDs; 03593 } 03594 } 03595 else { 03596 const bool notAlreadyFound = RemoteGIDSet.insert (gid).second; 03597 if (notAlreadyFound) { // gid did not exist in the set before 03598 if (! sortGhostsAssociatedWithEachProcessor_) { 03599 // The user doesn't want to sort remote GIDs (for 03600 // each remote process); they want us to keep remote 03601 // GIDs in their original order. We do this by 03602 // stuffing each remote GID into an array as we 03603 // encounter it for the first time. The std::set 03604 // helpfully tracks first encounters. 03605 RemoteGIDUnorderedVector.push_back (gid); 03606 } 03607 ++numRemoteColGIDs; 03608 } 03609 } 03610 } // for each entry k in row r 03611 } // if row r contains > 0 entries 03612 } // for each locally owned row r 03613 03614 // Possible short-circuit for serial scenario: 03615 // 03616 // If all domain GIDs are present as column indices, then set 03617 // ColMap=DomainMap. By construction, LocalGIDs is a subset of 03618 // DomainGIDs. 03619 // 03620 // If we have 03621 // * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs, 03622 // and 03623 // * Number of local GIDs is number of domain GIDs 03624 // then 03625 // * LocalGIDs \subset DomainGIDs && size(LocalGIDs) == 03626 // size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs 03627 // on the calling process. 03628 // 03629 // We will concern ourselves only with the special case of a 03630 // serial DomainMap, obviating the need for communication. 03631 // 03632 // If 03633 // * DomainMap has a serial communicator 03634 // then we can set the column Map as the domain Map 03635 // return. Benefit: this graph won't need an Import object 03636 // later. 03637 // 03638 // Note, for a serial domain map, there can be no RemoteGIDs, 03639 // because there are no remote processes. Likely explanations 03640 // for this are: 03641 // * user submitted erroneous column indices 03642 // * user submitted erroneous domain Map 03643 if (domainMap_->getComm ()->getSize () == 1) { 03644 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03645 numRemoteColGIDs != 0, std::runtime_error, 03646 ": " << numRemoteColGIDs << " column " 03647 << (numRemoteColGIDs != 1 ? "indices are" : "index is") 03648 << " not in the domain Map." << endl 03649 << "Either these indices are invalid or the domain Map is invalid." 03650 << endl << "Remember that nonsquare matrices, or matrices where the " 03651 "row and range Maps are different, require calling the version of " 03652 "fillComplete that takes the domain and range Maps as input."); 03653 if (numLocalColGIDs == domainMap_->getNodeNumElements()) { 03654 colMap_ = domainMap_; 03655 checkInternalState (); 03656 return; 03657 } 03658 } 03659 03660 // Populate myColumns with a list of all column GIDs. Put 03661 // locally owned (in the domain Map) GIDs at the front: they 03662 // correspond to "same" and "permuted" entries between the 03663 // column Map and the domain Map. Put remote GIDs at the back. 03664 myColumns.resize (numLocalColGIDs + numRemoteColGIDs); 03665 // get pointers into myColumns for each part 03666 ArrayView<GO> LocalColGIDs = myColumns (0, numLocalColGIDs); 03667 ArrayView<GO> RemoteColGIDs = myColumns (numLocalColGIDs, numRemoteColGIDs); 03668 03669 // Copy the remote GIDs into myColumns 03670 if (sortGhostsAssociatedWithEachProcessor_) { 03671 // The std::set puts GIDs in increasing order. 03672 std::copy (RemoteGIDSet.begin(), RemoteGIDSet.end(), 03673 RemoteColGIDs.begin()); 03674 } else { 03675 // Respect the originally encountered order. 03676 std::copy (RemoteGIDUnorderedVector.begin(), 03677 RemoteGIDUnorderedVector.end(), RemoteColGIDs.begin()); 03678 } 03679 03680 // Make a list of process ranks corresponding to the remote GIDs. 03681 Array<int> RemoteImageIDs (numRemoteColGIDs); 03682 // Look up the remote process' ranks in the domain Map. 03683 { 03684 const LookupStatus stat = 03685 domainMap_->getRemoteIndexList (RemoteColGIDs, RemoteImageIDs ()); 03686 #ifdef HAVE_TPETRA_DEBUG 03687 // If any process returns IDNotPresent, then at least one of 03688 // the remote indices was not present in the domain Map. This 03689 // means that the Import object cannot be constructed, because 03690 // of incongruity between the column Map and domain Map. 03691 // This has two likely causes: 03692 // - The user has made a mistake in the column indices 03693 // - The user has made a mistake with respect to the domain Map 03694 const int missingID_lcl = (stat == IDNotPresent ? 1 : 0); 03695 int missingID_gbl = 0; 03696 reduceAll<int, int> (*getComm (), REDUCE_MAX, missingID_lcl, 03697 outArg (missingID_gbl)); 03698 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03699 missingID_gbl == 1, std::runtime_error, 03700 ": Some column indices are not in the domain Map." << endl 03701 << "Either these column indices are invalid or the domain Map is " 03702 "invalid." << endl << "Likely cause: For a nonsquare matrix, you " 03703 "must give the domain and range Maps as input to fillComplete."); 03704 #else 03705 (void) stat; // forestall compiler warning for unused variable 03706 #endif // HAVE_TPETRA_DEBUG 03707 } 03708 // Sort incoming remote column indices by their owning process 03709 // rank, so that all columns coming from a given remote process 03710 // are contiguous. This means the Import's Distributor doesn't 03711 // need to reorder data. 03712 // 03713 // NOTE (mfh 02 Sep 2014) This needs to be a stable sort, so 03714 // that it respects either of the possible orderings of GIDs 03715 // (sorted, or original order) specified above. 03716 sort2 (RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin()); 03717 03718 // Copy the local GIDs into myColumns. Two cases: 03719 // 1. If the number of Local column GIDs is the same as the 03720 // number of Local domain GIDs, we can simply read the domain 03721 // GIDs into the front part of ColIndices (see logic above 03722 // from the serial short circuit case) 03723 // 2. We step through the GIDs of the DomainMap, checking to see 03724 // if each domain GID is a column GID. We want to do this to 03725 // maintain a consistent ordering of GIDs between the columns 03726 // and the domain. 03727 03728 const size_t numDomainElts = domainMap_->getNodeNumElements (); 03729 if (numLocalColGIDs == numDomainElts) { 03730 // If the number of locally owned GIDs are the same as the 03731 // number of local domain Map elements, then the local domain 03732 // Map elements are the same as the locally owned GIDs. 03733 if (domainMap_->isContiguous ()) { 03734 // NOTE (mfh 03 Mar 2013, 02 Sep 2014) In the common case 03735 // that the domain Map is contiguous, it's more efficient to 03736 // avoid calling getNodeElementList(), since that 03737 // permanently constructs and caches the GID list in the 03738 // contiguous Map. 03739 GO curColMapGid = domainMap_->getMinGlobalIndex (); 03740 for (size_t k = 0; k < numLocalColGIDs; ++k, ++curColMapGid) { 03741 LocalColGIDs[k] = curColMapGid; 03742 } 03743 } 03744 else { 03745 ArrayView<const GO> domainElts = domainMap_->getNodeElementList (); 03746 std::copy (domainElts.begin(), domainElts.end(), LocalColGIDs.begin()); 03747 } 03748 } 03749 else { 03750 // Count the number of locally owned GIDs, both to keep track 03751 // of the current array index, and as a sanity check. 03752 size_t numLocalCount = 0; 03753 if (domainMap_->isContiguous ()) { 03754 // NOTE (mfh 03 Mar 2013, 02 Sep 2014) In the common case 03755 // that the domain Map is contiguous, it's more efficient to 03756 // avoid calling getNodeElementList(), since that 03757 // permanently constructs and caches the GID list in the 03758 // contiguous Map. 03759 GO curColMapGid = domainMap_->getMinGlobalIndex (); 03760 for (size_t i = 0; i < numDomainElts; ++i, ++curColMapGid) { 03761 if (GIDisLocal[i]) { 03762 LocalColGIDs[numLocalCount++] = curColMapGid; 03763 } 03764 } 03765 } 03766 else { 03767 ArrayView<const GO> domainElts = domainMap_->getNodeElementList (); 03768 for (size_t i = 0; i < numDomainElts; ++i) { 03769 if (GIDisLocal[i]) { 03770 LocalColGIDs[numLocalCount++] = domainElts[i]; 03771 } 03772 } 03773 } 03774 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 03775 numLocalCount != numLocalColGIDs, std::logic_error, 03776 ": numLocalCount = " << numLocalCount << " != numLocalColGIDs = " 03777 << numLocalColGIDs << ". This should never happen. Please report " 03778 "this bug to the Tpetra developers."); 03779 } 03780 03781 // FIXME (mfh 03 Apr 2013) Now would be a good time to use the 03782 // information we collected above to construct the Import. In 03783 // particular, building an Import requires: 03784 // 03785 // 1. numSameIDs (length of initial contiguous sequence of GIDs 03786 // on this process that are the same in both Maps; this 03787 // equals the number of domain Map elements on this process) 03788 // 03789 // 2. permuteToLIDs and permuteFromLIDs (both empty in this 03790 // case, since there's no permutation going on; the column 03791 // Map starts with the domain Map's GIDs, and immediately 03792 // after them come the remote GIDs) 03793 // 03794 // 3. remoteGIDs (exactly those GIDs that we found out above 03795 // were not in the domain Map) and remoteLIDs (which we could 03796 // have gotten above by using the three-argument version of 03797 // getRemoteIndexList() that computes local indices as well 03798 // as process ranks, instead of the two-argument version that 03799 // was used above) 03800 // 03801 // 4. remotePIDs (which we have from the getRemoteIndexList() 03802 // call above) 03803 // 03804 // 5. Sorting remotePIDs, and applying that permutation to 03805 // remoteGIDs and remoteLIDs (by calling sort3 above instead 03806 // of sort2) 03807 // 03808 // 6. Everything after the sort3 call in Import::setupExport(): 03809 // a. Create the Distributor via createFromRecvs(), which 03810 // computes exportGIDs and exportPIDs 03811 // b. Compute exportLIDs from exportGIDs (by asking the 03812 // source Map, in this case the domain Map, to convert 03813 // global to local) 03814 // 03815 // Steps 1-5 come for free, since we must do that work anyway in 03816 // order to compute the column Map. In particular, Step 3 is 03817 // even more expensive than Step 6a, since it involves both 03818 // creating and using a new Distributor object. 03819 03820 } // if the graph is globally indexed 03821 03822 const global_size_t gstInv = 03823 Teuchos::OrdinalTraits<global_size_t>::invalid (); 03824 // FIXME (mfh 05 Mar 2014) Doesn't the index base of a Map have to 03825 // be the same as the Map's min GID? If the first column is empty 03826 // (contains no entries), then the column Map's min GID won't 03827 // necessarily be the same as the domain Map's index base. 03828 const GO indexBase = domainMap_->getIndexBase (); 03829 colMap_ = rcp (new map_type (gstInv, myColumns, indexBase, 03830 domainMap_->getComm (), 03831 domainMap_->getNode ())); 03832 checkInternalState (); 03833 } 03834 03835 03836 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03837 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::mergeAllIndices() 03838 { 03839 TEUCHOS_TEST_FOR_EXCEPT( isGloballyIndexed() ); // call only after makeIndicesLocal() 03840 TEUCHOS_TEST_FOR_EXCEPT( ! isSorted() ); // call only after sortIndices() 03841 if (! isMerged ()) { 03842 for (size_t row=0; row < getNodeNumRows(); ++row) { 03843 RowInfo rowInfo = getRowInfo(row); 03844 mergeRowIndices(rowInfo); 03845 } 03846 // we just merged every row 03847 noRedundancies_ = true; 03848 } 03849 } 03850 03851 03852 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03853 void 03854 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeImportExport() 03855 { 03856 TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap (), std::logic_error, "Tpetra::" 03857 "CrsGraph: It's not allowed to call makeImportExport() unless the graph " 03858 "has a column Map."); 03859 RCP<ParameterList> params = this->getNonconstParameterList (); // could be null 03860 03861 // Don't do any checks to see if we need to create the Import, if 03862 // it exists already. 03863 // 03864 // FIXME (mfh 25 Mar 2013) This will become incorrect if we 03865 // change CrsGraph in the future to allow changing the column 03866 // Map after fillComplete. For now, the column Map is fixed 03867 // after the first fillComplete call. 03868 if (importer_.is_null ()) { 03869 // Create the Import instance if necessary. 03870 if (domainMap_ != colMap_ && (! domainMap_->isSameAs (*colMap_))) { 03871 if (params.is_null () || ! params->isSublist ("Import")) { 03872 importer_ = rcp (new import_type (domainMap_, colMap_)); 03873 } else { 03874 RCP<ParameterList> importSublist = sublist (params, "Import", true); 03875 importer_ = rcp (new import_type (domainMap_, colMap_, importSublist)); 03876 } 03877 } 03878 } 03879 03880 // Don't do any checks to see if we need to create the Export, if 03881 // it exists already. 03882 if (exporter_.is_null ()) { 03883 // Create the Export instance if necessary. 03884 if (rangeMap_ != rowMap_ && ! rangeMap_->isSameAs (*rowMap_)) { 03885 if (params.is_null () || ! params->isSublist ("Export")) { 03886 exporter_ = rcp (new export_type (rowMap_, rangeMap_)); 03887 } 03888 else { 03889 RCP<ParameterList> exportSublist = sublist (params, "Export", true); 03890 exporter_ = rcp (new export_type (rowMap_, rangeMap_, exportSublist)); 03891 } 03892 } 03893 } 03894 } 03895 03896 03897 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03898 std::string 03899 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::description() const 03900 { 03901 std::ostringstream oss; 03902 oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description(); 03903 if (isFillComplete()) { 03904 oss << "{status = fill complete" 03905 << ", global rows = " << getGlobalNumRows() 03906 << ", global cols = " << getGlobalNumCols() 03907 << ", global num entries = " << getGlobalNumEntries() 03908 << "}"; 03909 } 03910 else { 03911 oss << "{status = fill not complete" 03912 << ", global rows = " << getGlobalNumRows() 03913 << "}"; 03914 } 03915 return oss.str(); 03916 } 03917 03918 03919 template <class LocalOrdinal, class GlobalOrdinal, class Node> 03920 void 03921 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 03922 describe (Teuchos::FancyOStream &out, 03923 const Teuchos::EVerbosityLevel verbLevel) const 03924 { 03925 const char tfecfFuncName[] = "describe()"; 03926 using std::endl; 03927 using std::setw; 03928 using Teuchos::VERB_DEFAULT; 03929 using Teuchos::VERB_NONE; 03930 using Teuchos::VERB_LOW; 03931 using Teuchos::VERB_MEDIUM; 03932 using Teuchos::VERB_HIGH; 03933 using Teuchos::VERB_EXTREME; 03934 Teuchos::EVerbosityLevel vl = verbLevel; 03935 if (vl == VERB_DEFAULT) vl = VERB_LOW; 03936 RCP<const Comm<int> > comm = this->getComm(); 03937 const int myImageID = comm->getRank(), 03938 numImages = comm->getSize(); 03939 size_t width = 1; 03940 for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) { 03941 ++width; 03942 } 03943 width = std::max<size_t>(width,Teuchos::as<size_t>(11)) + 2; 03944 Teuchos::OSTab tab(out); 03945 // none: print nothing 03946 // low: print O(1) info from node 0 03947 // medium: print O(P) info, num entries per node 03948 // high: print O(N) info, num entries per row 03949 // extreme: print O(NNZ) info: print graph indices 03950 // 03951 // for medium and higher, print constituent objects at specified verbLevel 03952 if (vl != VERB_NONE) { 03953 if (myImageID == 0) out << this->description() << std::endl; 03954 // O(1) globals, minus what was already printed by description() 03955 if (isFillComplete() && myImageID == 0) { 03956 out << "Global number of diagonals = " << globalNumDiags_ << std::endl; 03957 out << "Global max number of row entries = " << globalMaxNumRowEntries_ << std::endl; 03958 } 03959 // constituent objects 03960 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) { 03961 if (myImageID == 0) out << "\nRow map: " << std::endl; 03962 rowMap_->describe(out,vl); 03963 if (colMap_ != null) { 03964 if (myImageID == 0) out << "\nColumn map: " << std::endl; 03965 colMap_->describe(out,vl); 03966 } 03967 if (domainMap_ != null) { 03968 if (myImageID == 0) out << "\nDomain map: " << std::endl; 03969 domainMap_->describe(out,vl); 03970 } 03971 if (rangeMap_ != null) { 03972 if (myImageID == 0) out << "\nRange map: " << std::endl; 03973 rangeMap_->describe(out,vl); 03974 } 03975 } 03976 // O(P) data 03977 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) { 03978 for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) { 03979 if (myImageID == imageCtr) { 03980 out << "Node ID = " << imageCtr << std::endl 03981 << "Node number of entries = " << nodeNumEntries_ << std::endl 03982 << "Node number of diagonals = " << nodeNumDiags_ << std::endl 03983 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl; 03984 if (isUpperTriangular()) { 03985 out << "Graph is upper triangular" << std::endl; 03986 } 03987 if (isLowerTriangular()) { 03988 out << "Graph is lower triangular" << std::endl; 03989 } 03990 if (indicesAreAllocated()) { 03991 out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl; 03992 } 03993 else { 03994 out << "Indices are not allocated." << std::endl; 03995 } 03996 } 03997 comm->barrier(); 03998 comm->barrier(); 03999 comm->barrier(); 04000 } 04001 } 04002 // O(N) and O(NNZ) data 04003 if (vl == VERB_HIGH || vl == VERB_EXTREME) { 04004 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": reduce verbosity level; graph row information was deleted at fillComplete()."); 04005 for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) { 04006 if (myImageID == imageCtr) { 04007 out << std::setw(width) << "Node ID" 04008 << std::setw(width) << "Global Row" 04009 << std::setw(width) << "Num Entries"; 04010 if (vl == VERB_EXTREME) { 04011 out << "Entries"; 04012 } 04013 out << std::endl; 04014 for (size_t r=0; r < getNodeNumRows(); ++r) { 04015 RowInfo rowinfo = getRowInfo(r); 04016 GlobalOrdinal gid = rowMap_->getGlobalElement(r); 04017 out << std::setw(width) << myImageID 04018 << std::setw(width) << gid 04019 << std::setw(width) << rowinfo.numEntries; 04020 if (vl == VERB_EXTREME) { 04021 if (isGloballyIndexed()) { 04022 ArrayView<const GlobalOrdinal> rowview = getGlobalView(rowinfo); 04023 for (size_t j=0; j < rowinfo.numEntries; ++j) out << rowview[j] << " "; 04024 } 04025 else if (isLocallyIndexed()) { 04026 ArrayView<const LocalOrdinal> rowview = getLocalView(rowinfo); 04027 for (size_t j=0; j < rowinfo.numEntries; ++j) out << colMap_->getGlobalElement(rowview[j]) << " "; 04028 } 04029 } 04030 out << std::endl; 04031 } 04032 } 04033 comm->barrier(); 04034 comm->barrier(); 04035 comm->barrier(); 04036 } 04037 } 04038 } 04039 } 04040 04041 04042 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04043 bool 04044 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04045 checkSizes (const SrcDistObject& source) 04046 { 04047 (void) source; // forestall "unused variable" compiler warnings 04048 04049 // It's not clear what kind of compatibility checks on sizes can 04050 // be performed here. Epetra_CrsGraph doesn't check any sizes for 04051 // compatibility. 04052 return true; 04053 } 04054 04055 04056 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04057 void 04058 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04059 copyAndPermute (const SrcDistObject& source, 04060 size_t numSameIDs, 04061 const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs, 04062 const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs) 04063 { 04064 using Teuchos::Array; 04065 using Teuchos::ArrayView; 04066 typedef LocalOrdinal LO; 04067 typedef GlobalOrdinal GO; 04068 const char tfecfFuncName[] = "copyAndPermute"; 04069 typedef CrsGraph<LO, GO, Node> this_type; 04070 typedef RowGraph<LO, GO, Node> row_graph_type; 04071 04072 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04073 permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error, 04074 ": permuteToLIDs and permuteFromLIDs must have the same size."); 04075 // Make sure that the source object has the right type. We only 04076 // actually need it to be a RowGraph, with matching first three 04077 // template parameters. If it's a CrsGraph, we can use view mode 04078 // instead of copy mode to get each row's data. 04079 // 04080 // FIXME (mfh 07 Jul 2013) It should not be necessary for any of 04081 // the template parameters but GO to match. GO has to match 04082 // because the graph has to send indices as global ordinals, if 04083 // the source and target graphs do not have the same column Map. 04084 // If LO doesn't match, the graphs could communicate using global 04085 // indices. It could be possible that Node affects the graph's 04086 // storage format, but packAndPrepare should assume a common 04087 // communication format in any case. 04088 const row_graph_type* srcRowGraph = dynamic_cast<const row_graph_type*> (&source); 04089 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04090 srcRowGraph == NULL, std::invalid_argument, 04091 ": The source object must be a RowGraph with matching first three " 04092 "template parameters."); 04093 04094 // If the source object is actually a CrsGraph, we can use view 04095 // mode instead of copy mode to access the entries in each row, 04096 // if the graph is not fill complete. 04097 const this_type* srcCrsGraph = dynamic_cast<const this_type*> (&source); 04098 04099 const map_type& srcRowMap = * (srcRowGraph->getRowMap ()); 04100 const map_type& tgtRowMap = * (this->getRowMap ()); 04101 const bool src_filled = srcRowGraph->isFillComplete (); 04102 Array<GO> row_copy; 04103 LO myid = 0; 04104 04105 // 04106 // "Copy" part of "copy and permute." 04107 // 04108 if (src_filled || srcCrsGraph == NULL) { 04109 // If the source graph is fill complete, we can't use view mode, 04110 // because the data might be stored in a different format not 04111 // compatible with the expectations of view mode. Also, if the 04112 // source graph is not a CrsGraph, we can't use view mode, 04113 // because RowGraph only provides copy mode access to the data. 04114 for (size_t i = 0; i < numSameIDs; ++i, ++myid) { 04115 const GO gid = srcRowMap.getGlobalElement (myid); 04116 size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (gid); 04117 row_copy.resize (row_length); 04118 size_t check_row_length = 0; 04119 srcRowGraph->getGlobalRowCopy (gid, row_copy (), check_row_length); 04120 this->insertGlobalIndices (gid, row_copy ()); 04121 } 04122 } else { 04123 for (size_t i = 0; i < numSameIDs; ++i, ++myid) { 04124 const GO gid = srcRowMap.getGlobalElement (myid); 04125 ArrayView<const GO> row; 04126 srcCrsGraph->getGlobalRowView (gid, row); 04127 this->insertGlobalIndices (gid, row); 04128 } 04129 } 04130 04131 // 04132 // "Permute" part of "copy and permute." 04133 // 04134 if (src_filled || srcCrsGraph == NULL) { 04135 for (LO i = 0; i < permuteToLIDs.size (); ++i) { 04136 const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]); 04137 const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]); 04138 size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (srcgid); 04139 row_copy.resize (row_length); 04140 size_t check_row_length = 0; 04141 srcRowGraph->getGlobalRowCopy (srcgid, row_copy (), check_row_length); 04142 this->insertGlobalIndices (mygid, row_copy ()); 04143 } 04144 } else { 04145 for (LO i = 0; i < permuteToLIDs.size (); ++i) { 04146 const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]); 04147 const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]); 04148 ArrayView<const GO> row; 04149 srcCrsGraph->getGlobalRowView (srcgid, row); 04150 this->insertGlobalIndices (mygid, row); 04151 } 04152 } 04153 } 04154 04155 04156 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04157 void 04158 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04159 packAndPrepare (const SrcDistObject& source, 04160 const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs, 04161 Teuchos::Array<GlobalOrdinal> &exports, 04162 const Teuchos::ArrayView<size_t> & numPacketsPerLID, 04163 size_t& constantNumPackets, 04164 Distributor &distor) 04165 { 04166 using Teuchos::Array; 04167 const char tfecfFuncName[] = "packAndPrepare"; 04168 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04169 exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error, 04170 ": exportLIDs and numPacketsPerLID must have the same size."); 04171 typedef RowGraph<LocalOrdinal, GlobalOrdinal, Node> row_graph_type; 04172 const row_graph_type& srcGraph = dynamic_cast<const row_graph_type&> (source); 04173 04174 // We don't check whether src_graph has had fillComplete called, 04175 // because it doesn't matter whether the *source* graph has been 04176 // fillComplete'd. The target graph can not be fillComplete'd yet. 04177 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04178 this->isFillComplete (), std::runtime_error, 04179 ": The target graph of an Import or Export must not be fill complete."); 04180 04181 srcGraph.pack (exportLIDs, exports, numPacketsPerLID, constantNumPackets, distor); 04182 } 04183 04184 04185 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04186 void 04187 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04188 pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 04189 Teuchos::Array<GlobalOrdinal>& exports, 04190 const Teuchos::ArrayView<size_t>& numPacketsPerLID, 04191 size_t& constantNumPackets, 04192 Distributor& distor) const 04193 { 04194 using Teuchos::Array; 04195 typedef LocalOrdinal LO; 04196 typedef GlobalOrdinal GO; 04197 const char tfecfFuncName[] = "pack"; 04198 (void) distor; // forestall "unused argument" compiler warning 04199 04200 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04201 exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error, 04202 ": exportLIDs and numPacketsPerLID must have the same size."); 04203 04204 const map_type& srcMap = * (this->getMap ()); 04205 constantNumPackets = 0; 04206 04207 // Set numPacketsPerLID[i] to the number of entries owned by the 04208 // calling process in (local) row exportLIDs[i] of the graph, that 04209 // the caller wants us to send out. Compute the total number of 04210 // packets (that is, entries) owned by this process in all the 04211 // rows that the caller wants us to send out. 04212 size_t totalNumPackets = 0; 04213 Array<GO> row; 04214 for (LO i = 0; i < exportLIDs.size (); ++i) { 04215 const GO GID = srcMap.getGlobalElement (exportLIDs[i]); 04216 size_t row_length = this->getNumEntriesInGlobalRow (GID); 04217 numPacketsPerLID[i] = row_length; 04218 totalNumPackets += row_length; 04219 } 04220 04221 exports.resize (totalNumPackets); 04222 04223 // Loop again over the rows to export, and pack rows of indices 04224 // into the output buffer. 04225 size_t exportsOffset = 0; 04226 for (LO i = 0; i < exportLIDs.size (); ++i) { 04227 const GO GID = srcMap.getGlobalElement (exportLIDs[i]); 04228 size_t row_length = this->getNumEntriesInGlobalRow (GID); 04229 row.resize (row_length); 04230 size_t check_row_length = 0; 04231 this->getGlobalRowCopy (GID, row (), check_row_length); 04232 typename Array<GO>::const_iterator row_iter = row.begin(); 04233 typename Array<GO>::const_iterator row_end = row.end(); 04234 size_t j = 0; 04235 for (; row_iter != row_end; ++row_iter, ++j) { 04236 exports[exportsOffset+j] = *row_iter; 04237 } 04238 exportsOffset += row.size(); 04239 } 04240 } 04241 04242 04243 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04244 void 04245 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04246 unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs, 04247 const Teuchos::ArrayView<const GlobalOrdinal> &imports, 04248 const Teuchos::ArrayView<size_t> &numPacketsPerLID, 04249 size_t constantNumPackets, 04250 Distributor& /* distor */, 04251 CombineMode /* CM */) 04252 { 04253 // FIXME (mfh 02 Apr 2012) REPLACE combine mode has a perfectly 04254 // reasonable meaning, whether or not the matrix is fill complete. 04255 // It's just more work to implement. 04256 04257 // We are not checking the value of the CombineMode input 04258 // argument. For CrsGraph, we only support import/export 04259 // operations if fillComplete has not yet been called. Any 04260 // incoming column-indices are inserted into the target graph. In 04261 // this context, CombineMode values of ADD vs INSERT are 04262 // equivalent. What is the meaning of REPLACE for CrsGraph? If a 04263 // duplicate column-index is inserted, it will be compressed out 04264 // when fillComplete is called. 04265 // 04266 // Note: I think REPLACE means that an existing row is replaced by 04267 // the imported row, i.e., the existing indices are cleared. CGB, 04268 // 6/17/2010 04269 04270 const char tfecfFuncName[] = "unpackAndCombine"; 04271 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04272 importLIDs.size() != numPacketsPerLID.size(), std::runtime_error, 04273 ": importLIDs and numPacketsPerLID must have the same size."); 04274 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 04275 isFillComplete (), std::runtime_error, 04276 ": Import or Export operations are not allowed on the destination " 04277 "CrsGraph if it is fill complete."); 04278 size_t importsOffset = 0; 04279 04280 typedef typename ArrayView<const LocalOrdinal>::const_iterator iter_type; 04281 iter_type impLIDiter = importLIDs.begin(); 04282 iter_type impLIDend = importLIDs.end(); 04283 04284 for (size_t i = 0; impLIDiter != impLIDend; ++impLIDiter, ++i) { 04285 LocalOrdinal row_length = numPacketsPerLID[i]; 04286 ArrayView<const GlobalOrdinal> row (&imports[importsOffset], row_length); 04287 insertGlobalIndicesFiltered (this->getMap ()->getGlobalElement (*impLIDiter), row); 04288 importsOffset += row_length; 04289 } 04290 } 04291 04292 04293 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04294 void 04295 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04296 removeEmptyProcessesInPlace (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& newMap) 04297 { 04298 using Teuchos::Comm; 04299 using Teuchos::null; 04300 using Teuchos::ParameterList; 04301 using Teuchos::RCP; 04302 04303 // We'll set all the state "transactionally," so that this method 04304 // satisfies the strong exception guarantee. This object's state 04305 // won't be modified until the end of this method. 04306 RCP<const map_type> rowMap, domainMap, rangeMap, colMap; 04307 RCP<import_type> importer; 04308 RCP<export_type> exporter; 04309 04310 rowMap = newMap; 04311 RCP<const Comm<int> > newComm = 04312 (newMap.is_null ()) ? null : newMap->getComm (); 04313 04314 if (! domainMap_.is_null ()) { 04315 if (domainMap_.getRawPtr () == rowMap_.getRawPtr ()) { 04316 // Common case: original domain and row Maps are identical. 04317 // In that case, we need only replace the original domain Map 04318 // with the new Map. This ensures that the new domain and row 04319 // Maps _stay_ identical. 04320 domainMap = newMap; 04321 } else { 04322 domainMap = domainMap_->replaceCommWithSubset (newComm); 04323 } 04324 } 04325 if (! rangeMap_.is_null ()) { 04326 if (rangeMap_.getRawPtr () == rowMap_.getRawPtr ()) { 04327 // Common case: original range and row Maps are identical. In 04328 // that case, we need only replace the original range Map with 04329 // the new Map. This ensures that the new range and row Maps 04330 // _stay_ identical. 04331 rangeMap = newMap; 04332 } else { 04333 rangeMap = rangeMap_->replaceCommWithSubset (newComm); 04334 } 04335 } 04336 if (! colMap.is_null ()) { 04337 colMap = colMap_->replaceCommWithSubset (newComm); 04338 } 04339 04340 // (Re)create the Export and / or Import if necessary. 04341 if (! newComm.is_null ()) { 04342 RCP<ParameterList> params = this->getNonconstParameterList (); // could be null 04343 // 04344 // The operations below are collective on the new communicator. 04345 // 04346 // (Re)create the Export object if necessary. If I haven't 04347 // called fillComplete yet, I don't have a rangeMap, so I must 04348 // first check if the _original_ rangeMap is not null. Ditto 04349 // for the Import object and the domain Map. 04350 if (! rangeMap_.is_null () && 04351 rangeMap != rowMap && 04352 ! rangeMap->isSameAs (*rowMap)) { 04353 if (params.is_null () || ! params->isSublist ("Export")) { 04354 exporter = rcp (new export_type (rowMap, rangeMap)); 04355 } 04356 else { 04357 RCP<ParameterList> exportSublist = sublist (params, "Export", true); 04358 exporter = rcp (new export_type (rowMap, rangeMap, exportSublist)); 04359 } 04360 } 04361 // (Re)create the Import object if necessary. 04362 if (! domainMap_.is_null () && 04363 domainMap != colMap && 04364 ! domainMap->isSameAs (*colMap)) { 04365 if (params.is_null () || ! params->isSublist ("Import")) { 04366 importer = rcp (new import_type (domainMap, colMap)); 04367 } else { 04368 RCP<ParameterList> importSublist = sublist (params, "Import", true); 04369 importer = rcp (new import_type (domainMap, colMap, importSublist)); 04370 } 04371 } 04372 } // if newComm is not null 04373 04374 // Defer side effects until the end. If no destructors throw 04375 // exceptions (they shouldn't anyway), then this method satisfies 04376 // the strong exception guarantee. 04377 exporter_ = exporter; 04378 importer_ = importer; 04379 rowMap_ = rowMap; 04380 // mfh 31 Mar 2013: DistObject's map_ is the row Map of a CrsGraph 04381 // or CrsMatrix. CrsGraph keeps a redundant pointer (rowMap_) to 04382 // the same object. We might want to get rid of this redundant 04383 // pointer sometime, but for now, we'll leave it alone and just 04384 // set map_ to the same object. 04385 this->map_ = rowMap; 04386 domainMap_ = domainMap; 04387 rangeMap_ = rangeMap; 04388 colMap_ = colMap; 04389 } 04390 04391 04392 template <class LocalOrdinal, class GlobalOrdinal, class Node> 04393 bool 04394 CrsGraph<LocalOrdinal,GlobalOrdinal,Node>:: 04395 hasRowInfo () const 04396 { 04397 if (indicesAreAllocated () && 04398 getProfileType () == StaticProfile && 04399 rowPtrs_.is_null ()) { 04400 return false; 04401 } else { 04402 return true; 04403 } 04404 } 04405 04406 } // namespace Tpetra 04407 04408 // Include KokkosRefactor partial specialisation if enabled 04409 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR) 04410 #include "Tpetra_KokkosRefactor_CrsGraph_def.hpp" 04411 #endif 04412 04413 04414 // 04415 // Explicit instantiation macro 04416 // 04417 // Must be expanded from within the Tpetra namespace! 04418 // 04419 #define TPETRA_CRSGRAPH_GRAPH_INSTANT(LO,GO,NODE) template class CrsGraph< LO , GO , NODE >; 04420 #define TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::sortRowIndicesAndValues< S >(const RowInfo, const Teuchos::ArrayView< S >& ); 04421 #define TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::mergeRowIndicesAndValues< S >(RowInfo, const ArrayView< S >& ); 04422 #define TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE) template ArrayRCP< S > CrsGraph< LO , GO , NODE >::allocateValues1D< S >() const; 04423 #define TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE) template ArrayRCP< Array< S > > CrsGraph< LO , GO , NODE >::allocateValues2D< S >() const; 04424 04425 #define TPETRA_CRSGRAPH_INSTANT(S,LO,GO,NODE) \ 04426 TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) \ 04427 TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) \ 04428 TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE) \ 04429 TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE) 04430 04431 #endif // TPETRA_CRSGRAPH_DEF_HPP
1.7.6.1
