|
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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00025 // 00026 // ************************************************************************ 00027 //@HEADER 00028 00029 // FINISH: need to check that fill is active before performing a number of the methods here; adding to the tests currently 00030 00031 #ifndef TPETRA_CRSGRAPH_DEF_HPP 00032 #define TPETRA_CRSGRAPH_DEF_HPP 00033 00034 #include <Kokkos_NodeTrace.hpp> 00035 #include <Teuchos_TestForException.hpp> 00036 #include <Teuchos_NullIteratorTraits.hpp> 00037 #include <algorithm> 00038 #include <string> 00039 00040 #ifdef DOXYGEN_USE_ONLY 00041 #include "Tpetra_CrsGraph_decl.hpp" 00042 #endif 00043 00044 namespace Tpetra { 00045 00046 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00047 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsGraph(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00048 size_t maxNumEntriesPerRow, ProfileType pftype) 00049 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00050 , rowMap_(rowMap) 00051 , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode()) 00052 , nodeNumEntries_(0) 00053 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00054 , pftype_(pftype) 00055 , numAllocForAllRows_(maxNumEntriesPerRow) 00056 , indicesAreAllocated_(false) 00057 , indicesAreLocal_(false) 00058 , indicesAreGlobal_(false) 00059 { 00060 staticAssertions(); 00061 std::string tfecfFuncName("CrsGraph(rowMap,maxNumEntriesPerRow)"); 00062 TEST_FOR_EXCEPTION_CLASS_FUNC(maxNumEntriesPerRow > OrdinalTraits<size_t>::max() || (maxNumEntriesPerRow < 1 && maxNumEntriesPerRow != 0), std::runtime_error, 00063 ": maxNumEntriesPerRow must be non-negative."); 00064 resumeFill(); 00065 checkInternalState(); 00066 } 00067 00068 00069 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00070 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsGraph(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00071 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 00072 size_t maxNumEntriesPerRow, ProfileType pftype) 00073 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00074 , rowMap_(rowMap) 00075 , colMap_(colMap) 00076 , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode()) 00077 , nodeNumEntries_(0) 00078 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00079 , pftype_(pftype) 00080 , numAllocForAllRows_(maxNumEntriesPerRow) 00081 , indicesAreAllocated_(false) 00082 , indicesAreLocal_(false) 00083 , indicesAreGlobal_(false) 00084 { 00085 staticAssertions(); 00086 std::string tfecfFuncName("CrsGraph(rowMap,colMap,maxNumEntriesPerRow)"); 00087 TEST_FOR_EXCEPTION_CLASS_FUNC(maxNumEntriesPerRow > OrdinalTraits<size_t>::max() || (maxNumEntriesPerRow < 1 && maxNumEntriesPerRow != 0), std::runtime_error, 00088 ": maxNumEntriesPerRow must be non-negative."); 00089 resumeFill(); 00090 checkInternalState(); 00091 } 00092 00093 00094 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00095 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsGraph(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00096 const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype) 00097 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00098 , rowMap_(rowMap) 00099 , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode()) 00100 , nodeNumEntries_(0) 00101 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00102 , pftype_(pftype) 00103 , numAllocPerRow_(NumEntriesPerRowToAlloc) 00104 , numAllocForAllRows_(0) 00105 , indicesAreAllocated_(false) 00106 , indicesAreLocal_(false) 00107 , indicesAreGlobal_(false) 00108 { 00109 staticAssertions(); 00110 std::string tfecfFuncName("CrsGraph(rowMap,NumEntriesPerRowToAlloc)"); 00111 TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::runtime_error, 00112 ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node."); 00113 size_t numMin = OrdinalTraits<size_t>::max(), 00114 numMax = OrdinalTraits<size_t>::zero(); 00115 for (size_t r=0; r < getNodeNumRows(); ++r) { 00116 numMin = std::min<size_t>( numMin, NumEntriesPerRowToAlloc[r] ); 00117 numMax = std::max<size_t>( numMax, NumEntriesPerRowToAlloc[r] ); 00118 } 00119 TEST_FOR_EXCEPTION_CLASS_FUNC((numMin < OrdinalTraits<size_t>::one() && numMin != OrdinalTraits<size_t>::zero()) || numMax > OrdinalTraits<size_t>::max(), std::runtime_error, ": Invalid user-specified number of non-zeros per row."); 00120 resumeFill(); 00121 checkInternalState(); 00122 } 00123 00124 00125 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00126 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsGraph(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 00127 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 00128 const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype) 00129 : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap) 00130 , rowMap_(rowMap) 00131 , colMap_(colMap) 00132 , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode()) 00133 , nodeNumEntries_(0) 00134 , nodeNumAllocated_(OrdinalTraits<size_t>::invalid()) 00135 , pftype_(pftype) 00136 , numAllocPerRow_(NumEntriesPerRowToAlloc) 00137 , numAllocForAllRows_(0) 00138 , indicesAreAllocated_(false) 00139 , indicesAreLocal_(false) 00140 , indicesAreGlobal_(false) 00141 { 00142 staticAssertions(); 00143 std::string tfecfFuncName("CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc)"); 00144 TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::runtime_error, 00145 ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node."); 00146 resumeFill(); 00147 checkInternalState(); 00148 } 00149 00150 00151 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00152 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~CrsGraph() 00153 {} 00154 00155 00156 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00157 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumRows() const 00158 { 00159 return rowMap_->getGlobalNumElements(); 00160 } 00161 00162 00163 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00164 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumCols() const 00165 { 00166 std::string tfecfFuncName("getGlobalNumCols()"); 00167 TEST_FOR_EXCEPTION(hasColMap() == false, std::runtime_error, ": requires column map."); 00168 return colMap_->getGlobalNumElements(); 00169 } 00170 00171 00172 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00173 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumRows() const 00174 { 00175 return rowMap_->getNodeNumElements(); 00176 } 00177 00178 00179 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00180 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumCols() const 00181 { 00182 std::string tfecfFuncName("getNodeNumCols()"); 00183 TEST_FOR_EXCEPTION_CLASS_FUNC(hasColMap() != true, std::runtime_error, ": requires column map."); 00184 return colMap_->getNodeNumElements(); 00185 } 00186 00187 00188 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00189 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumDiags() const 00190 { 00191 return nodeNumDiags_; 00192 } 00193 00194 00195 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00196 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumDiags() const 00197 { 00198 return globalNumDiags_; 00199 } 00200 00201 00202 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00203 RCP<Node> 00204 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const 00205 { 00206 return rowMap_->getNode(); 00207 } 00208 00209 00210 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00211 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 00212 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowMap() const 00213 { 00214 return rowMap_; 00215 } 00216 00217 00218 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00219 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 00220 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getColMap() const 00221 { 00222 return colMap_; 00223 } 00224 00225 00226 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00227 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 00228 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const 00229 { 00230 return domainMap_; 00231 } 00232 00233 00234 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00235 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 00236 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const 00237 { 00238 return rangeMap_; 00239 } 00240 00241 00242 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00243 RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > 00244 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getImporter() const 00245 { 00246 return importer_; 00247 } 00248 00249 00250 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00251 RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > 00252 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getExporter() const 00253 { 00254 return exporter_; 00255 } 00256 00257 00258 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00259 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasColMap() const 00260 { 00261 return (colMap_ != null); 00262 } 00263 00264 00265 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00266 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStorageOptimized() const 00267 { 00268 bool isOpt = lclGraph_.isOptimized(); 00269 #ifdef HAVE_TPETRA_DEBUG 00270 std::string tfecfFuncName("isStorageOptimized()"); 00271 TEST_FOR_EXCEPTION_CLASS_FUNC( (isOpt == true) && (getProfileType() == DynamicProfile), std::logic_error, 00272 ": Violated stated post-conditions. Please contact Tpetra team."); 00273 #endif 00274 return isOpt; 00275 } 00276 00277 00278 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00279 ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getProfileType() const 00280 { 00281 return pftype_; 00282 } 00283 00284 00285 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00286 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumEntries() const 00287 { 00288 return globalNumEntries_; 00289 } 00290 00291 00292 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00293 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumEntries() const 00294 { 00295 return nodeNumEntries_; 00296 } 00297 00298 00299 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00300 global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalMaxNumRowEntries() const 00301 { 00302 return globalMaxNumRowEntries_; 00303 } 00304 00305 00306 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00307 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeMaxNumRowEntries() const 00308 { 00309 return nodeMaxNumRowEntries_; 00310 } 00311 00312 00313 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00314 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const 00315 { 00316 return fillComplete_; 00317 } 00318 00319 00320 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00321 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillActive() const 00322 { 00323 return !fillComplete_; 00324 } 00325 00326 00327 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00328 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isUpperTriangular() const 00329 { 00330 return upperTriangular_; 00331 } 00332 00333 00334 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00335 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLowerTriangular() const 00336 { 00337 return lowerTriangular_; 00338 } 00339 00340 00341 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00342 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLocallyIndexed() const 00343 { 00344 return indicesAreLocal_; 00345 } 00346 00347 00348 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00349 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isGloballyIndexed() const 00350 { 00351 return indicesAreGlobal_; 00352 } 00353 00354 00355 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00356 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeAllocationSize() const 00357 { 00358 return nodeNumAllocated_; 00359 } 00360 00361 00362 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00363 const RCP<const Comm<int> > & 00364 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getComm() const 00365 { 00366 return rowMap_->getComm(); 00367 } 00368 00369 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00370 GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getIndexBase() const 00371 { 00372 return rowMap_->getIndexBase(); 00373 } 00374 00375 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00376 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::indicesAreAllocated() const 00377 { 00378 return indicesAreAllocated_; 00379 } 00380 00381 00384 // // 00385 // Internal utility methods // 00386 // // 00389 00390 00393 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00394 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isSorted() const 00395 { 00396 return indicesAreSorted_; 00397 } 00398 00399 00402 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00403 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isMerged() const 00404 { 00405 return noRedundancies_; 00406 } 00407 00408 00411 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00412 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setSorted(bool sorted) 00413 { 00414 indicesAreSorted_ = sorted; 00415 } 00416 00417 00420 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00421 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setMerged(bool merged) 00422 { 00423 noRedundancies_ = merged; 00424 } 00425 00426 00429 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00430 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateIndices(ELocalGlobal lg) 00431 { 00432 // this is a protected function, only callable by us. if it was called incorrectly, it is our fault. 00433 std::string tfecfFuncName("allocateIndices()"); 00434 TEST_FOR_EXCEPTION_CLASS_FUNC( isLocallyIndexed() && lg==GlobalIndices, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 00435 TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() && lg==LocalIndices, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 00436 TEST_FOR_EXCEPTION_CLASS_FUNC( indicesAreAllocated() == true, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 00437 const size_t numRows = getNodeNumRows(); 00438 indicesAreLocal_ = (lg == LocalIndices); 00439 indicesAreGlobal_ = (lg == GlobalIndices); 00440 nodeNumAllocated_ = 0; 00441 if (getProfileType() == StaticProfile) { 00442 // 00443 // STATIC ALLOCATION PROFILE 00444 // 00445 // determine how many entries to allocate and setup offsets into 1D arrays 00446 rowBegs_ = arcp<size_t>(numRows+1); 00447 rowBegs_[numRows] = 0; 00448 if (numRows > 0) { 00449 rowEnds_ = arcp<size_t>(numRows); 00450 // allocate offsets 00451 if (numAllocPerRow_ != null) { 00452 nodeNumAllocated_ = 0; 00453 for (size_t i=0; i < numRows; ++i) { 00454 rowBegs_[i] = nodeNumAllocated_; 00455 rowEnds_[i] = rowBegs_[i]; 00456 nodeNumAllocated_ += numAllocPerRow_[i]; 00457 } 00458 rowBegs_[numRows] = nodeNumAllocated_; 00459 } 00460 else { 00461 nodeNumAllocated_ = numAllocForAllRows_ * numRows; 00462 rowBegs_[0] = 0; 00463 for (size_t i=0; i < numRows; ++i) { 00464 rowEnds_[i] = rowBegs_[i]; 00465 rowBegs_[i+1] = rowBegs_[i] + numAllocForAllRows_; 00466 } 00467 } 00468 // allocate the indices 00469 if (nodeNumAllocated_ > 0) { 00470 if (lg == LocalIndices) { 00471 lclInds1D_ = arcp<LocalOrdinal>(nodeNumAllocated_); 00472 } 00473 else { 00474 gblInds1D_ = arcp<GlobalOrdinal>(nodeNumAllocated_); 00475 } 00476 } 00477 } 00478 } 00479 else { 00480 if (numRows > 0) { 00481 numEntriesPerRow_ = arcp<size_t>(numRows); 00482 std::fill( numEntriesPerRow_.begin(), numEntriesPerRow_.end(), (size_t)0 ); 00483 // 00484 // DYNAMIC ALLOCATION PROFILE 00485 // 00486 ArrayRCP<const size_t>::iterator numalloc = numAllocPerRow_.begin(); 00487 size_t howmany = numAllocForAllRows_; 00488 if (lg == LocalIndices) { 00489 lclInds2D_ = arcp< ArrayRCP<LocalOrdinal> >(numRows); 00490 nodeNumAllocated_ = 0; 00491 for (size_t i=0; i < numRows; ++i) { 00492 if (numAllocPerRow_ != null) howmany = *numalloc++; 00493 nodeNumAllocated_ += howmany; 00494 if (howmany > 0) lclInds2D_[i] = arcp<LocalOrdinal>(howmany); 00495 } 00496 } 00497 else { // allocate global indices 00498 gblInds2D_ = arcp< ArrayRCP<GlobalOrdinal> >(numRows); 00499 nodeNumAllocated_ = 0; 00500 for (size_t i=0; i < numRows; ++i) { 00501 if (numAllocPerRow_ != null) howmany = *numalloc++; 00502 nodeNumAllocated_ += howmany; 00503 if (howmany > 0) gblInds2D_[i] = arcp<GlobalOrdinal>(howmany); 00504 } 00505 } 00506 } 00507 } // if numRows > 0 00508 // done with these 00509 numAllocForAllRows_ = 0; 00510 numAllocPerRow_ = null; 00511 indicesAreAllocated_ = true; 00512 checkInternalState(); 00513 } 00514 00515 00518 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00519 template <class T> 00520 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues(ArrayRCP<T> &values1D, ArrayRCP<ArrayRCP<T> > &values2D) const 00521 { 00522 std::string tfecfFuncName("allocateValues()"); 00523 TEST_FOR_EXCEPTION_CLASS_FUNC( indicesAreAllocated() == false, std::runtime_error, ": graph indices must be allocated before values."); 00524 values1D = null; 00525 values2D = null; 00526 if (getProfileType() == StaticProfile) { 00527 if (lclInds1D_ != null) { 00528 values1D = arcp<T>(lclInds1D_.size()); 00529 } 00530 else if (gblInds1D_ != null) { 00531 values1D = arcp<T>(gblInds1D_.size()); 00532 } 00533 } 00534 else { 00535 values2D = arcp<ArrayRCP<T> >(getNodeNumRows()); 00536 if (lclInds2D_ != null) { 00537 for (size_t r=0; r < (size_t)lclInds2D_.size(); ++r) { 00538 if (lclInds2D_[r] != null) { 00539 values2D[r] = arcp<T>(lclInds2D_[r].size()); 00540 } 00541 } 00542 } 00543 else if (gblInds2D_ != null) { 00544 for (size_t r=0; r < (size_t)gblInds2D_.size(); ++r) { 00545 if (gblInds2D_[r] != null) { 00546 values2D[r] = arcp<T>(gblInds2D_[r].size()); 00547 } 00548 } 00549 } 00550 } 00551 } 00552 00553 00556 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00557 template <ELocalGlobal lg> 00558 RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateAlloc(RowInfo rowinfo, size_t newAllocSize) 00559 { 00560 #ifdef HAVE_TPETRA_DEBUG 00561 TEST_FOR_EXCEPT( rowMap_->isNodeLocalElement(rowinfo.localRow) == false ); 00562 TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize ); 00563 TEST_FOR_EXCEPT( (lg == LocalIndices && isLocallyIndexed() == false) || (lg == GlobalIndices && isGloballyIndexed() == false) ); 00564 TEST_FOR_EXCEPT( newAllocSize == 0 ); 00565 TEST_FOR_EXCEPT( indicesAreAllocated() == false ); 00566 #endif 00567 // allocate a larger space for row "lrow" 00568 // copy any existing data from previous allocation to new allocation 00569 // update sizes 00570 if (lg == LocalIndices) { 00571 ArrayRCP<LocalOrdinal> old_alloc = lclInds2D_[rowinfo.localRow]; 00572 lclInds2D_[rowinfo.localRow] = arcp<LocalOrdinal>(newAllocSize); 00573 std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, lclInds2D_[rowinfo.localRow].begin()); 00574 } 00575 else /* if lg == GlobalIndices */ { 00576 ArrayRCP<GlobalOrdinal> old_alloc = gblInds2D_[rowinfo.localRow]; 00577 gblInds2D_[rowinfo.localRow] = arcp<GlobalOrdinal>(newAllocSize); 00578 std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, gblInds2D_[rowinfo.localRow].begin()); 00579 } 00580 // 00581 nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize); 00582 rowinfo.allocSize = newAllocSize; 00583 return rowinfo; 00584 } 00585 00586 00589 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00590 template <ELocalGlobal lg, class T> 00591 RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateAllocAndValues(RowInfo rowinfo, size_t newAllocSize, ArrayRCP<T> &rowVals) 00592 { 00593 #ifdef HAVE_TPETRA_DEBUG 00594 TEST_FOR_EXCEPT( rowMap_->isNodeLocalElement(rowinfo.localRow) == false ); 00595 TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize ); 00596 TEST_FOR_EXCEPT( (lg == LocalIndices && isLocallyIndexed() == false) || (lg == GlobalIndices && isGloballyIndexed() == false) ); 00597 TEST_FOR_EXCEPT( newAllocSize == 0 ); 00598 TEST_FOR_EXCEPT( indicesAreAllocated() == false ); 00599 #endif 00600 // allocate a larger space for row "lrow" 00601 // copy any existing data from previous allocation to new allocation 00602 // update sizes 00603 if (lg == LocalIndices) { 00604 ArrayRCP<LocalOrdinal> old_alloc = lclInds2D_[rowinfo.localRow]; 00605 lclInds2D_[rowinfo.localRow] = arcp<LocalOrdinal>(newAllocSize); 00606 std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, lclInds2D_[rowinfo.localRow].begin()); 00607 } 00608 else /* if lg == GlobalIndices */ { 00609 ArrayRCP<GlobalOrdinal> old_alloc = gblInds2D_[rowinfo.localRow]; 00610 gblInds2D_[rowinfo.localRow] = arcp<GlobalOrdinal>(newAllocSize); 00611 std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, gblInds2D_[rowinfo.localRow].begin()); 00612 } 00613 ArrayRCP<const T> oldVals = rowVals; 00614 rowVals = arcp<T>(newAllocSize); 00615 std::copy(oldVals.begin(), oldVals.begin() + rowinfo.numEntries, rowVals.begin()); 00616 // 00617 nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize); 00618 rowinfo.allocSize = newAllocSize; 00619 return rowinfo; 00620 } 00621 00622 00625 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00626 ArrayView<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalView(RowInfo rowinfo) const 00627 { 00628 ArrayView<const LocalOrdinal> view; 00629 if (rowinfo.allocSize > 0) { 00630 if (lclInds1D_ != null) { 00631 view = lclInds1D_(rowinfo.offset1D,rowinfo.allocSize); 00632 } 00633 else if (lclInds2D_[rowinfo.localRow] != null) { 00634 view = lclInds2D_[rowinfo.localRow](); 00635 } 00636 } 00637 return view; 00638 } 00639 00640 00643 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00644 ArrayView<LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalViewNonConst(RowInfo rowinfo) 00645 { 00646 ArrayView<LocalOrdinal> view; 00647 if (rowinfo.allocSize > 0) { 00648 if (lclInds1D_ != null) { 00649 view = lclInds1D_(rowinfo.offset1D,rowinfo.allocSize); 00650 } 00651 else if (lclInds2D_[rowinfo.localRow] != null) { 00652 view = lclInds2D_[rowinfo.localRow](); 00653 } 00654 } 00655 return view; 00656 } 00657 00658 00661 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00662 ArrayView<const GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalView(RowInfo rowinfo) const 00663 { 00664 ArrayView<const GlobalOrdinal> view; 00665 if (rowinfo.allocSize > 0) { 00666 if (gblInds1D_ != null) { 00667 view = gblInds1D_(rowinfo.offset1D,rowinfo.allocSize); 00668 } 00669 else if (gblInds2D_[rowinfo.localRow] != null) { 00670 view = gblInds2D_[rowinfo.localRow](); 00671 } 00672 } 00673 return view; 00674 } 00675 00676 00679 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00680 ArrayView<GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalViewNonConst(RowInfo rowinfo) 00681 { 00682 ArrayView<GlobalOrdinal> view; 00683 if (rowinfo.allocSize > 0) { 00684 if (gblInds1D_ != null) { 00685 view = gblInds1D_(rowinfo.offset1D,rowinfo.allocSize); 00686 } 00687 else if (gblInds2D_[rowinfo.localRow] != null) { 00688 view = gblInds2D_[rowinfo.localRow](); 00689 } 00690 } 00691 return view; 00692 } 00693 00694 00697 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00698 RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowInfo(size_t myRow) const 00699 { 00700 #ifdef HAVE_TPETRA_DEBUG 00701 std::string tfecfFuncName("getRowInfo()"); 00702 TEST_FOR_EXCEPTION_CLASS_FUNC(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 00703 #endif 00704 const size_t STINV = OrdinalTraits<size_t>::invalid(); 00705 RowInfo ret; 00706 ret.localRow = myRow; 00707 if (nodeNumAllocated_ != 0 && nodeNumAllocated_ != STINV) { 00708 // graph data structures have the info that we need 00709 // 00710 // if static graph, offsets tell us the allocation size 00711 if (getProfileType() == StaticProfile) { 00712 ret.offset1D = rowBegs_[myRow]; 00713 ret.numEntries = rowEnds_[myRow] - rowBegs_[myRow]; 00714 ret.allocSize = rowBegs_[myRow+1] - rowBegs_[myRow]; 00715 } 00716 else { 00717 ret.offset1D = STINV; 00718 if (isLocallyIndexed()) { 00719 ret.allocSize = lclInds2D_[myRow].size(); 00720 } 00721 else { 00722 ret.allocSize = gblInds2D_[myRow].size(); 00723 } 00724 ret.numEntries = numEntriesPerRow_[myRow]; 00725 } 00726 } 00727 else if (nodeNumAllocated_ == 0) { 00728 // have performed allocation, but the graph has no allocation or entries 00729 ret.allocSize = 0; 00730 ret.numEntries = 0; 00731 ret.offset1D = STINV; 00732 } 00733 else if (indicesAreAllocated() == false) { 00734 // haven't performed allocation yet; probably won't hit this code 00735 if (numAllocPerRow_ == null) { 00736 ret.allocSize = numAllocForAllRows_; 00737 } 00738 else { 00739 ret.allocSize = numAllocPerRow_[myRow]; 00740 } 00741 ret.numEntries = 0; 00742 ret.offset1D = STINV; 00743 } 00744 else { 00745 // don't know how we ended up here... 00746 TEST_FOR_EXCEPT(true); 00747 } 00748 return ret; 00749 } 00750 00751 00754 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00755 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::staticAssertions() const 00756 { 00757 // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal) 00758 // This is so that we can store LocalOrdinals in the memory formerly occupied by GlobalOrdinals 00759 // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal) and max(size_t) >= max(LocalOrdinal) 00760 // This is so that we can represent any LocalOrdinal as a size_t, and any LocalOrdinal as a GlobalOrdinal 00761 Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size1; (void)cta_size1; 00762 Teuchos::CompileTimeAssert<sizeof(global_size_t) < sizeof(size_t) > cta_size2; (void)cta_size2; 00763 // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check 00764 std::string msg = typeName(*this) + ": Object cannot be allocated with stated template arguments: size assumptions are not valid."; 00765 TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(), std::runtime_error, msg); 00766 TEST_FOR_EXCEPTION( (global_size_t)OrdinalTraits<LocalOrdinal>::max() > (global_size_t)OrdinalTraits<GlobalOrdinal>::max(), std::runtime_error, msg); 00767 TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(), std::runtime_error, msg); 00768 TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(), std::runtime_error, msg); 00769 } 00770 00771 00774 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00775 template <ELocalGlobal lg> 00776 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndices(const SLocalGlobalNCViews &inds) const 00777 { 00778 const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *colMap_; 00779 Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; 00780 (void)cta_lg; 00781 size_t numFiltered = 0; 00782 #ifdef HAVE_TPETRA_DEBUG 00783 size_t numFiltered_debug = 0; 00784 #endif 00785 if (lg == GlobalIndices) { 00786 ArrayView<GlobalOrdinal> ginds = inds.ginds; 00787 typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(), 00788 cptr = ginds.begin(); 00789 while (cptr != ginds.end()) { 00790 if (cmap.isNodeGlobalElement(*cptr)) { 00791 *fend++ = *cptr; 00792 #ifdef HAVE_TPETRA_DEBUG 00793 ++numFiltered_debug; 00794 #endif 00795 } 00796 ++cptr; 00797 } 00798 numFiltered = fend - ginds.begin(); 00799 } 00800 else if (lg == LocalIndices) { 00801 ArrayView<LocalOrdinal> linds = inds.linds; 00802 typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(), 00803 cptr = linds.begin(); 00804 while (cptr != linds.end()) { 00805 if (cmap.isNodeLocalElement(*cptr)) { 00806 *fend++ = *cptr; 00807 #ifdef HAVE_TPETRA_DEBUG 00808 ++numFiltered_debug; 00809 #endif 00810 } 00811 ++cptr; 00812 } 00813 numFiltered = fend - linds.begin(); 00814 } 00815 #ifdef HAVE_TPETRA_DEBUG 00816 TEST_FOR_EXCEPT( numFiltered != numFiltered_debug ); 00817 #endif 00818 return numFiltered; 00819 } 00820 00821 00824 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00825 template <ELocalGlobal lg, class T> 00826 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndicesAndValues(const SLocalGlobalNCViews &inds, const ArrayView<T> &vals) const 00827 { 00828 const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *colMap_; 00829 Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; 00830 (void)cta_lg; 00831 size_t numFiltered = 0; 00832 typename ArrayView<T>::iterator fvalsend = vals.begin(), 00833 valscptr = vals.begin(); 00834 #ifdef HAVE_TPETRA_DEBUG 00835 size_t numFiltered_debug = 0; 00836 #endif 00837 if (lg == GlobalIndices) { 00838 ArrayView<GlobalOrdinal> ginds = inds.ginds; 00839 typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(), 00840 cptr = ginds.begin(); 00841 while (cptr != ginds.end()) { 00842 if (cmap.isNodeGlobalElement(*cptr)) { 00843 *fend++ = *cptr; 00844 *fvalsend++ = *valscptr; 00845 #ifdef HAVE_TPETRA_DEBUG 00846 ++numFiltered_debug; 00847 #endif 00848 } 00849 ++cptr; 00850 ++valscptr; 00851 } 00852 numFiltered = fend - ginds.begin(); 00853 } 00854 else if (lg == LocalIndices) { 00855 ArrayView<LocalOrdinal> linds = inds.linds; 00856 typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(), 00857 cptr = linds.begin(); 00858 while (cptr != linds.end()) { 00859 if (cmap.isNodeLocalElement(*cptr)) { 00860 *fend++ = *cptr; 00861 *fvalsend++ = *valscptr; 00862 #ifdef HAVE_TPETRA_DEBUG 00863 ++numFiltered_debug; 00864 #endif 00865 } 00866 ++cptr; 00867 ++valscptr; 00868 } 00869 numFiltered = fend - linds.begin(); 00870 } 00871 #ifdef HAVE_TPETRA_DEBUG 00872 TEST_FOR_EXCEPT( numFiltered != numFiltered_debug ); 00873 TEST_FOR_EXCEPT( valscptr != vals.end() ); 00874 TEST_FOR_EXCEPT( numFiltered != (size_t)(fvalsend - vals.begin()) ); 00875 #endif 00876 return numFiltered; 00877 } 00878 00879 00882 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00883 template <ELocalGlobal lg> 00884 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertIndices(RowInfo rowinfo, const SLocalGlobalViews &newInds) 00885 { 00886 Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; (void)cta_lg; 00887 size_t numNewInds = 0; 00888 if (lg == GlobalIndices) { 00889 ArrayView<const GlobalOrdinal> new_ginds = newInds.ginds; 00890 ArrayView<GlobalOrdinal> gind_view = getGlobalViewNonConst(rowinfo); 00891 numNewInds = new_ginds.size(); 00892 std::copy(new_ginds.begin(), new_ginds.end(), gind_view.begin()+rowinfo.numEntries); 00893 } 00894 else if (lg == LocalIndices) { 00895 ArrayView<const LocalOrdinal> new_linds = newInds.linds; 00896 ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo); 00897 numNewInds = new_linds.size(); 00898 std::copy(new_linds.begin(), new_linds.end(), lind_view.begin()+rowinfo.numEntries); 00899 } 00900 if (getProfileType() == StaticProfile) { 00901 rowEnds_[rowinfo.localRow] += numNewInds; 00902 } 00903 else { 00904 numEntriesPerRow_[rowinfo.localRow] += numNewInds; 00905 } 00906 nodeNumEntries_ += numNewInds; 00907 indicesAreSorted_ = false; 00908 noRedundancies_ = false; 00909 return numNewInds; 00910 } 00911 00912 00915 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00916 template <ELocalGlobal lg, class IterO, class IterN> 00917 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertIndicesAndValues(RowInfo rowinfo, const SLocalGlobalViews &newInds, IterO rowVals, IterN newVals) 00918 { 00919 size_t numNewInds = insertIndices<lg>(rowinfo,newInds); 00920 std::copy( newVals, newVals + numNewInds, rowVals + rowinfo.numEntries ); 00921 } 00922 00923 00926 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00927 template <ELocalGlobal lg, class IterO, class IterN, class BinaryFunction> 00928 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::transformValues(RowInfo rowinfo, const SLocalGlobalViews &inds, IterO rowVals, IterN newVals, BinaryFunction f) const 00929 { 00930 Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; (void)cta_lg; 00931 const size_t STINV = OrdinalTraits<size_t>::invalid(); 00932 if (lg == GlobalIndices) { 00933 ArrayView<const GlobalOrdinal> search_ginds = inds.ginds; 00934 for (size_t j=0; j < (size_t)search_ginds.size(); ++j) { 00935 const size_t k = findGlobalIndex(rowinfo, search_ginds[j]); 00936 if (k != STINV) { 00937 rowVals[k] = f( rowVals[k], newVals[j] ); 00938 } 00939 } 00940 } 00941 else if (lg == LocalIndices) { 00942 ArrayView<const LocalOrdinal> search_linds = inds.linds; 00943 for (size_t j=0; j < (size_t)search_linds.size(); ++j) { 00944 const size_t k = findLocalIndex(rowinfo, search_linds[j]); 00945 if (k != STINV) { 00946 rowVals[k] = f( rowVals[k], newVals[j] ); 00947 } 00948 } 00949 } 00950 } 00951 00952 00955 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00956 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndices(RowInfo rowinfo) 00957 { 00958 if (rowinfo.numEntries > 0) { 00959 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo); 00960 std::sort(inds_view.begin(), inds_view.begin() + rowinfo.numEntries); 00961 } 00962 } 00963 00964 00967 // in the future, perhaps this could use std::sort with a boost::zip_iterator 00968 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00969 template <class Iter> 00970 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndicesAndValues(RowInfo rowinfo, Iter rowValueIters) 00971 { 00972 if (rowinfo.numEntries > 0) { 00973 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo); 00974 sort2(inds_view.begin(), inds_view.begin() + rowinfo.numEntries, rowValueIters); 00975 } 00976 } 00977 00978 00981 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 00982 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRowIndices(RowInfo rowinfo) 00983 { 00984 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo); 00985 typename ArrayView<LocalOrdinal>::iterator beg, end, newend; 00986 beg = inds_view.begin(); 00987 end = inds_view.begin() + rowinfo.numEntries; 00988 newend = std::unique(beg,end); 00989 const size_t mergedEntries = newend - beg; 00990 #ifdef HAVE_TPETRA_DEBUG 00991 // merge should not have eliminated any entries; if so, the assignment below will destory the packed structure 00992 TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries ); 00993 #endif 00994 if (getProfileType() == StaticProfile) { 00995 rowEnds_[rowinfo.localRow] = rowBegs_[rowinfo.localRow] + mergedEntries; 00996 } 00997 else { 00998 numEntriesPerRow_[rowinfo.localRow] = mergedEntries; 00999 } 01000 nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries); 01001 } 01002 01003 01006 // in the future, this could use std::unique with a boost::zip_iterator 01007 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01008 template <class Iter, class BinaryFunction> 01009 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRowIndicesAndValues(RowInfo rowinfo, Iter rowValueIter, BinaryFunction f) 01010 { 01011 ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo); 01012 typename ArrayView<LocalOrdinal>::iterator beg, end, newend; 01013 beg = inds_view.begin(); 01014 end = inds_view.begin() + rowinfo.numEntries; 01015 newend = beg; 01016 if (beg != end) { 01017 typename ArrayView<LocalOrdinal>::iterator cur = beg + 1; 01018 Iter vcur = rowValueIter + 1, 01019 vend = rowValueIter; 01020 cur = beg+1; 01021 while (cur != end) { 01022 if (*cur != *newend) { 01023 // new entry; save it 01024 ++newend; 01025 ++vend; 01026 (*newend) = (*cur); 01027 (*vend) = (*vcur); 01028 } 01029 else { 01030 // old entry; merge it 01031 (*vend) = f(*vend,*vcur); 01032 } 01033 ++cur; 01034 ++vcur; 01035 } 01036 ++newend; // point one past the last entry, per typical [beg,end) semantics 01037 } 01038 const size_t mergedEntries = newend - beg; 01039 #ifdef HAVE_TPETRA_DEBUG 01040 // merge should not have eliminated any entries; if so, the assignment below will destory the packed structure 01041 TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries ); 01042 #endif 01043 if (getProfileType() == StaticProfile) { 01044 rowEnds_[rowinfo.localRow] = rowBegs_[rowinfo.localRow] + mergedEntries; 01045 } 01046 else { 01047 numEntriesPerRow_[rowinfo.localRow] = mergedEntries; 01048 } 01049 nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries); 01050 } 01051 01052 01055 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01056 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setDomainRangeMaps( 01057 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 01058 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap) 01059 { 01060 // simple pointer comparison for equality 01061 if (domainMap_ != domainMap) { 01062 domainMap_ = domainMap; 01063 importer_ = null; 01064 } 01065 if (rangeMap_ != rangeMap) { 01066 rangeMap_ = rangeMap; 01067 exporter_ = null; 01068 } 01069 } 01070 01071 01074 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01075 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::findLocalIndex(RowInfo rowinfo, LocalOrdinal ind) const 01076 { 01077 typedef typename ArrayView<const LocalOrdinal>::iterator IT; 01078 bool found = true; 01079 // get a view of the row, if it wasn't passed by the caller 01080 ArrayView<const LocalOrdinal> rowinds = getLocalView(rowinfo); 01081 IT rptr, locptr = Teuchos::NullIteratorTraits<IT>::getNull(); 01082 rptr = rowinds.begin(); 01083 if (isSorted()) { 01084 // binary search 01085 std::pair<IT,IT> p = std::equal_range(rptr,rptr+rowinfo.numEntries,ind); 01086 if (p.first == p.second) found = false; 01087 else locptr = p.first; 01088 } 01089 else { 01090 // direct search 01091 locptr = std::find(rptr,rptr+rowinfo.numEntries,ind); 01092 if (locptr == rptr+rowinfo.numEntries) found = false; 01093 } 01094 size_t ret = OrdinalTraits<size_t>::invalid(); 01095 if (found) { 01096 ret = (locptr - rptr); 01097 } 01098 return ret; 01099 } 01100 01101 01104 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01105 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::findGlobalIndex(RowInfo rowinfo, GlobalOrdinal ind) const 01106 { 01107 typedef typename ArrayView<const GlobalOrdinal>::iterator IT; 01108 bool found = true; 01109 // get a view of the row, if it wasn't passed by the caller 01110 ArrayView<const GlobalOrdinal> rowinds = getGlobalView(rowinfo); 01111 IT rptr, locptr; 01112 rptr = rowinds.begin(); 01113 if (isSorted()) { 01114 // binary search 01115 std::pair<IT,IT> p = std::equal_range(rptr,rptr+rowinfo.numEntries,ind); 01116 if (p.first == p.second) found = false; 01117 else locptr = p.first; 01118 } 01119 else { 01120 // direct search 01121 locptr = std::find(rptr,rptr+rowinfo.numEntries,ind); 01122 if (locptr == rptr+rowinfo.numEntries) found = false; 01123 } 01124 size_t ret = OrdinalTraits<size_t>::invalid(); 01125 if (found) { 01126 ret = (locptr - rptr); 01127 } 01128 return ret; 01129 } 01130 01131 01134 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01135 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::clearGlobalConstants() 01136 { 01137 globalNumEntries_ = OrdinalTraits<global_size_t>::invalid(); 01138 globalNumDiags_ = OrdinalTraits<global_size_t>::invalid(); 01139 globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid(); 01140 nodeNumDiags_ = OrdinalTraits< size_t>::invalid(); 01141 nodeMaxNumRowEntries_ = OrdinalTraits< size_t>::invalid(); 01142 haveGlobalConstants_ = false; 01143 } 01144 01145 01148 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01149 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkInternalState() const 01150 { 01151 #ifdef HAVE_TPETRA_DEBUG 01152 const global_size_t GSTI = OrdinalTraits<global_size_t>::invalid(); 01153 const size_t STI = OrdinalTraits<size_t>::invalid(); 01154 std::string err = typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team."; 01155 // check the internal state of this data structure 01156 // this is called by numerous state-changing methods, in a debug build, to ensure that the object 01157 // always remains in a valid state 01158 // the graph should have been allocated with a row map 01159 TEST_FOR_EXCEPTION( rowMap_ == null, std::logic_error, err ); 01160 // am either complete or active 01161 TEST_FOR_EXCEPTION( isFillActive() == isFillComplete(), std::logic_error, err ); 01162 // if the graph has been fill completed, then all maps should be present 01163 TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err ); 01164 // if storage has been optimized, then indices should have been allocated (even if trivially so) 01165 TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err ); 01166 // if storage has been optimized, then number of allocated is now the number of entries 01167 TEST_FOR_EXCEPTION( isStorageOptimized() == true && getNodeAllocationSize() != getNodeNumEntries(), std::logic_error, err ); 01168 // if graph doesn't have the global constants, then they should all be marked as invalid 01169 TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != GSTI || globalNumDiags_ != GSTI || globalMaxNumRowEntries_ != GSTI ), std::logic_error, err ); 01170 // if the graph has global cosntants, then they should be valid. 01171 TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == GSTI || globalNumDiags_ == GSTI || globalMaxNumRowEntries_ == GSTI ), std::logic_error, err ); 01172 TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ), 01173 std::logic_error, err ); 01174 // if indices are allocated, then the allocation specifications should have been released 01175 TEST_FOR_EXCEPTION( indicesAreAllocated() == true && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null), std::logic_error, err ); 01176 // if indices are not allocated, then information dictating allocation quantities should be present 01177 TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != STI || nodeNumEntries_ != 0), std::logic_error, err ); 01178 // if storage is optimized, then profile should be static 01179 TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile, std::logic_error, err ); 01180 // rowBegs_ is required to have N+1 rows; rowBegs_[N] == gblInds1D_.size()/lclInds1D_.size() 01181 TEST_FOR_EXCEPTION( isGloballyIndexed() && rowBegs_ != null && ((size_t)rowBegs_.size() != getNodeNumRows()+1 || rowBegs_[getNodeNumRows()] != (size_t)gblInds1D_.size()), std::logic_error, err ); 01182 TEST_FOR_EXCEPTION( isLocallyIndexed() && rowBegs_ != null && ((size_t)rowBegs_.size() != getNodeNumRows()+1 || rowBegs_[getNodeNumRows()] != (size_t)lclInds1D_.size()), std::logic_error, err ); 01183 // rowEnds_ is required to have N rows 01184 TEST_FOR_EXCEPTION( rowEnds_ != null && (size_t)rowEnds_.size() != getNodeNumRows(), std::logic_error, err ); 01185 // if profile is dynamic and we have allocated, then 2D allocations should be present 01186 TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && lclInds2D_ == null && gblInds2D_ == null, 01187 std::logic_error, err ); 01188 // if profile is dynamic, then numentries and 2D indices are needed and should be present 01189 TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (numEntriesPerRow_ == null || (lclInds2D_ == null && gblInds2D_ == null)), 01190 std::logic_error, err ); 01191 // if profile is dynamic, then 1D allocations should not be present 01192 TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (lclInds1D_ != null || gblInds1D_ != null), std::logic_error, err ); 01193 // if profile is dynamic, then row offsets should not be present 01194 TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (rowBegs_ != null || rowEnds_ != null), std::logic_error, err ); 01195 // if profile is static and we have allocated non-trivially, then 1D allocations should be present 01196 TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeAllocationSize() > 0 && lclInds1D_ == null && gblInds1D_ == null, 01197 std::logic_error, err ); 01198 // if profile is static and we have a non-trivial application, then row offsets should be allocated 01199 TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (rowBegs_ == null || rowEnds_ == null), 01200 std::logic_error, err ); 01201 // if profile is static, then 2D allocations should not be present 01202 TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (lclInds2D_ != null || gblInds2D_ != null), std::logic_error, err ); 01203 // if profile is static, then we have no need for numentries and it should not be present 01204 TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeNumRows() > 0 && numEntriesPerRow_ != null, std::logic_error, err ); 01205 // if indices are not allocated, then none of the buffers should be. 01206 TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (rowBegs_ != null || rowEnds_ != null || numEntriesPerRow_ != null || 01207 lclInds1D_ != null || lclInds2D_ != null || 01208 gblInds1D_ != null || gblInds2D_ != null), std::logic_error, err ); 01209 // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_) 01210 TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false, std::logic_error, err ); 01211 // indices may be local or global, but not both 01212 TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true, std::logic_error, err ); 01213 // if indices are local, then global allocations should not be present 01214 TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (gblInds1D_ != null || gblInds2D_ != null), std::logic_error, err ); 01215 // if indices are global, then local allocations should not be present 01216 TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (lclInds1D_ != null || lclInds2D_ != null), std::logic_error, err ); 01217 // if indices are local, then local allocations should be present 01218 TEST_FOR_EXCEPTION( indicesAreLocal_ == true && getNodeAllocationSize() > 0 && lclInds1D_ == null && getNodeNumRows() > 0 && lclInds2D_ == null, 01219 std::logic_error, err ); 01220 // if indices are global, then global allocations should be present 01221 TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && getNodeAllocationSize() > 0 && gblInds1D_ == null && getNodeNumRows() > 0 && gblInds2D_ == null, 01222 std::logic_error, err ); 01223 // if indices are allocated, then we should have recorded how many were allocated 01224 TEST_FOR_EXCEPTION( indicesAreAllocated() == true && nodeNumAllocated_ == STI, std::logic_error, err ); 01225 // if indices are not allocated, then the allocation size should be marked invalid 01226 TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != STI, std::logic_error, err ); 01227 // check the actual allocations 01228 if (indicesAreAllocated()) { 01229 size_t actualNumAllocated = 0; 01230 if (pftype_ == DynamicProfile) { 01231 if (isGloballyIndexed() && gblInds2D_ != null) { 01232 for (size_t r = 0; r < getNodeNumRows(); ++r) { 01233 actualNumAllocated += gblInds2D_[r].size(); 01234 } 01235 } 01236 else if (isLocallyIndexed() && lclInds2D_ != null) { 01237 for (size_t r = 0; r < getNodeNumRows(); ++r) { 01238 actualNumAllocated += lclInds2D_[r].size(); 01239 } 01240 } 01241 } 01242 else { // pftype_ == StaticProfile) 01243 actualNumAllocated = rowBegs_[getNodeNumRows()]; 01244 TEST_FOR_EXCEPTION( isLocallyIndexed() == true && (size_t)lclInds1D_.size() != actualNumAllocated, std::logic_error, err ); 01245 TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)gblInds1D_.size() != actualNumAllocated, std::logic_error, err ); 01246 } 01247 TEST_FOR_EXCEPTION(indicesAreAllocated() == true && actualNumAllocated != nodeNumAllocated_, std::logic_error, err ); 01248 } 01249 #endif 01250 } 01251 01252 01255 // // 01256 // User-visible class methods // 01257 // // 01260 01261 01264 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01265 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const 01266 { 01267 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01268 size_t ret = OrdinalTraits<size_t>::invalid(); 01269 if (lrow != OrdinalTraits<LocalOrdinal>::invalid()) 01270 { 01271 RowInfo rowinfo = getRowInfo(lrow); 01272 ret = rowinfo.numEntries; 01273 } 01274 return ret; 01275 } 01276 01277 01280 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01281 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInLocalRow(LocalOrdinal localRow) const 01282 { 01283 size_t ret = OrdinalTraits<size_t>::invalid(); 01284 if (rowMap_->isNodeLocalElement(localRow)) { 01285 RowInfo rowinfo = getRowInfo(localRow); 01286 ret = rowinfo.numEntries; 01287 } 01288 return ret; 01289 } 01290 01291 01294 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01295 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const 01296 { 01297 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01298 size_t ret = OrdinalTraits<size_t>::invalid(); 01299 if (lrow != OrdinalTraits<LocalOrdinal>::invalid()) 01300 { 01301 RowInfo rowinfo = getRowInfo(lrow); 01302 ret = rowinfo.allocSize; 01303 } 01304 return ret; 01305 } 01306 01307 01310 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01311 size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const 01312 { 01313 size_t ret = OrdinalTraits<size_t>::invalid(); 01314 if (rowMap_->isNodeLocalElement(localRow)) { 01315 RowInfo rowinfo = getRowInfo(localRow); 01316 ret = rowinfo.allocSize; 01317 } 01318 return ret; 01319 } 01320 01321 01324 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01325 ArrayRCP<const size_t> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeRowBegs() const 01326 { 01327 return rowBegs_; 01328 } 01329 01330 01333 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01334 ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodePackedIndices() const 01335 { 01336 return lclInds1D_; 01337 } 01338 01339 01342 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01343 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowCopy(LocalOrdinal localRow, const ArrayView<LocalOrdinal> &indices, size_t &NumIndices) const 01344 { 01345 // can only do this if 01346 // * we have local indices: isLocallyIndexed() 01347 // or 01348 // * we are capable of producing them: isGloballyIndexed() && hasColMap() 01349 // short circuit if we aren't allocated 01350 std::string tfecfFuncName("getLocalRowCopy(localRow,...)"); 01351 TEST_FOR_EXCEPTION_CLASS_FUNC(isGloballyIndexed() == true && hasColMap() == false, std::runtime_error, ": local indices cannot be produced."); 01352 TEST_FOR_EXCEPTION_CLASS_FUNC(rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, 01353 ": localRow (== " << localRow << ") is not valid on this node."); 01354 const RowInfo rowinfo = getRowInfo(localRow); 01355 NumIndices = rowinfo.numEntries; 01356 TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)indices.size() < NumIndices, std::runtime_error, 01357 ": specified storage (size==" << indices.size() 01358 << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ")."); 01359 if (isLocallyIndexed()) { 01360 ArrayView<const LocalOrdinal> lview = getLocalView(rowinfo); 01361 std::copy( lview.begin(), lview.begin() + NumIndices, indices.begin()); 01362 } 01363 else if (isGloballyIndexed()) { 01364 ArrayView<const GlobalOrdinal> gview = getGlobalView(rowinfo); 01365 for (size_t j=0; j < NumIndices; ++j) { 01366 indices[j] = colMap_->getLocalElement(gview[j]); 01367 } 01368 } 01369 else { 01370 #ifdef HAVE_TPETRA_DEBUG 01371 // should have fallen in one of the above 01372 TEST_FOR_EXCEPTION_CLASS_FUNC( indicesAreAllocated() == true, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 01373 #endif 01374 NumIndices = 0; 01375 } 01376 return; 01377 } 01378 01379 01382 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01383 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowCopy(GlobalOrdinal globalRow, const ArrayView<GlobalOrdinal> &indices, size_t &NumIndices) const 01384 { 01385 // we either currently store global indices, or we have a column map with which to transcribe our local indices for the user 01386 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01387 std::string tfecfFuncName("getGlobalRowCopy(globalRow,...)"); 01388 TEST_FOR_EXCEPTION_CLASS_FUNC(lrow == OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error, 01389 ": globalRow (== " << globalRow << ") does not belong to this node."); 01390 const RowInfo rowinfo = getRowInfo((size_t)lrow); 01391 NumIndices = rowinfo.numEntries; 01392 TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)indices.size() < NumIndices, std::runtime_error, 01393 ": specified storage (size==" << indices.size() 01394 << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ")."); 01395 if (isLocallyIndexed()) { 01396 ArrayView<const LocalOrdinal> lview = getLocalView(rowinfo); 01397 for (size_t j=0; j < NumIndices; ++j) { 01398 indices[j] = colMap_->getGlobalElement(lview[j]); 01399 } 01400 } 01401 else if (isGloballyIndexed()) { 01402 ArrayView<const GlobalOrdinal> gview = getGlobalView(rowinfo); 01403 std::copy(gview.begin(), gview.begin() + NumIndices, indices.begin()); 01404 } 01405 return; 01406 } 01407 01408 01411 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01412 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowView(LocalOrdinal localRow, ArrayView<const LocalOrdinal> &indices) const 01413 { 01414 std::string tfecfFuncName("getLocalRowView()"); 01415 TEST_FOR_EXCEPTION_CLASS_FUNC(isGloballyIndexed() == true, std::runtime_error, ": local indices cannot be provided."); 01416 indices = null; 01417 if (rowMap_->isNodeLocalElement(localRow) == true) { 01418 const RowInfo rowinfo = getRowInfo(localRow); 01419 if (rowinfo.numEntries > 0) { 01420 indices = getLocalView(rowinfo); 01421 indices = indices(0,rowinfo.numEntries); 01422 } 01423 } 01424 #ifdef HAVE_TPETRA_DEBUG 01425 TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)indices.size() != getNumEntriesInLocalRow(localRow), std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team."); 01426 #endif 01427 return; 01428 } 01429 01430 01433 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01434 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowView(GlobalOrdinal globalRow, ArrayView<const GlobalOrdinal> &indices) const 01435 { 01436 std::string tfecfFuncName("getGlobalRowView()"); 01437 TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": global indices cannot be provided."); 01438 indices = null; 01439 if (rowMap_->isNodeGlobalElement(globalRow) == true) { 01440 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 01441 const RowInfo rowinfo = getRowInfo(lrow); 01442 if (rowinfo.numEntries > 0) { 01443 indices = getGlobalView(rowinfo); 01444 indices = indices(0,rowinfo.numEntries); 01445 } 01446 } 01447 #ifdef HAVE_TPETRA_DEBUG 01448 TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)indices.size() != getNumEntriesInGlobalRow(globalRow), std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team."); 01449 #endif 01450 return; 01451 } 01452 01453 01456 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01457 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertLocalIndices( 01458 LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &indices) 01459 { 01460 std::string tfecfFuncName("insertLocalIndices()"); 01461 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active."); 01462 TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true, std::runtime_error, ": graph indices are global; use insertGlobalIndices()."); 01463 TEST_FOR_EXCEPTION_CLASS_FUNC( hasColMap() == false, std::runtime_error, ": cannot insert local indices without a column map."); 01464 TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, ": row does not belong to this node."); 01465 if (indicesAreAllocated() == false) { 01466 allocateIndices(LocalIndices); 01467 } 01468 // use column map to filter the entries 01469 Array<LocalOrdinal> f_inds(indices); 01470 SLocalGlobalNCViews inds_ncview; 01471 inds_ncview.linds = f_inds(); 01472 const size_t numFilteredEntries = filterIndices<LocalIndices>(inds_ncview); 01473 if (numFilteredEntries > 0) { 01474 RowInfo rowInfo = getRowInfo(localRow); 01475 const size_t curNumEntries = rowInfo.numEntries; 01476 const size_t newNumEntries = curNumEntries + numFilteredEntries; 01477 if (newNumEntries > rowInfo.allocSize) { 01478 TEST_FOR_EXCEPTION_CLASS_FUNC(getProfileType() == StaticProfile, std::runtime_error, ": new indices exceed statically allocated graph structure."); 01479 TPETRA_EFFICIENCY_WARNING(true,std::runtime_error, 01480 "::insertLocalIndices(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation."); 01481 // update allocation only as much as necessary 01482 rowInfo = updateAlloc<LocalIndices>(rowInfo, newNumEntries); 01483 } 01484 SLocalGlobalViews inds_view; 01485 inds_view.linds = f_inds(0,numFilteredEntries); 01486 insertIndices<LocalIndices>(rowInfo, inds_view); 01487 } 01488 #ifdef HAVE_TPETRA_DEBUG 01489 TEST_FOR_EXCEPTION_CLASS_FUNC(indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error, 01490 ": Violated stated post-conditions. Please contact Tpetra team."); 01491 #endif 01492 } 01493 01494 01497 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01498 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertGlobalIndices( 01499 GlobalOrdinal grow, const ArrayView<const GlobalOrdinal> &indices) 01500 { 01501 std::string tfecfFuncName("insertGlobalIndices()"); 01502 TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": graph indices are local; use insertLocalIndices()."); 01503 // this can't really be satisfied for now, because if we are fillComplete(), then we are local 01504 // in the future, this may change. however, the rule that modification require active fill will not change. 01505 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active."); 01506 if (indicesAreAllocated() == false) { 01507 allocateIndices(GlobalIndices); 01508 } 01509 const LocalOrdinal myRow = rowMap_->getLocalElement(grow); 01510 if (myRow != OrdinalTraits<LocalOrdinal>::invalid()) 01511 { 01512 // if we have a column map, use it to filter the entries. 01513 Array<GlobalOrdinal> filtered_indices; 01514 SLocalGlobalViews inds_view; 01515 if (hasColMap()) { 01516 SLocalGlobalNCViews inds_ncview; 01517 // filter indices and values through the column map 01518 filtered_indices.assign(indices.begin(), indices.end()); 01519 inds_ncview.ginds = filtered_indices(); 01520 const size_t numFilteredEntries = filterIndices<GlobalIndices>(inds_ncview); 01521 inds_view.ginds = filtered_indices(0,numFilteredEntries); 01522 } 01523 else { 01524 inds_view.ginds = indices; 01525 } 01526 const size_t numFilteredEntries = inds_view.ginds.size(); 01527 // add the new indices and values 01528 if (numFilteredEntries > 0) { 01529 RowInfo rowInfo = getRowInfo(myRow); 01530 const size_t curNumEntries = rowInfo.numEntries; 01531 const size_t newNumEntries = curNumEntries + numFilteredEntries; 01532 if (newNumEntries > rowInfo.allocSize) { 01533 TEST_FOR_EXCEPTION_CLASS_FUNC(getProfileType() == StaticProfile, std::runtime_error, ": new indices exceed statically allocated graph structure."); 01534 TPETRA_EFFICIENCY_WARNING(true,std::runtime_error, 01535 "::insertGlobalValues(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation."); 01536 // update allocation only as much as necessary 01537 rowInfo = updateAlloc<GlobalIndices>(rowInfo, newNumEntries); 01538 } 01539 insertIndices<GlobalIndices>(rowInfo, inds_view); 01540 #ifdef HAVE_TPETRA_DEBUG 01541 { 01542 const size_t chkNewNumEntries = getNumEntriesInLocalRow(myRow); 01543 TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 01544 } 01545 #endif 01546 } 01547 } 01548 else { 01549 // a nonlocal row 01550 for (typename ArrayView<const GlobalOrdinal>::iterator i=indices.begin(); i != indices.end(); ++i) { 01551 nonlocals_[grow].push_back(*i); 01552 } 01553 } 01554 #ifdef HAVE_TPETRA_DEBUG 01555 TEST_FOR_EXCEPTION_CLASS_FUNC(indicesAreAllocated() == false || isGloballyIndexed() == false, std::logic_error, 01556 ": Violated stated post-conditions. Please contact Tpetra team."); 01557 #endif 01558 } 01559 01560 01563 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01564 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::removeLocalIndices(LocalOrdinal lrow) 01565 { 01566 std::string tfecfFuncName("removeLocalIndices()"); 01567 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active."); 01568 TEST_FOR_EXCEPTION_CLASS_FUNC( isStorageOptimized() == true, std::runtime_error, ": cannot remove indices after optimizeStorage() has been called."); 01569 TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true, std::runtime_error, ": graph indices are global; use removeGlobalIndices()."); 01570 TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error, ": row does not belong to this node."); 01571 if (indicesAreAllocated() == false) { 01572 allocateIndices(LocalIndices); 01573 } 01574 // 01575 clearGlobalConstants(); 01576 // 01577 RowInfo sizeInfo = getRowInfo(lrow); 01578 if (sizeInfo.allocSize > 0 && numEntriesPerRow_ != null) { 01579 numEntriesPerRow_[lrow] = 0; 01580 } 01581 #ifdef HAVE_TPETRA_DEBUG 01582 TEST_FOR_EXCEPTION_CLASS_FUNC(getNumEntriesInLocalRow(lrow) != 0 || indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error, 01583 ": Violated stated post-conditions. Please contact Tpetra team."); 01584 #endif 01585 } 01586 01587 01588 // TODO: in the future, globalAssemble() should use import/export functionality 01591 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01592 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::globalAssemble() 01593 { 01594 using std::deque; 01595 using std::pair; 01596 using std::make_pair; 01597 typedef typename std::map<GlobalOrdinal,std::deque<GlobalOrdinal> >::const_iterator NLITER; 01598 int numImages = Teuchos::size(*getComm()); 01599 int myImageID = Teuchos::rank(*getComm()); 01600 #ifdef HAVE_TPETRA_DEBUG 01601 Teuchos::barrier( *rowMap_->getComm() ); 01602 #endif 01603 std::string tfecfFuncName("globalAssemble()"); 01604 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active."); 01605 // Determine if any nodes have global entries to share 01606 { 01607 size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals; 01608 Teuchos::reduceAll<int,size_t>(*getComm(),Teuchos::REDUCE_MAX,MyNonlocals, 01609 outArg(MaxGlobalNonlocals)); 01610 if (MaxGlobalNonlocals == 0) return; // no entries to share 01611 } 01612 01613 // compute a list of NLRs from nonlocals_ and use it to compute: 01614 // IdsAndRows: a vector of (id,row) pairs 01615 // NLR2Id: a map from NLR to the Id that owns it 01616 // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i 01617 // sendIDs: a list of all images I send to 01618 // recvIDs: a list of all images I receive from (constructed later) 01619 Array<pair<int,GlobalOrdinal> > IdsAndRows; 01620 std::map<GlobalOrdinal,int> NLR2Id; 01621 Teuchos::SerialDenseMatrix<int,char> globalNeighbors; 01622 Array<int> sendIDs, recvIDs; 01623 { 01624 // nonlocals_ contains the entries we are holding for all non-local rows 01625 // we want a list of the rows for which we have data 01626 Array<GlobalOrdinal> NLRs; 01627 std::set<GlobalOrdinal> setOfRows; 01628 for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) { 01629 setOfRows.insert(iter->first); 01630 } 01631 // copy the elements in the set into an Array 01632 NLRs.resize(setOfRows.size()); 01633 std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin()); 01634 01635 // get a list of ImageIDs for the non-local rows (NLRs) 01636 Array<int> NLRIds(NLRs.size()); 01637 { 01638 LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds()); 01639 char lclerror = ( stat == IDNotPresent ? 1 : 0 ); 01640 char gblerror; 01641 Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,lclerror,outArg(gblerror)); 01642 std::string tfecfFuncName("globalAssemble()"); 01643 TEST_FOR_EXCEPTION_CLASS_FUNC(gblerror != 0, std::runtime_error, ": non-local entries correspond to invalid rows."); 01644 } 01645 01646 // build up a list of neighbors, as well as a map between NLRs and Ids 01647 // localNeighbors[i] != 0 iff I have data to send to image i 01648 // put NLRs,Ids into an array of pairs 01649 IdsAndRows.reserve(NLRs.size()); 01650 Array<char> localNeighbors(numImages,0); 01651 typename Array<GlobalOrdinal>::const_iterator nlr; 01652 typename Array<int>::const_iterator id; 01653 for (nlr = NLRs.begin(), id = NLRIds.begin(); 01654 nlr != NLRs.end(); ++nlr, ++id) { 01655 NLR2Id[*nlr] = *id; 01656 localNeighbors[*id] = 1; 01657 IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr)); 01658 } 01659 for (int j=0; j<numImages; ++j) { 01660 if (localNeighbors[j]) { 01661 sendIDs.push_back(j); 01662 } 01663 } 01664 // sort IdsAndRows, by Ids first, then rows 01665 std::sort(IdsAndRows.begin(),IdsAndRows.end()); 01666 // gather from other nodes to form the full graph 01667 globalNeighbors.shapeUninitialized(numImages,numImages); 01668 Teuchos::gatherAll(*getComm(),numImages,localNeighbors.getRawPtr(),numImages*numImages,globalNeighbors.values()); 01669 // globalNeighbors at this point contains (on all images) the 01670 // connectivity between the images. 01671 // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j 01672 } 01673 01675 // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH 01676 // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID 01678 01679 // loop over all columns to know from which images I can expect to receive something 01680 for (int j=0; j<numImages; ++j) { 01681 if (globalNeighbors(myImageID,j)) { 01682 recvIDs.push_back(j); 01683 } 01684 } 01685 const size_t numRecvs = recvIDs.size(); 01686 01687 // we know how many we're sending to already 01688 // form a contiguous list of all data to be sent 01689 // track the number of entries for each ID 01690 Array<pair<GlobalOrdinal,GlobalOrdinal> > IJSendBuffer; 01691 Array<size_t> sendSizes(sendIDs.size(), 0); 01692 size_t numSends = 0; 01693 for (typename Array<pair<int,GlobalOrdinal> >::const_iterator IdAndRow = IdsAndRows.begin(); 01694 IdAndRow != IdsAndRows.end(); ++IdAndRow) { 01695 int id = IdAndRow->first; 01696 GlobalOrdinal row = IdAndRow->second; 01697 // have we advanced to a new send? 01698 if (sendIDs[numSends] != id) { 01699 numSends++; 01700 TEST_FOR_EXCEPTION_CLASS_FUNC(sendIDs[numSends] != id, std::logic_error, ": internal logic error. Contact Tpetra team."); 01701 } 01702 // copy data for row into contiguous storage 01703 for (typename deque<GlobalOrdinal>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j) 01704 { 01705 IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) ); 01706 sendSizes[numSends]++; 01707 } 01708 } 01709 if (IdsAndRows.size() > 0) { 01710 numSends++; // one last increment, to make it a count instead of an index 01711 } 01712 TEST_FOR_EXCEPTION_CLASS_FUNC(Teuchos::as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, ": internal logic error. Contact Tpetra team."); 01713 01714 // don't need this data anymore 01715 nonlocals_.clear(); 01716 01718 // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS 01720 // perform non-blocking sends: send sizes to our recipients 01721 Array<RCP<Teuchos::CommRequest> > sendRequests; 01722 for (size_t s=0; s < numSends ; ++s) { 01723 // we'll fake the memory management, because all communication will be local to this method and the scope of our data 01724 sendRequests.push_back( Teuchos::isend<int,size_t>(*getComm(),rcp<size_t>(&sendSizes[s],false),sendIDs[s]) ); 01725 } 01726 // perform non-blocking receives: receive sizes from our senders 01727 Array<RCP<Teuchos::CommRequest> > recvRequests; 01728 Array<size_t> recvSizes(numRecvs); 01729 for (size_t r=0; r < numRecvs; ++r) { 01730 // we'll fake the memory management, because all communication will be local to this method and the scope of our data 01731 recvRequests.push_back( Teuchos::ireceive(*getComm(),rcp(&recvSizes[r],false),recvIDs[r]) ); 01732 } 01733 // wait on all 01734 if (!sendRequests.empty()) { 01735 Teuchos::waitAll(*getComm(),sendRequests()); 01736 } 01737 if (!recvRequests.empty()) { 01738 Teuchos::waitAll(*getComm(),recvRequests()); 01739 } 01740 Teuchos::barrier(*getComm()); 01741 sendRequests.clear(); 01742 recvRequests.clear(); 01743 01745 // NOW SEND/RECEIVE ALL ROW DATA 01747 // from the size info, build the ArrayViews into IJSendBuffer 01748 Array<ArrayView<pair<GlobalOrdinal,GlobalOrdinal> > > sendBuffers(numSends,null); 01749 { 01750 size_t cur = 0; 01751 for (size_t s=0; s<numSends; ++s) { 01752 sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]); 01753 cur += sendSizes[s]; 01754 } 01755 } 01756 // perform non-blocking sends 01757 for (size_t s=0; s < numSends ; ++s) 01758 { 01759 // we'll fake the memory management, because all communication will be local to this method and the scope of our data 01760 ArrayRCP<pair<GlobalOrdinal,GlobalOrdinal> > tmparcp = arcp(sendBuffers[s].getRawPtr(),0,sendBuffers[s].size(),false); 01761 sendRequests.push_back( Teuchos::isend<int,pair<GlobalOrdinal,GlobalOrdinal> >(*getComm(),tmparcp,sendIDs[s]) ); 01762 } 01763 // calculate amount of storage needed for receives 01764 // setup pointers for the receives as well 01765 size_t totalRecvSize = std::accumulate(recvSizes.begin(),recvSizes.end(),0); 01766 Array<pair<GlobalOrdinal,GlobalOrdinal> > IJRecvBuffer(totalRecvSize); 01767 // from the size info, build the ArrayViews into IJRecvBuffer 01768 Array<ArrayView<pair<GlobalOrdinal,GlobalOrdinal> > > recvBuffers(numRecvs,null); 01769 { 01770 size_t cur = 0; 01771 for (size_t r=0; r<numRecvs; ++r) { 01772 recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]); 01773 cur += recvSizes[r]; 01774 } 01775 } 01776 // perform non-blocking recvs 01777 for (size_t r=0; r < numRecvs ; ++r) { 01778 // we'll fake the memory management, because all communication will be local to this method and the scope of our data 01779 ArrayRCP<pair<GlobalOrdinal,GlobalOrdinal> > tmparcp = arcp(recvBuffers[r].getRawPtr(),0,recvBuffers[r].size(),false); 01780 recvRequests.push_back( Teuchos::ireceive(*getComm(),tmparcp,recvIDs[r]) ); 01781 } 01782 // perform waits 01783 if (!sendRequests.empty()) { 01784 Teuchos::waitAll(*getComm(),sendRequests()); 01785 } 01786 if (!recvRequests.empty()) { 01787 Teuchos::waitAll(*getComm(),recvRequests()); 01788 } 01789 Teuchos::barrier(*getComm()); 01790 sendRequests.clear(); 01791 recvRequests.clear(); 01792 01794 // NOW PROCESS THE RECEIVED ROW DATA 01796 // TODO: instead of adding one entry at a time, add one row at a time. 01797 // this requires resorting; they arrived sorted by sending node, so that entries could be non-contiguous if we received 01798 // multiple entries for a particular row from different processors. 01799 // it also requires restoring the data, which may make it not worth the trouble. 01800 for (typename Array<pair<GlobalOrdinal,GlobalOrdinal> >::const_iterator ij = IJRecvBuffer.begin(); ij != IJRecvBuffer.end(); ++ij) 01801 { 01802 insertGlobalIndices(ij->first, tuple<GlobalOrdinal>(ij->second)); 01803 } 01804 checkInternalState(); 01805 } 01806 01807 01810 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01811 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::resumeFill() 01812 { 01813 #ifdef HAVE_TPETRA_DEBUG 01814 Teuchos::barrier( *rowMap_->getComm() ); 01815 #endif 01816 clearGlobalConstants(); 01817 lclGraph_.clear(); 01818 lowerTriangular_ = false; 01819 upperTriangular_ = false; 01820 indicesAreSorted_ = false; 01821 noRedundancies_ = false; 01822 fillComplete_ = false; 01823 #ifdef HAVE_TPETRA_DEBUG 01824 std::string tfecfFuncName("resumeFill()"); 01825 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false || isFillComplete() == true, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team."); 01826 #endif 01827 } 01828 01829 01832 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01833 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(OptimizeOption os) 01834 { 01835 fillComplete(rowMap_,rowMap_,os); 01836 } 01837 01838 01841 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01842 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete( 01843 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 01844 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, 01845 OptimizeOption os) 01846 { 01847 #ifdef HAVE_TPETRA_DEBUG 01848 Teuchos::barrier( *rowMap_->getComm() ); 01849 #endif 01850 std::string tfecfFuncName("fillComplete()"); 01851 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false || isFillComplete() == true, std::runtime_error, ": Graph fill state must be active."); 01852 // allocate if unallocated 01853 if (indicesAreAllocated() == false) { 01854 // allocate global, in case we do not have a column map 01855 allocateIndices( GlobalIndices ); 01856 } 01857 // global assemble 01858 if (Teuchos::size(*getComm()) > 1) { 01859 globalAssemble(); 01860 } 01861 else { 01862 TEST_FOR_EXCEPTION_CLASS_FUNC(nonlocals_.size() > 0, std::runtime_error, ": cannot have non-local entries on a serial run. Invalid entries were submitted to the CrsMatrix."); 01863 } 01864 // set domain/range map: may clear the import/export objects 01865 setDomainRangeMaps(domainMap,rangeMap); 01866 // make column map 01867 if (hasColMap() == false) { 01868 makeColMap(); 01869 } 01870 // make indices local 01871 if (isGloballyIndexed() == true) { 01872 makeIndicesLocal(); 01873 } 01874 // sort entries 01875 if (isSorted() == false) { 01876 sortAllIndices(); 01877 } 01878 // merge entries 01879 if (isMerged() == false) { 01880 mergeAllIndices(); 01881 } 01882 // make import/export objects 01883 makeImportExport(); 01884 // compute global constants 01885 computeGlobalConstants(); 01886 // fill local objects 01887 fillLocalGraph(os); 01888 // 01889 fillComplete_ = true; 01890 #ifdef HAVE_TPETRA_DEBUG 01891 TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team."); 01892 #endif 01893 // 01894 checkInternalState(); 01895 } 01896 01897 01900 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01901 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::pushToLocalGraph() 01902 { 01903 std::string tfecfFuncName("pushToLocalGraph()"); 01904 TEST_FOR_EXCEPTION_CLASS_FUNC(lclGraph_.isFinalized() == true, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 01905 // fill local graph 01906 if (getProfileType() == StaticProfile) { 01907 lclGraph_.set1DStructure(lclInds1D_,rowBegs_,rowEnds_); 01908 lclInds1D_ = null; 01909 rowBegs_ = null; 01910 rowEnds_ = null; 01911 } 01912 else { 01913 lclGraph_.set2DStructure(lclInds2D_,numEntriesPerRow_); 01914 lclInds2D_ = null; 01915 numEntriesPerRow_ = null; 01916 } 01917 nodeNumAllocated_ = 0; 01918 nodeNumEntries_ = 0; 01919 } 01920 01921 01924 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01925 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::pullFromLocalGraph() 01926 { 01927 std::string tfecfFuncName("pullFromLocalGraph()"); 01928 TEST_FOR_EXCEPTION_CLASS_FUNC(lclGraph_.isFinalized() == false, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 01929 // get new data from local graph 01930 // this requires updating the allocation size, but not necessarily the 01931 if (lclGraph_.is1DStructure()) { 01932 lclGraph_.get1DStructure(lclInds1D_,rowBegs_,rowEnds_); 01933 pftype_ = StaticProfile; 01934 #ifdef HAVE_TPETRA_DEBUG 01935 TEST_FOR_EXCEPT( rowBegs_ == null ); 01936 #endif 01937 nodeNumAllocated_ = rowBegs_[getNodeNumRows()]; 01938 if (nodeNumAllocated_ > 0) { 01939 lclInds1D_ = lclInds1D_.persistingView(0,nodeNumAllocated_); 01940 } 01941 nodeNumEntries_ = 0; 01942 for (size_t r=0; r < getNodeNumRows(); ++r) { 01943 nodeNumEntries_ += rowEnds_[r] - rowBegs_[r]; 01944 } 01945 } 01946 else { 01947 lclGraph_.get2DStructure(lclInds2D_,numEntriesPerRow_); 01948 pftype_ = DynamicProfile; 01949 nodeNumAllocated_ = 0; 01950 nodeNumEntries_ = 0; 01951 for (size_t r=0; r < getNodeNumRows(); ++r) { 01952 nodeNumAllocated_ += lclInds2D_[r].size(); 01953 nodeNumEntries_ += numEntriesPerRow_[r]; 01954 } 01955 } 01956 } 01957 01958 01961 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01962 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalGraph(OptimizeOption os) 01963 { 01964 pushToLocalGraph(); 01965 // finalize local graph(with os) 01966 const bool optStorage = (os == DoOptimizeStorage); 01967 lclGraph_.finalize( optStorage ); 01968 // get the data back from the local objects 01969 pullFromLocalGraph(); 01970 } 01971 01972 01975 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01976 const Kokkos::CrsGraph<LocalOrdinal,Node,LocalMatOps> & 01977 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraph() const 01978 { 01979 return lclGraph_; 01980 } 01981 01982 01985 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01986 Kokkos::CrsGraph<LocalOrdinal,Node,LocalMatOps> & 01987 CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraphNonConst() 01988 { 01989 return lclGraph_; 01990 } 01991 01992 01995 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 01996 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeGlobalConstants() 01997 { 01998 // compute the local constants first 01999 const size_t nlrs = getNodeNumRows(); 02000 // reset all local properties 02001 upperTriangular_ = true; 02002 lowerTriangular_ = true; 02003 nodeMaxNumRowEntries_ = 0; 02004 nodeNumDiags_ = 0; 02005 // indices are already sorted in each row 02006 const Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap = *rowMap_; 02007 if (indicesAreAllocated() == true && nodeNumAllocated_ > 0) { 02008 for (size_t r=0; r < nlrs; ++r) { 02009 GlobalOrdinal rgid = rowMap.getGlobalElement(r); 02010 // determine the local column index for this row, used for delimiting the diagonal 02011 const LocalOrdinal rlcid = colMap_->getLocalElement(rgid); 02012 RowInfo rowinfo = getRowInfo(r); 02013 ArrayView<const LocalOrdinal> rview = getLocalView(rowinfo); 02014 typename ArrayRCP<const LocalOrdinal>::iterator beg, end, cur; 02015 beg = rview.begin(); 02016 end = beg + rowinfo.numEntries; 02017 if (beg != end) { 02018 for (cur = beg; cur != end; ++cur) { 02019 // is this the diagonal? 02020 if (rlcid == (*cur)) ++nodeNumDiags_; 02021 } 02022 // because of sorting, smallest column index is (*beg); it indicates upper triangularity 02023 if (Teuchos::as<size_t>(beg[0]) < r) upperTriangular_ = false; 02024 // because of sorting, largest column index is (*newend); it indicates lower triangularity 02025 if (r < Teuchos::as<size_t>(end[-1])) lowerTriangular_ = false; 02026 } 02027 // compute num entries for this row, accumulate into nodeNumEntries_, update nodeMaxNumRowEntries_ 02028 nodeMaxNumRowEntries_ = std::max( nodeMaxNumRowEntries_, rowinfo.numEntries ); 02029 } 02030 } 02031 02032 // compute global constants using computed local constants 02033 if (haveGlobalConstants_ == false) { 02034 global_size_t lcl[2], gbl[2]; 02035 lcl[0] = nodeNumEntries_; 02036 lcl[1] = nodeNumDiags_; 02037 Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_SUM,2,lcl,gbl); 02038 globalNumEntries_ = gbl[0]; 02039 globalNumDiags_ = gbl[1]; 02040 Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_MAX,nodeMaxNumRowEntries_, 02041 outArg(globalMaxNumRowEntries_)); 02042 haveGlobalConstants_ = true; 02043 } 02044 } 02045 02046 02049 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02050 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeIndicesLocal() 02051 { 02052 // All nodes must be in the same index state. 02053 // Update index state by checking isLocallyIndexed/Global on all nodes 02054 computeIndexState(); 02055 std::string tfecfFuncName("makeIndicesLocal()"); 02056 TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() && isGloballyIndexed(), std::logic_error, 02057 ": inconsistent index state. Indices must be local on all nodes or global on all nodes."); 02058 // If user has not prescribed column map, create one from indices 02059 makeColMap(); 02060 // Transform indices to local index space 02061 const size_t nlrs = getNodeNumRows(); 02062 // 02063 if (isGloballyIndexed() && nlrs > 0) { 02064 // allocate data for local indices 02065 if (getProfileType() == StaticProfile) { 02066 // reinterpret the the compute buffer as LocalOrdinal 02067 if (nodeNumAllocated_) { 02068 lclInds1D_ = arcp_reinterpret_cast<LocalOrdinal>(gblInds1D_).persistingView(0,nodeNumAllocated_); 02069 } 02070 for (size_t r=0; r < getNodeNumRows(); ++r) { 02071 const size_t offset = rowBegs_[r], 02072 numentry = rowEnds_[r] - rowBegs_[r]; 02073 for (size_t j=0; j<numentry; ++j) { 02074 GlobalOrdinal gid = gblInds1D_[offset + j]; 02075 LocalOrdinal lid = colMap_->getLocalElement(gid); 02076 lclInds1D_[offset + j] = lid; 02077 #ifdef HAVE_TPETRA_DEBUG 02078 TEST_FOR_EXCEPTION_CLASS_FUNC(lclInds1D_[offset + j] == OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error, 02079 ": Internal error in fillComplete(). Please contact Tpetra team."); 02080 #endif 02081 } 02082 } 02083 // done with this pointer (allocation will persist in lclInds1D_) 02084 gblInds1D_ = null; 02085 } 02086 else { // getProfileType() == DynamicProfile 02087 lclInds2D_ = arcp< ArrayRCP<LocalOrdinal> >(nlrs); 02088 // if we have views, then make views 02089 for (size_t r=0; r < getNodeNumRows(); ++r) { 02090 if (gblInds2D_[r] != null) { 02091 ArrayRCP<GlobalOrdinal> ginds = gblInds2D_[r]; 02092 const size_t rna = gblInds2D_[r].size(); 02093 ArrayRCP< LocalOrdinal> linds = arcp_reinterpret_cast<LocalOrdinal>(ginds).persistingView(0,rna); 02094 // do the conversion in situ. this must be done from front to back. 02095 const size_t numentry = numEntriesPerRow_[r]; 02096 for (size_t j=0; j < numentry; ++j) { 02097 GlobalOrdinal gid = ginds[j]; 02098 LocalOrdinal lid = colMap_->getLocalElement(gid); 02099 linds[j] = lid; 02100 #ifdef HAVE_TPETRA_DEBUG 02101 TEST_FOR_EXCEPTION_CLASS_FUNC(linds[j] == OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error, 02102 ": Internal error in makeIndicesLocal(). Please contact Tpetra team."); 02103 #endif 02104 } 02105 // reinterpret the the compute buffer as LocalOrdinal; keep the original size 02106 lclInds2D_[r] = linds; 02107 gblInds2D_[r] = null; 02108 } 02109 } 02110 gblInds2D_ = null; 02111 } 02112 } 02113 indicesAreLocal_ = true; 02114 indicesAreGlobal_ = false; 02115 checkInternalState(); 02116 } 02117 02118 02121 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02122 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeIndexState() 02123 { 02124 char myIndices[2] = {0,0}; 02125 if (indicesAreLocal_) myIndices[0] = 1; 02126 if (indicesAreGlobal_) myIndices[1] = 1; 02127 char allIndices[2]; 02128 Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,2,myIndices,allIndices); 02129 indicesAreLocal_ = (allIndices[0]==1); // If indices are local on one PE, should be local on all 02130 indicesAreGlobal_ = (allIndices[1]==1); // If indices are global on one PE should be local on all 02131 } 02132 02133 02134 02137 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02138 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortAllIndices() 02139 { 02140 TEST_FOR_EXCEPT(isGloballyIndexed()==true); // this should be called only after makeIndicesLocal() 02141 if (isSorted() == false) { 02142 for (size_t row=0; row < getNodeNumRows(); ++row) { 02143 RowInfo rowInfo = getRowInfo(row); 02144 sortRowIndices(rowInfo); 02145 } 02146 // we just sorted every row 02147 setSorted(true); 02148 } 02149 } 02150 02151 02154 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02155 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeColMap() 02156 { 02157 std::string tfecfFuncName("makeColMap()"); 02158 typedef OrdinalTraits<GlobalOrdinal> GOT; 02159 // 02160 if (hasColMap()) return; 02161 const size_t nlrs = getNodeNumRows(); 02162 // 02163 computeIndexState(); 02164 TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": indices must still be global when making the column map."); 02165 // ultimate goal: list of indices for column map 02166 Array<GlobalOrdinal> myColumns; 02167 // if isGloballyIndexed() == false and isLocallyIndexed() == false, then we have nothing to do 02168 if (isGloballyIndexed() == true) { 02169 // Construct two lists of columns for which we have a non-zero 02170 // Local GIDs: Column GIDs which are locally present on the Domain map 02171 // Remote GIDs: Column GIDs which are not locally present present on the Domain map 02172 // 02173 // instead of list of local GIDs, we will use a set<LocalOrdinal> LocalGID. 02174 // 02175 const LocalOrdinal LINV = OrdinalTraits<LocalOrdinal>::invalid(); 02176 size_t numLocalColGIDs = 0, numRemoteColGIDs = 0; 02177 // 02178 // intitial: partitioning into local and remote 02179 Array<char> GIDisLocal(domainMap_->getNodeNumElements(),0); 02180 std::set<GlobalOrdinal> RemoteGIDSet; 02181 for (size_t r=0; r < nlrs; ++r) { 02182 RowInfo rowinfo = getRowInfo(r); 02183 if (rowinfo.numEntries > 0) { 02184 ArrayView<const GlobalOrdinal> rowgids = getGlobalView(rowinfo); 02185 rowgids = rowgids(0,rowinfo.numEntries); 02186 for (typename ArrayView<const GlobalOrdinal>::iterator cind = rowgids.begin(); cind != rowgids.end(); ++cind) 02187 { 02188 GlobalOrdinal gid = (*cind); 02189 LocalOrdinal lid = domainMap_->getLocalElement(gid); 02190 if (lid != LINV) { 02191 char alreadyFound = GIDisLocal[lid]; 02192 if (alreadyFound == 0) { 02193 GIDisLocal[lid] = 1; 02194 ++numLocalColGIDs; 02195 } 02196 } 02197 else { 02198 std::pair<typename std::set<GlobalOrdinal>::iterator, bool> ip; 02199 ip = RemoteGIDSet.insert(gid); 02200 if (ip.second == true) { // gid did not exist in the set and was actually inserted 02201 ++numRemoteColGIDs; 02202 } 02203 } 02204 } 02205 } 02206 } 02207 02208 // Possible short-circuit for serial scenario 02209 // If the all domain GIDs are present as column indices, then set ColMap=DomainMap 02210 // By construction, LocalGIDs \subset DomainGIDs 02211 // If we have 02212 // * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs, 02213 // and 02214 // * Number of local GIDs is number of domain GIDs 02215 // then 02216 // * LocalGIDs \subset DomainGIDs && size(LocalGIDs) == size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs 02217 // on this node. 02218 // We will concern ourselves only with the special case of a serial DomainMap, obliviating the need for a communication. 02219 // If 02220 // * DomainMap has a serial comm 02221 // then we can set Column map as Domain map and return. Benefit: this graph won't need an Import object later. 02222 // 02223 // Note, for a serial domain map, there can be no RemoteGIDs, because there are no remote nodes. 02224 // Likely explanations for this are: 02225 // * user submitted erroneous column indices 02226 // * user submitted erroneous domain map 02227 if (Teuchos::size(*domainMap_->getComm()) == 1) 02228 { 02229 TEST_FOR_EXCEPTION_CLASS_FUNC(numRemoteColGIDs != 0, std::runtime_error, 02230 ": Some column IDs are not in the domain map." << std::endl 02231 << "Either these column IDs are invalid or the domain map is invalid." << std::endl 02232 << "Remember, for a rectangular matrix, the domain map must be passed to fillComplete()."); 02233 if (numLocalColGIDs == domainMap_->getNodeNumElements()) { 02234 colMap_ = domainMap_; 02235 checkInternalState(); 02236 return; 02237 } 02238 } 02239 02240 // Now, populate myColumns() with a list of all column GIDs. 02241 // Put local GIDs at the front: they correspond to "same" and "permuted" entries between the column map and the domain map 02242 // Put remote GIDs at the back 02243 myColumns.resize(numLocalColGIDs + numRemoteColGIDs); 02244 // get pointers into myColumns for each part 02245 ArrayView<GlobalOrdinal> LocalColGIDs, RemoteColGIDs; 02246 LocalColGIDs = myColumns(0,numLocalColGIDs); 02247 RemoteColGIDs = myColumns(numLocalColGIDs,numRemoteColGIDs); 02248 02249 // Copy the Remote GIDs into myColumns 02250 std::copy(RemoteGIDSet.begin(), RemoteGIDSet.end(), RemoteColGIDs.begin()); 02251 // We will make a list of corresponding node IDs 02252 Array<int> RemoteImageIDs(numRemoteColGIDs); 02253 // Lookup the Remote Nodes IDs in the Domain map 02254 { 02255 LookupStatus stat = domainMap_->getRemoteIndexList(RemoteColGIDs, RemoteImageIDs()); 02256 // This error check is crucial: this tells us that one of the remote indices was not present in the domain map. 02257 // This means that the Import object cannot be constructed, because of incongruity between the column map and domain map. 02258 // * The user has made a mistake in the column indices 02259 // * The user has made a mistake w.r.t. the domain map 02260 // Same error message as above for serial case. 02261 char missingID_lcl = (stat == IDNotPresent ? 1 : 0); 02262 char missingID_gbl; 02263 Teuchos::reduceAll<int,char>(*getComm(),Teuchos::REDUCE_MAX,missingID_lcl, 02264 outArg(missingID_gbl)); 02265 TEST_FOR_EXCEPTION_CLASS_FUNC(missingID_gbl == 1, std::runtime_error, 02266 ": Some column IDs are not in the domain map." << std::endl 02267 << "Either these column IDs are invalid or the domain map is invalid." << std::endl 02268 << "Common cause: for a rectangular matrix, the domain map must be passed to fillComplete()."); 02269 } 02270 // Sort External column indices so that all columns coming from a given remote processor are contiguous 02271 // This obliviates the need for the Distributor associated with the Import from having to reorder data. 02272 sort2(RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin()); 02273 02274 // Copy the Local GIDs into myColumns. Two cases: 02275 // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we 02276 // can simply read the domain GIDs into the front part of ColIndices (see logic above from the serial short circuit case) 02277 // (2) We step through the GIDs of the DomainMap, checking to see if each domain GID is a column GID. 02278 // we want to do this to maintain a consistent ordering of GIDs between the columns and the domain. 02279 ArrayView<const GlobalOrdinal> mge = domainMap_->getNodeElementList(); 02280 if (numLocalColGIDs == domainMap_->getNodeNumElements()) { 02281 std::copy(mge.begin(), mge.end(), LocalColGIDs.begin()); 02282 } 02283 else { 02284 size_t numlocalagain = 0; 02285 for (size_t i=0; i < domainMap_->getNodeNumElements(); ++i) { 02286 if (GIDisLocal[i]) { 02287 LocalColGIDs[numlocalagain++] = mge[i]; 02288 } 02289 } 02290 TEST_FOR_EXCEPTION_CLASS_FUNC(numlocalagain != numLocalColGIDs, std::logic_error, ": Internal logic error. Please contact Tpetra team."); 02291 } 02292 } 02293 colMap_ = rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(GOT::invalid(), myColumns, domainMap_->getIndexBase(), domainMap_->getComm(), domainMap_->getNode()) ); 02294 checkInternalState(); 02295 return; 02296 } 02297 02298 02301 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02302 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeAllIndices() 02303 { 02304 TEST_FOR_EXCEPT( isGloballyIndexed() != false ); // this should be called only after makeIndicesLocal() 02305 TEST_FOR_EXCEPT( isSorted() != true ); // this should be called only after sortIndices() 02306 if ( isMerged() == false ) { 02307 for (size_t row=0; row < getNodeNumRows(); ++row) { 02308 RowInfo rowInfo = getRowInfo(row); 02309 mergeRowIndices(rowInfo); 02310 } 02311 // we just merged every row 02312 setMerged(true); 02313 } 02314 } 02315 02316 02319 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02320 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeImportExport() 02321 { 02322 TEST_FOR_EXCEPT(hasColMap()==false); // must have column map 02323 // create import, export 02324 if (domainMap_ != colMap_ && (!domainMap_->isSameAs(*colMap_))) { 02325 importer_ = rcp( new Import<LocalOrdinal,GlobalOrdinal,Node>(domainMap_,colMap_) ); 02326 } 02327 else { 02328 importer_ = null; 02329 } 02330 if (rangeMap_ != rowMap_ && (!rangeMap_->isSameAs(*rowMap_))) { 02331 exporter_ = rcp( new Export<LocalOrdinal,GlobalOrdinal,Node>(rowMap_,rangeMap_) ); 02332 } 02333 else { 02334 exporter_ = null; 02335 } 02336 } 02337 02338 02341 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02342 std::string CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const 02343 { 02344 std::ostringstream oss; 02345 oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description(); 02346 if (isFillComplete()) { 02347 oss << "{status = fill complete" 02348 << ", global rows = " << getGlobalNumRows() 02349 << ", global cols = " << getGlobalNumCols() 02350 << ", global num entries = " << getGlobalNumEntries() 02351 << "}"; 02352 } 02353 else { 02354 oss << "{status = fill not complete" 02355 << ", global rows = " << getGlobalNumRows() 02356 << "}"; 02357 } 02358 return oss.str(); 02359 } 02360 02361 02364 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02365 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const 02366 { 02367 using std::endl; 02368 using std::setw; 02369 using Teuchos::VERB_DEFAULT; 02370 using Teuchos::VERB_NONE; 02371 using Teuchos::VERB_LOW; 02372 using Teuchos::VERB_MEDIUM; 02373 using Teuchos::VERB_HIGH; 02374 using Teuchos::VERB_EXTREME; 02375 Teuchos::EVerbosityLevel vl = verbLevel; 02376 if (vl == VERB_DEFAULT) vl = VERB_LOW; 02377 RCP<const Comm<int> > comm = this->getComm(); 02378 const int myImageID = comm->getRank(), 02379 numImages = comm->getSize(); 02380 size_t width = 1; 02381 for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) { 02382 ++width; 02383 } 02384 width = std::max<size_t>(width,11) + 2; 02385 Teuchos::OSTab tab(out); 02386 // none: print nothing 02387 // low: print O(1) info from node 0 02388 // medium: print O(P) info, num entries per node 02389 // high: print O(N) info, num entries per row 02390 // extreme: print O(NNZ) info: print graph indices 02391 // 02392 // for medium and higher, print constituent objects at specified verbLevel 02393 if (vl != VERB_NONE) { 02394 if (myImageID == 0) out << this->description() << std::endl; 02395 // O(1) globals, minus what was already printed by description() 02396 if (isFillComplete() && myImageID == 0) { 02397 out << "Global number of diagonals = " << globalNumDiags_ << std::endl; 02398 out << "Global max number of entries = " << globalMaxNumRowEntries_ << std::endl; 02399 } 02400 // constituent objects 02401 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) { 02402 if (myImageID == 0) out << "\nRow map: " << std::endl; 02403 rowMap_->describe(out,vl); 02404 if (colMap_ != null) { 02405 if (myImageID == 0) out << "\nColumn map: " << std::endl; 02406 colMap_->describe(out,vl); 02407 } 02408 if (domainMap_ != null) { 02409 if (myImageID == 0) out << "\nDomain map: " << std::endl; 02410 domainMap_->describe(out,vl); 02411 } 02412 if (rangeMap_ != null) { 02413 if (myImageID == 0) out << "\nRange map: " << std::endl; 02414 rangeMap_->describe(out,vl); 02415 } 02416 } 02417 // O(P) data 02418 if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) { 02419 for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) { 02420 if (myImageID == imageCtr) { 02421 out << "Node ID = " << imageCtr << std::endl 02422 << "Node number of entries = " << nodeNumEntries_ << std::endl 02423 << "Node number of diagonals = " << nodeNumDiags_ << std::endl 02424 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl; 02425 if (indicesAreAllocated()) { 02426 out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl; 02427 } 02428 else { 02429 out << "Indices are not allocated." << std::endl; 02430 } 02431 } 02432 comm->barrier(); 02433 comm->barrier(); 02434 comm->barrier(); 02435 } 02436 } 02437 // O(N) and O(NNZ) data 02438 if (vl == VERB_HIGH || vl == VERB_EXTREME) { 02439 for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) { 02440 if (myImageID == imageCtr) { 02441 out << std::setw(width) << "Node ID" 02442 << std::setw(width) << "Global Row" 02443 << std::setw(width) << "Num Entries"; 02444 if (vl == VERB_EXTREME) { 02445 out << "Entries"; 02446 } 02447 out << std::endl; 02448 for (size_t r=0; r < getNodeNumRows(); ++r) { 02449 RowInfo rowinfo = getRowInfo(r); 02450 GlobalOrdinal gid = rowMap_->getGlobalElement(r); 02451 out << std::setw(width) << myImageID 02452 << std::setw(width) << gid 02453 << std::setw(width) << rowinfo.numEntries; 02454 if (vl == VERB_EXTREME) { 02455 if (isGloballyIndexed()) { 02456 ArrayView<const GlobalOrdinal> rowview = getGlobalView(rowinfo); 02457 for (size_t j=0; j < rowinfo.numEntries; ++j) out << rowview[j] << " "; 02458 } 02459 else if (isLocallyIndexed()) { 02460 ArrayView<const LocalOrdinal> rowview = getLocalView(rowinfo); 02461 for (size_t j=0; j < rowinfo.numEntries; ++j) out << colMap_->getGlobalElement(rowview[j]) << " "; 02462 } 02463 } 02464 out << std::endl; 02465 } 02466 } 02467 comm->barrier(); 02468 comm->barrier(); 02469 comm->barrier(); 02470 } 02471 } 02472 } 02473 } 02474 02475 02478 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02479 bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkSizes(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>& source) 02480 { 02481 // It's not clear what kind of compatibility checks on sizes can be performed here. 02482 // Epetra_CrsGraph doesn't check any sizes for compatibility. 02483 return true; 02484 } 02485 02486 02489 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02490 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::copyAndPermute( 02491 const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source, 02492 size_t numSameIDs, 02493 const ArrayView<const LocalOrdinal> &permuteToLIDs, 02494 const ArrayView<const LocalOrdinal> &permuteFromLIDs) 02495 { 02496 const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&>(source); 02497 std::string tfecfFuncName("copyAndPermute()"); 02498 TEST_FOR_EXCEPTION_CLASS_FUNC(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error, ": permuteToLIDs and permuteFromLIDs must have the same size."); 02499 bool src_filled = src_graph.isFillComplete(); 02500 02501 Array<GlobalOrdinal> row_copy; 02502 LocalOrdinal myid = 0; 02503 for (size_t i=0; i<numSameIDs; ++i, ++myid) { 02504 GlobalOrdinal gid = src_graph.getMap()->getGlobalElement(myid); 02505 if (src_filled) { 02506 size_t row_length = src_graph.getNumEntriesInGlobalRow(gid); 02507 row_copy.resize(row_length); 02508 size_t check_row_length = 0; 02509 src_graph.getGlobalRowCopy(gid, row_copy(), check_row_length); 02510 insertGlobalIndices( gid, row_copy() ); 02511 } 02512 else { 02513 ArrayView<const GlobalOrdinal> row; 02514 src_graph.getGlobalRowView( gid,row ); 02515 insertGlobalIndices( gid, row ); 02516 } 02517 } 02518 02519 for (LocalOrdinal i=0; i<permuteToLIDs.size(); ++i) { 02520 GlobalOrdinal mygid = this->getMap()->getGlobalElement(permuteToLIDs[i]); 02521 GlobalOrdinal srcgid= src_graph.getMap()->getGlobalElement(permuteFromLIDs[i]); 02522 if (src_filled) { 02523 size_t row_length = src_graph.getNumEntriesInGlobalRow(srcgid); 02524 row_copy.resize(row_length); 02525 size_t check_row_length = 0; 02526 src_graph.getGlobalRowCopy(srcgid, row_copy(), check_row_length); 02527 insertGlobalIndices( mygid, row_copy() ); 02528 } 02529 else { 02530 ArrayView<const GlobalOrdinal> row; 02531 src_graph.getGlobalRowView( srcgid, row ); 02532 insertGlobalIndices( mygid, row ); 02533 } 02534 } 02535 } 02536 02537 02540 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02541 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::packAndPrepare( 02542 const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source, 02543 const ArrayView<const LocalOrdinal> &exportLIDs, 02544 Array<GlobalOrdinal> &exports, 02545 const ArrayView<size_t> & numPacketsPerLID, 02546 size_t& constantNumPackets, 02547 Distributor &distor) 02548 { 02549 std::string tfecfFuncName("packAndPrepare()"); 02550 TEST_FOR_EXCEPTION_CLASS_FUNC(exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error, 02551 ": exportLIDs and numPacketsPerLID must have the same size."); 02552 const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&>(source); 02553 // We don't check whether src_graph has had fillComplete called, because it doesn't matter whether the 02554 // *source* graph has been fillComplete'd. The target graph can not be fillComplete'd yet. 02555 TEST_FOR_EXCEPTION_CLASS_FUNC(this->isFillComplete() == true, std::runtime_error, 02556 ": import/export operations are not allowed on destination CrsGraph after fillComplete has been called."); 02557 constantNumPackets = 0; 02558 // first set the contents of numPacketsPerLID, and accumulate a total-num-packets: 02559 size_t totalNumPackets = 0; 02560 Array<GlobalOrdinal> row; 02561 for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) { 02562 GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]); 02563 size_t row_length = src_graph.getNumEntriesInGlobalRow(GID); 02564 numPacketsPerLID[i] = row_length; 02565 totalNumPackets += row_length; 02566 } 02567 02568 exports.resize(totalNumPackets); 02569 02570 // now loop again and pack rows of indices into exports: 02571 size_t exportsOffset = 0; 02572 for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) { 02573 GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]); 02574 size_t row_length = src_graph.getNumEntriesInGlobalRow(GID); 02575 row.resize(row_length); 02576 size_t check_row_length = 0; 02577 src_graph.getGlobalRowCopy(GID, row(), check_row_length); 02578 typename Array<GlobalOrdinal>::const_iterator 02579 row_iter = row.begin(), row_end = row.end(); 02580 size_t j = 0; 02581 for (; row_iter != row_end; ++row_iter, ++j) { 02582 exports[exportsOffset+j] = *row_iter; 02583 } 02584 exportsOffset += row.size(); 02585 } 02586 } 02587 02588 02591 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02592 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::unpackAndCombine( 02593 const ArrayView<const LocalOrdinal> &importLIDs, 02594 const ArrayView<const GlobalOrdinal> &imports, 02595 const ArrayView<size_t> &numPacketsPerLID, 02596 size_t constantNumPackets, 02597 Distributor & /* distor */, 02598 CombineMode /* CM */) 02599 { 02600 // We are not checking the value of the CombineMode input-argument. 02601 // For CrsGraph, we only support import/export operations if fillComplete has not yet been called. 02602 // Any incoming column-indices are inserted into the target graph. In this context, CombineMode values 02603 // of ADD vs INSERT are equivalent. What is the meaning of REPLACE for CrsGraph? If a duplicate column-index 02604 // is inserted, it will be compressed out when fillComplete is called. 02605 // 02606 // Note: I think REPLACE means that an existing row is replaced by the imported row, i.e., the existing indices are cleared. CGB, 6/17/2010 02607 02608 std::string tfecfFuncName("unpackAndCombine()"); 02609 TEST_FOR_EXCEPTION_CLASS_FUNC(importLIDs.size() != numPacketsPerLID.size(), std::runtime_error, 02610 ": importLIDs and numPacketsPerLID must have the same size."); 02611 TEST_FOR_EXCEPTION_CLASS_FUNC(this->isFillComplete() == true, std::runtime_error, 02612 ": import/export operations are not allowed on destination CrsGraph after fillComplete has been called."); 02613 size_t importsOffset = 0; 02614 typename ArrayView<const LocalOrdinal>::iterator 02615 impLIDiter = importLIDs.begin(), impLIDend = importLIDs.end(); 02616 size_t i = 0; 02617 for (; impLIDiter != impLIDend; ++impLIDiter, ++i) { 02618 LocalOrdinal row_length = numPacketsPerLID[i]; 02619 const ArrayView<const GlobalOrdinal> row(&imports[importsOffset], row_length); 02620 insertGlobalIndices(this->getMap()->getGlobalElement(*impLIDiter), row); 02621 importsOffset += row_length; 02622 } 02623 } 02624 02625 02628 // // 02629 // Deprecated methods // 02630 // // 02633 02634 02636 // DEPRECATED 02638 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02639 ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowView(LocalOrdinal localRow) const 02640 { 02641 std::string tfecfFuncName("getLocalRowView(localRow)"); 02642 TEST_FOR_EXCEPTION_CLASS_FUNC(isGloballyIndexed() == true, std::runtime_error, 02643 ": local indices do not exist."); 02644 TEST_FOR_EXCEPTION_CLASS_FUNC(rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, 02645 ": localRow (== " << localRow << ") is not valid on this node."); 02646 ArrayRCP<const LocalOrdinal> ret; 02647 RowInfo sizeInfo = getRowInfo(localRow); 02648 if (getProfileType() == StaticProfile && sizeInfo.numEntries > 0) { 02649 ret = lclInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.numEntries); 02650 } 02651 else if (getProfileType() == DynamicProfile && sizeInfo.numEntries > 0) { 02652 ret = lclInds2D_[localRow].persistingView(0,sizeInfo.numEntries); 02653 } 02654 return ret; 02655 } 02656 02657 02659 // DEPRECATED 02661 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02662 ArrayRCP<const GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowView(GlobalOrdinal globalRow) const 02663 { 02664 std::string tfecfFuncName("getGlobalRowView(globalRow)"); 02665 TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, 02666 ": global indices do not exist."); 02667 const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow); 02668 TEST_FOR_EXCEPTION_CLASS_FUNC(lrow == OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error, 02669 ": globalRow (== " << globalRow << ") does not belong to this node."); 02670 ArrayRCP<const GlobalOrdinal> ret; 02671 RowInfo sizeInfo = getRowInfo(lrow); 02672 if (getProfileType() == StaticProfile && sizeInfo.numEntries > 0) { 02673 ret = gblInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.numEntries); 02674 } 02675 else if (getProfileType() == DynamicProfile && sizeInfo.numEntries > 0) { 02676 ret = gblInds2D_[lrow].persistingView(0,sizeInfo.numEntries); 02677 } 02678 return ret; 02679 } 02680 02681 02683 // DEPRECATED 02685 template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps> 02686 void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::optimizeStorage() 02687 { 02688 // provided only for backwards compatibility 02689 // previous semantics required that fillComplete() had been called. 02690 std::string tfecfFuncName("optimizeStorage()"); 02691 TEST_FOR_EXCEPTION_CLASS_FUNC(isFillComplete() == false, std::runtime_error, ": requires that fillComplete() has already been called."); 02692 if (isStorageOptimized() == false) { 02693 resumeFill(); 02694 fillComplete(DoOptimizeStorage); 02695 } 02696 } 02697 02698 02699 } // namespace Tpetra 02700 02701 02702 // 02703 // Explicit instantiation macro 02704 // 02705 // Must be expanded from within the Tpetra namespace! 02706 // 02707 02708 #define TPETRA_CRSGRAPH_INSTANT(LO,GO,NODE) \ 02709 \ 02710 template class CrsGraph< LO , GO , NODE >; \ 02711 02712 02713 #endif // TPETRA_CRSGRAPH_DEF_HPP
1.7.4
