|
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 #ifndef TPETRA_MULTIVECTOR_DEF_HPP 00030 #define TPETRA_MULTIVECTOR_DEF_HPP 00031 00032 #include <Kokkos_NodeTrace.hpp> 00033 00034 #include <Teuchos_TestForException.hpp> 00035 #include <Teuchos_as.hpp> 00036 00037 #include "Tpetra_Vector.hpp" 00038 00039 #ifdef DOXYGEN_USE_ONLY 00040 #include "Tpetra_MultiVector_decl.hpp" 00041 #endif 00042 00043 namespace Tpetra { 00044 00045 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00046 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector( 00047 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 00048 size_t NumVectors, 00049 bool zeroOut /* default is true */ 00050 ) 00051 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), 00052 lclMV_(map->getNode()) { 00053 TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument, 00054 "Tpetra::MultiVector::MultiVector(): NumVectors must be strictly positive."); 00055 const size_t myLen = getLocalLength(); 00056 if (myLen > 0) { 00057 Teuchos::RCP<Node> node = map->getNode(); 00058 Teuchos::ArrayRCP<Scalar> data = node->template allocBuffer<Scalar>(myLen*NumVectors); 00059 MVT::initializeValues(lclMV_,myLen,NumVectors,data,myLen); 00060 if (zeroOut) { 00061 MVT::Init(lclMV_, Teuchos::ScalarTraits<Scalar>::zero()); 00062 } 00063 } 00064 else { 00065 MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0); 00066 } 00067 } 00068 00069 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00070 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source) 00071 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(source), lclMV_(MVT::getNode(source.lclMV_)) { 00072 // copy data from the source MultiVector into this multivector 00073 Teuchos::RCP<Node> node = MVT::getNode(source.lclMV_); 00074 const LocalOrdinal myLen = getLocalLength(); 00075 const size_t numVecs = source.getNumVectors(); 00076 if (myLen > 0) { 00077 // allocate data 00078 Teuchos::ArrayRCP<Scalar> data = node->template allocBuffer<Scalar>(myLen*numVecs); 00079 MVT::initializeValues(lclMV_,myLen,numVecs,data,myLen); 00080 // copy data 00081 { 00082 Teuchos::ArrayRCP<Scalar> dstdata = data; 00083 Teuchos::ArrayRCP<const Scalar> srcdata = MVT::getValues(source.lclMV_); 00084 for (size_t j = 0; j < numVecs; ++j) { 00085 Teuchos::ArrayRCP<const Scalar> srcj = source.getSubArrayRCP(srcdata,j); 00086 KOKKOS_NODE_TRACE("MultiVector::MultiVector(MV)") 00087 node->template copyBuffers<Scalar>(myLen,srcj,dstdata); 00088 dstdata += myLen; 00089 } 00090 } 00091 } 00092 else { 00093 MVT::initializeValues(lclMV_,0,numVecs,Teuchos::null,0); 00094 } 00095 } 00096 00097 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00098 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector( 00099 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 00100 const Teuchos::ArrayView<const Scalar> &A, size_t LDA, 00101 size_t NumVectors) 00102 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) { 00103 TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument, 00104 "Tpetra::MultiVector::MultiVector(A,LDA): NumVectors must be strictly positive."); 00105 const size_t myLen = getLocalLength(); 00106 using Teuchos::ArrayRCP; 00107 TEST_FOR_EXCEPTION(LDA < myLen, std::runtime_error, 00108 "Tpetra::MultiVector::MultiVector(A,LDA): LDA must be large enough to accomodate the local entries."); 00109 #ifdef HAVE_TPETRA_DEBUG 00110 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(A.size()) < LDA*(NumVectors-1)+myLen, std::runtime_error, 00111 "Tpetra::MultiVector::MultiVector(A,LDA): A does not contain enough data to specify the entries in this."); 00112 #endif 00113 if (myLen > 0) { 00114 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00115 Teuchos::ArrayRCP<Scalar> mydata = node->template allocBuffer<Scalar>(myLen*NumVectors); 00116 MVT::initializeValues(lclMV_,myLen,NumVectors,mydata,myLen); 00117 KOKKOS_NODE_TRACE("MultiVector::MultiVector(1D)") 00118 Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::WriteOnly,myLen*NumVectors,mydata); 00119 typename Teuchos::ArrayView<const Scalar>::iterator srcit = A.begin(); 00120 for (size_t j = 0; j < NumVectors; ++j) { 00121 std::copy(srcit,srcit+myLen,myview); 00122 srcit += LDA; 00123 myview += myLen; 00124 } 00125 mydata = Teuchos::null; 00126 myview = Teuchos::null; 00127 } 00128 else { 00129 MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0); 00130 } 00131 } 00132 00133 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00134 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 00135 const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> > &ArrayOfPtrs, 00136 size_t NumVectors) 00137 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) { 00138 TEST_FOR_EXCEPTION(NumVectors < 1 || NumVectors != Teuchos::as<size_t>(ArrayOfPtrs.size()), std::runtime_error, 00139 "Tpetra::MultiVector::MultiVector(ArrayOfPtrs): ArrayOfPtrs.size() must be strictly positive and as large as ArrayOfPtrs."); 00140 const size_t myLen = getLocalLength(); 00141 if (myLen > 0) { 00142 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00143 Teuchos::ArrayRCP<Scalar> mydata = node->template allocBuffer<Scalar>(myLen*NumVectors); 00144 MVT::initializeValues(lclMV_,myLen,NumVectors,mydata,myLen); 00145 KOKKOS_NODE_TRACE("MultiVector::MultiVector(2D)") 00146 Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::WriteOnly,myLen*NumVectors,mydata); 00147 for (size_t j = 0; j < NumVectors; ++j) { 00148 #ifdef HAVE_TPETRA_DEBUG 00149 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != getLocalLength(), std::runtime_error, 00150 "Tpetra::MultiVector::MultiVector(ArrayOfPtrs): ArrayOfPtrs[" << j << "].size() (== " << ArrayOfPtrs[j].size() << 00151 ") is not equal to getLocalLength() (== " << getLocalLength()); 00152 #endif 00153 typename Teuchos::ArrayView<const Scalar>::iterator src = ArrayOfPtrs[j].begin(); 00154 std::copy(src,src+myLen,myview); 00155 myview += myLen; 00156 } 00157 myview = Teuchos::null; 00158 mydata = Teuchos::null; 00159 } 00160 else { 00161 MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0); 00162 } 00163 } 00164 00165 00166 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00167 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 00168 Teuchos::ArrayRCP<Scalar> data, size_t LDA, Teuchos::ArrayView<const size_t> WhichVectors) 00169 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()), whichVectors_(WhichVectors) { 00170 const size_t myLen = getLocalLength(); 00171 size_t maxVector = *std::max_element(WhichVectors.begin(), WhichVectors.end()); 00172 TEST_FOR_EXCEPTION(LDA < myLen, std::runtime_error, 00173 "Tpetra::MultiVector::MultiVector(data,LDA,WhichVectors): LDA must be large enough to accomodate the local entries."); 00174 #ifdef HAVE_TPETRA_DEBUG 00175 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(data.size()) < LDA * maxVector + myLen, std::runtime_error, 00176 "Tpetra::MultiVector::MultiVector(data,LDA,WhichVectors): data does not contain enough data to specify the entries in this."); 00177 #endif 00178 if (WhichVectors.size() == 1) { 00179 // shift data so that desired vector is vector 0 00180 maxVector = 0; 00181 data += LDA*WhichVectors[0]; 00182 // kill whichVectors_; we are constant stride 00183 whichVectors_.clear(); 00184 } 00185 if (myLen > 0) { 00186 MVT::initializeValues(lclMV_,myLen,maxVector+1,data,LDA); 00187 } 00188 else { 00189 MVT::initializeValues(lclMV_,0,WhichVectors.size(),Teuchos::null,0); 00190 } 00191 } 00192 00193 00194 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00195 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 00196 Teuchos::ArrayRCP<Scalar> data, size_t LDA, size_t NumVectors) 00197 : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) { 00198 TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument, 00199 "Tpetra::MultiVector::MultiVector(data,LDA,NumVector): NumVectors must be strictly positive."); 00200 const LocalOrdinal myLen = getLocalLength(); 00201 #ifdef HAVE_TPETRA_DEBUG 00202 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(data.size()) < LDA*(NumVectors-1)+myLen, std::runtime_error, 00203 "Tpetra::MultiVector::MultiVector(data,LDA,NumVector): data does not contain enough data to specify the entries in this."); 00204 #endif 00205 MVT::initializeValues(lclMV_,myLen,NumVectors,data,LDA); 00206 } 00207 00208 00209 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00210 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~MultiVector() { 00211 } 00212 00213 00214 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00215 bool MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isConstantStride() const { 00216 return whichVectors_.empty(); 00217 } 00218 00219 00220 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00221 size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalLength() const { 00222 return this->getMap()->getNodeNumElements(); 00223 } 00224 00225 00226 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00227 global_size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getGlobalLength() const { 00228 return this->getMap()->getGlobalNumElements(); 00229 } 00230 00231 00232 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00233 size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getStride() const { 00234 if (isConstantStride()) { 00235 return MVT::getStride(lclMV_); 00236 } 00237 return 0; 00238 } 00239 00240 00241 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00242 bool MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::checkSizes(const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceObj) 00243 { 00244 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(sourceObj); 00245 // objects maps have already been checked. simply check the number of vectors. 00246 bool compat = (A.getNumVectors() == this->getNumVectors()); 00247 return compat; 00248 } 00249 00250 00251 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00252 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::copyAndPermute( 00253 const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> & sourceObj, 00254 size_t numSameIDs, 00255 const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs, 00256 const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs) { 00257 using Teuchos::ArrayView; 00258 using Teuchos::ArrayRCP; 00259 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceMV = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &>(sourceObj); 00260 typename ArrayView<const LocalOrdinal>::iterator pTo, pFrom; 00261 // any other error will be caught by Teuchos 00262 TEST_FOR_EXCEPTION(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error, 00263 "Tpetra::MultiVector::copyAndPermute(): permuteToLIDs and permuteFromLIDs must have the same size."); 00264 // one vector at a time 00265 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00266 const size_t numCols = getNumVectors(); 00267 // Get a host view of the local multivector data. 00268 for (size_t j = 0; j < numCols; ++j) { 00269 // The first numImportIDs GIDs are the same between source and target, 00270 // We can just copy them 00271 ArrayRCP<const Scalar> srcptr = sourceMV.getSubArrayRCP(sourceMV.cview_,j); 00272 ArrayRCP< Scalar> dstptr = getSubArrayRCP(ncview_,j); 00273 std::copy(srcptr,srcptr+numSameIDs,dstptr); 00274 // next, do permutations 00275 for (pTo = permuteToLIDs.begin(), pFrom = permuteFromLIDs.begin(); 00276 pTo != permuteToLIDs.end(); ++pTo, ++pFrom) { 00277 dstptr[*pTo] = srcptr[*pFrom]; 00278 } 00279 } 00280 } 00281 00282 00283 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00284 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::packAndPrepare( 00285 const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> & sourceObj, 00286 const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs, 00287 Teuchos::Array<Scalar> &exports, 00288 const Teuchos::ArrayView<size_t> &numExportPacketsPerLID, 00289 size_t& constantNumPackets, 00290 Distributor & /* distor */ ) { 00291 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceMV = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &>(sourceObj); 00292 using Teuchos::ArrayView; 00293 using Teuchos::ArrayRCP; 00294 /* The layout in the export for MultiVectors is as follows: 00295 exports = { all of the data from row exportLIDs.front() ; 00296 .... 00297 all of the data from row exportLIDs.back() } 00298 this doesn't have the best locality, but is necessary because the data for a Packet 00299 (all data associated with an LID) is required to be contiguous */ 00300 TEST_FOR_EXCEPTION(Teuchos::as<int>(numExportPacketsPerLID.size()) != exportLIDs.size(), std::runtime_error, 00301 "Tpetra::MultiVector::packAndPrepare(): size of numExportPacketsPerLID buffer should be the same as exportLIDs."); 00302 const KMV &srcData = sourceMV.lclMV_; 00303 const size_t numCols = sourceMV.getNumVectors(), 00304 stride = MVT::getStride(srcData); 00305 constantNumPackets = numCols; 00306 exports.resize(numCols*exportLIDs.size()); 00307 typename ArrayView<const LocalOrdinal>::iterator idptr; 00308 typename Teuchos::Array<Scalar>::iterator expptr; 00309 expptr = exports.begin(); 00310 00311 Teuchos::RCP<Node> node = MVT::getNode(srcData); 00312 if (sourceMV.isConstantStride()) { 00313 size_t i = 0; 00314 for (idptr = exportLIDs.begin(); idptr != exportLIDs.end(); ++idptr, ++i) { 00315 for (size_t j = 0; j < numCols; ++j) { 00316 *expptr++ = sourceMV.cview_[j*stride + (*idptr)]; 00317 } 00318 //we shouldn't need to set numExportPacketsPerLID[i] since we have set 00319 //constantNumPackets to a nonzero value. But we'll set it anyway, since 00320 //I'm not sure if the API will remain the way it is. 00321 numExportPacketsPerLID[i] = numCols; 00322 } 00323 } 00324 else { 00325 size_t i = 0; 00326 for (idptr = exportLIDs.begin(); idptr != exportLIDs.end(); ++idptr, ++i) { 00327 for (size_t j = 0; j < numCols; ++j) { 00328 *expptr++ = sourceMV.cview_[sourceMV.whichVectors_[j]*stride + (*idptr)]; 00329 } 00330 numExportPacketsPerLID[i] = numCols; 00331 } 00332 } 00333 } 00334 00335 00336 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00337 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unpackAndCombine( 00338 const Teuchos::ArrayView<const LocalOrdinal> &importLIDs, 00339 const Teuchos::ArrayView<const Scalar> &imports, 00340 const Teuchos::ArrayView<size_t> &numPacketsPerLID, 00341 size_t constantNumPackets, 00342 Distributor & /* distor */, 00343 CombineMode CM) { 00344 using Teuchos::ArrayView; 00345 using Teuchos::ArrayRCP; 00346 /* The layout in the export for MultiVectors is as follows: 00347 imports = { all of the data from row exportLIDs.front() ; 00348 .... 00349 all of the data from row exportLIDs.back() } 00350 this doesn't have the best locality, but is necessary because the data for a Packet 00351 (all data associated with an LID) is required to be contiguous */ 00352 #ifdef HAVE_TPETRA_DEBUG 00353 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(imports.size()) != getNumVectors()*importLIDs.size(), std::runtime_error, 00354 "Tpetra::MultiVector::unpackAndCombine(): sizing of imports buffer should be appropriate for the amount of data to be exported."); 00355 #endif 00356 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(constantNumPackets) == 0u, std::runtime_error, 00357 "Tpetra::MultiVector::unpackAndCombine(): 'constantNumPackets' input argument should be nonzero."); 00358 00359 const size_t myStride = MVT::getStride(lclMV_), 00360 numVecs = getNumVectors(); 00361 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(numPacketsPerLID.size()) != Teuchos::as<size_t>(importLIDs.size()), std::runtime_error, 00362 "Tpetra::MultiVector::unpackAndCombine(): 'numPacketsPerLID' must have same length as importLIDs."); 00363 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(numVecs) != Teuchos::as<size_t>(constantNumPackets), std::runtime_error, 00364 "Tpetra::MultiVector::unpackAndCombine(): 'constantNumPackets' must equal numVecs."); 00365 00366 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00367 if (numVecs > 0 && importLIDs.size()) { 00368 typename ArrayView<const Scalar>::iterator impptr; 00369 typename ArrayView<const LocalOrdinal>::iterator idptr; 00370 impptr = imports.begin(); 00371 if (CM == INSERT || CM == REPLACE) { 00372 if (isConstantStride()) { 00373 for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) { 00374 for (size_t j = 0; j < numVecs; ++j) { 00375 ncview_[myStride*j + *idptr] = *impptr++; 00376 } 00377 } 00378 } 00379 else { 00380 for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) { 00381 for (size_t j = 0; j < numVecs; ++j) { 00382 ncview_[myStride*whichVectors_[j] + *idptr] = *impptr++; 00383 } 00384 } 00385 } 00386 } 00387 else if (CM == ADD) { 00388 if (isConstantStride()) { 00389 for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) { 00390 for (size_t j = 0; j < numVecs; ++j) { 00391 ncview_[myStride*j + *idptr] += *impptr++; 00392 } 00393 } 00394 } 00395 else { 00396 for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) { 00397 for (size_t j = 0; j < numVecs; ++j) { 00398 ncview_[myStride*whichVectors_[j] + *idptr] += *impptr++; 00399 } 00400 } 00401 } 00402 } 00403 else { 00404 TEST_FOR_EXCEPTION(CM != ADD && CM != REPLACE && CM != INSERT, std::invalid_argument, 00405 "Tpetra::MultiVector::unpackAndCombine(): Invalid CombineMode: " << CM); 00406 } 00407 } 00408 } 00409 00410 00411 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00412 inline size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNumVectors() const { 00413 size_t ret; 00414 if (isConstantStride()) { 00415 ret = MVT::getNumCols(lclMV_); 00416 } 00417 else { 00418 ret = whichVectors_.size(); 00419 } 00420 return ret; 00421 } 00422 00423 00424 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00425 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dot( 00426 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 00427 const Teuchos::ArrayView<Scalar> &dots) const 00428 { 00429 using Teuchos::ArrayRCP; 00430 using Teuchos::arcp_const_cast; 00431 const size_t myLen = getLocalLength(), 00432 numVecs = getNumVectors(); 00433 #ifdef HAVE_TPETRA_DEBUG 00434 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error, 00435 "Tpetra::MultiVector::dot(): MultiVectors do not have compatible Maps:" << std::endl 00436 << "this->getMap(): " << std::endl << *this->getMap() 00437 << "A.getMap(): " << std::endl << *A.getMap() << std::endl); 00438 #else 00439 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error, 00440 "Tpetra::MultiVector::dot(): MultiVectors do not have the same local length."); 00441 #endif 00442 TEST_FOR_EXCEPTION(A.getNumVectors() != numVecs, std::runtime_error, 00443 "Tpetra::MultiVector::dot(): MultiVectors must have the same number of vectors."); 00444 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dots.size()) != numVecs, std::runtime_error, 00445 "Tpetra::MultiVector::dot(): dots.size() must be as large as the number of vectors in *this and A."); 00446 if (isConstantStride() && A.isConstantStride()) { 00447 MVT::Dot(lclMV_,A.lclMV_,dots); 00448 } 00449 else { 00450 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)); 00451 ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_)), 00452 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)); 00453 for (size_t j=0; j < numVecs; ++j) { 00454 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00455 avj = A.getSubArrayRCP(avptr,j); 00456 MVT::initializeValues(a,myLen, 1, avj, myLen); 00457 MVT::initializeValues(v,myLen, 1, vj, myLen); 00458 dots[j] = MVT::Dot((const KMV&)v,(const KMV &)a); 00459 } 00460 } 00461 if (this->isDistributed()) { 00462 Teuchos::Array<Scalar> ldots(dots); 00463 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),ldots.getRawPtr(),dots.getRawPtr()); 00464 } 00465 } 00466 00467 00468 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00469 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm2( 00470 const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const { 00471 using Teuchos::ScalarTraits; 00472 using Teuchos::ArrayView; 00473 typedef typename ScalarTraits<Scalar>::magnitudeType Mag; 00474 const size_t numVecs = this->getNumVectors(); 00475 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error, 00476 "Tpetra::MultiVector::norm2(norms): norms.size() must be as large as the number of vectors in *this."); 00477 if (isConstantStride()) { 00478 MVT::Norm2Squared(lclMV_,norms); 00479 } 00480 else { 00481 KMV v(MVT::getNode(lclMV_)); 00482 Teuchos::ArrayRCP<Scalar> vi; 00483 for (size_t i=0; i < numVecs; ++i) { 00484 vi = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[i]) ); 00485 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vi, MVT::getStride(lclMV_)); 00486 norms[i] = MVT::Norm2Squared(v); 00487 } 00488 } 00489 if (this->isDistributed()) { 00490 Teuchos::Array<Mag> lnorms(norms); 00491 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr()); 00492 } 00493 for (typename ArrayView<Mag>::iterator n = norms.begin(); n != norms.begin()+numVecs; ++n) { 00494 (*n) = ScalarTraits<Mag>::squareroot(*n); 00495 } 00496 } 00497 00498 00499 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00500 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normWeighted( 00501 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights, 00502 const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const { 00503 using Teuchos::ScalarTraits; 00504 using Teuchos::ArrayView; 00505 using Teuchos::ArrayRCP; 00506 using Teuchos::arcp_const_cast; 00507 typedef ScalarTraits<Scalar> SCT; 00508 typedef typename SCT::magnitudeType Mag; 00509 const Mag OneOverN = ScalarTraits<Mag>::one() / Teuchos::as<Mag>(getGlobalLength()); 00510 bool OneW = false; 00511 const size_t numVecs = this->getNumVectors(); 00512 if (weights.getNumVectors() == 1) { 00513 OneW = true; 00514 } 00515 else { 00516 TEST_FOR_EXCEPTION(weights.getNumVectors() != numVecs, std::runtime_error, 00517 "Tpetra::MultiVector::normWeighted(): MultiVector of weights must contain either one vector or the same number of vectors as this."); 00518 } 00519 #ifdef HAVE_TPETRA_DEBUG 00520 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*weights.getMap()), std::runtime_error, 00521 "Tpetra::MultiVector::normWeighted(): MultiVectors do not have compatible Maps:" << std::endl 00522 << "this->getMap(): " << std::endl << *this->getMap() 00523 << "weights.getMap(): " << std::endl << *weights.getMap() << std::endl); 00524 #else 00525 TEST_FOR_EXCEPTION( getLocalLength() != weights.getLocalLength(), std::runtime_error, 00526 "Tpetra::MultiVector::normWeighted(): MultiVectors do not have the same local length."); 00527 #endif 00528 const size_t myLen = getLocalLength(); 00529 // 00530 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error, 00531 "Tpetra::MultiVector::normWeighted(): norms.size() must be as large as the number of vectors in *this."); 00532 if (isConstantStride() && weights.isConstantStride()) { 00533 MVT::WeightedNorm(lclMV_,weights.lclMV_,norms); 00534 } 00535 else { 00536 KMV v(MVT::getNode(lclMV_)), w(MVT::getNode(lclMV_)); 00537 ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_)), 00538 wptr = arcp_const_cast<Scalar>(MVT::getValues(weights.lclMV_)); 00539 ArrayRCP<Scalar> wj = wptr.persistingView(0,myLen); 00540 MVT::initializeValues(w,myLen, 1, wj, myLen); 00541 for (size_t j=0; j < numVecs; ++j) { 00542 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j); 00543 MVT::initializeValues(v,myLen, 1, vj, myLen); 00544 if (!OneW) { 00545 wj = weights.getSubArrayRCP(wptr,j); 00546 MVT::initializeValues(w,myLen, 1, wj, myLen); 00547 } 00548 norms[j] = MVT::WeightedNorm((const KMV&)v,(const KMV &)w); 00549 } 00550 } 00551 if (this->isDistributed()) { 00552 Teuchos::Array<Mag> lnorms(norms); 00553 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr()); 00554 } 00555 for (typename ArrayView<Mag>::iterator n = norms.begin(); n != norms.begin()+numVecs; ++n) { 00556 *n = ScalarTraits<Mag>::squareroot(*n * OneOverN); 00557 } 00558 } 00559 00560 00561 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00562 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm1( 00563 const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const { 00564 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag; 00565 const size_t numVecs = this->getNumVectors(); 00566 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error, 00567 "Tpetra::MultiVector::norm1(norms): norms.size() must be as large as the number of vectors in *this."); 00568 if (isConstantStride()) { 00569 MVT::Norm1(lclMV_,norms); 00570 } 00571 else { 00572 KMV v(MVT::getNode(lclMV_)); 00573 Teuchos::ArrayRCP<Scalar> vj; 00574 for (size_t j=0; j < numVecs; ++j) { 00575 vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) ); 00576 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_)); 00577 norms[j] = MVT::Norm1((const KMV&)v); 00578 } 00579 } 00580 if (this->isDistributed()) { 00581 Teuchos::Array<Mag> lnorms(norms); 00582 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr()); 00583 } 00584 } 00585 00586 00587 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00588 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normInf( 00589 const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const { 00590 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag; 00591 const size_t numVecs = this->getNumVectors(); 00592 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error, 00593 "Tpetra::MultiVector::normInf(norms): norms.size() must be as large as the number of vectors in *this."); 00594 if (isConstantStride()) { 00595 MVT::NormInf(lclMV_,norms); 00596 } 00597 else { 00598 KMV v(MVT::getNode(lclMV_)); 00599 Teuchos::ArrayRCP<Scalar> vj; 00600 for (size_t j=0; j < numVecs; ++j) { 00601 vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) ); 00602 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_)); 00603 norms[j] = MVT::NormInf((const KMV&)v); 00604 } 00605 } 00606 if (this->isDistributed()) { 00607 Teuchos::Array<Mag> lnorms(norms); 00608 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_MAX,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr()); 00609 } 00610 } 00611 00612 00613 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00614 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::meanValue(const Teuchos::ArrayView<Scalar> &means) const 00615 { 00616 typedef Teuchos::ScalarTraits<Scalar> SCT; 00617 using Teuchos::ArrayRCP; 00618 using Teuchos::arcp_const_cast; 00619 const size_t numVecs = getNumVectors(); 00620 const size_t myLen = getLocalLength(); 00621 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(means.size()) != numVecs, std::runtime_error, 00622 "Tpetra::MultiVector::meanValue(): means.size() must be as large as the number of vectors in *this."); 00623 // compute local components of the means 00624 // sum these across all nodes 00625 if (isConstantStride()) { 00626 MVT::Sum(lclMV_,means); 00627 } 00628 else { 00629 KMV v(MVT::getNode(lclMV_)); 00630 ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_)); 00631 for (size_t j=0; j < numVecs; ++j) { 00632 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j); 00633 MVT::initializeValues(v,myLen, 1, vj, myLen); 00634 means[j] = MVT::Sum((const KMV &)v); 00635 } 00636 } 00637 if (this->isDistributed()) { 00638 Teuchos::Array<Scalar> lmeans(means); 00639 // only combine if we are a distributed MV 00640 Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lmeans.getRawPtr(),means.getRawPtr()); 00641 } 00642 const Scalar OneOverN = Teuchos::ScalarTraits<Scalar>::one() / Teuchos::as<Scalar>(getGlobalLength()); 00643 for (typename Teuchos::ArrayView<Scalar>::iterator i = means.begin(); i != means.begin()+numVecs; ++i) { 00644 (*i) = (*i)*OneOverN; 00645 } 00646 } 00647 00648 00649 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00650 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::randomize() { 00651 if (isConstantStride()) { 00652 MVT::Random(lclMV_); 00653 } 00654 else { 00655 const size_t numVecs = this->getNumVectors(); 00656 KMV v(MVT::getNode(lclMV_)); 00657 Teuchos::ArrayRCP<Scalar> vj; 00658 for (size_t j=0; j < numVecs; ++j) { 00659 vj = MVT::getValuesNonConst(lclMV_,whichVectors_[j]); 00660 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_)); 00661 MVT::Random(v); 00662 } 00663 } 00664 } 00665 00666 00667 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00668 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::putScalar(const Scalar &alpha) { 00669 const size_t numVecs = getNumVectors(); 00670 if (isConstantStride()) { 00671 MVT::Init(lclMV_,alpha); 00672 } 00673 else { 00674 KMV v(MVT::getNode(lclMV_)); 00675 Teuchos::ArrayRCP<Scalar> vj; 00676 for (size_t j=0; j < numVecs; ++j) { 00677 vj = MVT::getValuesNonConst(lclMV_,whichVectors_[j]); 00678 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_)); 00679 MVT::Init(v,alpha); 00680 } 00681 } 00682 } 00683 00684 00685 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00686 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(const Scalar &alpha) { 00687 // NOTE: can't substitute putScalar(0.0) for scale(0.0), because 00688 // the former will overwrite NaNs present in the MultiVector, while the 00689 // semantics of this call require multiplying them by 0, which IEEE requires to be NaN 00690 const size_t numVecs = getNumVectors(); 00691 if (alpha == Teuchos::ScalarTraits<Scalar>::one()) { 00692 // do nothing 00693 } 00694 else if (isConstantStride()) { 00695 MVT::Scale(lclMV_,alpha); 00696 } 00697 else { 00698 KMV v(MVT::getNode(lclMV_)); 00699 Teuchos::ArrayRCP<Scalar> vj; 00700 for (size_t j=0; j < numVecs; ++j) { 00701 vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) ); 00702 MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_)); 00703 MVT::Scale(v,alpha); 00704 } 00705 } 00706 } 00707 00708 00709 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00710 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(Teuchos::ArrayView<const Scalar> alphas) 00711 { 00712 using Teuchos::ArrayRCP; 00713 const size_t numVecs = this->getNumVectors(); 00714 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(alphas.size()) != numVecs, std::runtime_error, 00715 "Tpetra::MultiVector::scale(alphas): alphas.size() must be as large as the number of vectors in *this."); 00716 KMV vec(MVT::getNode(lclMV_)); 00717 const size_t myLen = MVT::getNumRows(lclMV_); 00718 if (myLen == 0) return; 00719 ArrayRCP<Scalar> mybuf = MVT::getValuesNonConst(lclMV_); 00720 for (size_t j = 0; j < numVecs; ++j) { 00721 if (alphas[j] == Teuchos::ScalarTraits<Scalar>::one()) { 00722 // do nothing: NaN * 1.0 == NaN, Number*1.0 == Number 00723 } 00724 else { 00725 ArrayRCP<Scalar> mybufj = getSubArrayRCP(mybuf,j); 00726 MVT::initializeValues(vec,myLen,1,mybufj,myLen); 00727 MVT::Scale(vec,alphas[j]); 00728 } 00729 } 00730 } 00731 00732 00733 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00734 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) { 00735 using Teuchos::ArrayRCP; 00736 using Teuchos::arcp_const_cast; 00737 const size_t numVecs = getNumVectors(), 00738 myLen = getLocalLength(); 00739 #ifdef HAVE_TPETRA_DEBUG 00740 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error, 00741 "Tpetra::MultiVector::scale(alpha,A): MultiVectors do not have compatible Maps:" << std::endl 00742 << "this->getMap(): " << std::endl << *this->getMap() 00743 << "A.getMap(): " << std::endl << *A.getMap() << std::endl); 00744 #else 00745 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error, 00746 "Tpetra::MultiVector::scale(alpha,A): MultiVectors do not have the same local length."); 00747 #endif 00748 TEST_FOR_EXCEPTION(A.getNumVectors() != numVecs, std::runtime_error, 00749 "Tpetra::MultiVector::scale(alpha,A): MultiVectors must have the same number of vectors."); 00750 if (isConstantStride() && A.isConstantStride()) { 00751 // set me == alpha*A 00752 MVT::Scale(lclMV_,alpha,(const KMV&)A.lclMV_); 00753 } 00754 else { 00755 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)); 00756 ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_), 00757 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)); 00758 for (size_t j=0; j < numVecs; ++j) { 00759 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00760 avj = A.getSubArrayRCP(avptr,j); 00761 MVT::initializeValues(a,myLen, 1, avj, myLen); 00762 MVT::initializeValues(v,myLen, 1, vj, myLen); 00763 MVT::Scale(v,alpha,(const KMV &)a); 00764 } 00765 } 00766 } 00767 00768 00769 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00770 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::reciprocal(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) { 00771 using Teuchos::ArrayRCP; 00772 using Teuchos::arcp_const_cast; 00773 #ifdef HAVE_TPETRA_DEBUG 00774 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error, 00775 "Tpetra::MultiVector::reciprocal(): MultiVectors do not have compatible Maps:" << std::endl 00776 << "this->getMap(): " << std::endl << *this->getMap() 00777 << "A.getMap(): " << std::endl << *A.getMap() << std::endl); 00778 #else 00779 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error, 00780 "Tpetra::MultiVector::reciprocal(): MultiVectors do not have the same local length."); 00781 #endif 00782 TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error, 00783 "Tpetra::MultiVector::reciprocal(): MultiVectors must have the same number of vectors."); 00784 const size_t numVecs = getNumVectors(); 00785 const size_t myLen = getLocalLength(); 00786 try { 00787 if (isConstantStride() && A.isConstantStride()) { 00788 MVT::Recip(lclMV_,(const KMV&)A.lclMV_); 00789 } 00790 else { 00791 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)); 00792 ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_), 00793 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)); 00794 for (size_t j=0; j < numVecs; ++j) { 00795 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00796 avj = A.getSubArrayRCP(avptr,j); 00797 MVT::initializeValues(a,myLen, 1, avj, myLen); 00798 MVT::initializeValues(v,myLen, 1, vj, myLen); 00799 MVT::Recip(v,(const KMV &)a); 00800 } 00801 } 00802 } 00803 catch (std::runtime_error &e) { 00804 TEST_FOR_EXCEPTION(true,std::runtime_error, 00805 "Tpetra::MultiVector::reciprocal(): caught exception from Kokkos:" << std::endl 00806 << e.what() << std::endl); 00807 } 00808 } 00809 00810 00811 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00812 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::abs(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) { 00813 using Teuchos::ArrayRCP; 00814 using Teuchos::arcp_const_cast; 00815 #ifdef HAVE_TPETRA_DEBUG 00816 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error, 00817 "Tpetra::MultiVector::abs(): MultiVectors do not have compatible Maps:" << std::endl 00818 << "this->getMap(): " << std::endl << *this->getMap() 00819 << "A.getMap(): " << std::endl << *A.getMap() << std::endl); 00820 #else 00821 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error, 00822 "Tpetra::MultiVector::abs(): MultiVectors do not have the same local length."); 00823 #endif 00824 TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error, 00825 "Tpetra::MultiVector::abs(): MultiVectors must have the same number of vectors."); 00826 const size_t myLen = getLocalLength(); 00827 const size_t numVecs = getNumVectors(); 00828 if (isConstantStride() && A.isConstantStride()) { 00829 MVT::Abs(lclMV_,(const KMV&)A.lclMV_); 00830 } 00831 else { 00832 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)); 00833 ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_), 00834 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)); 00835 for (size_t j=0; j < numVecs; ++j) { 00836 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00837 avj = A.getSubArrayRCP(avptr,j); 00838 MVT::initializeValues(a,myLen, 1, avj, myLen); 00839 MVT::initializeValues(v,myLen, 1, vj, myLen); 00840 MVT::Abs(v,(const KMV &)a); 00841 } 00842 } 00843 } 00844 00845 00846 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00847 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::update( 00848 const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 00849 const Scalar &beta) { 00850 // this = beta*this + alpha*A 00851 // must support case where &this == &A 00852 // can't short circuit on alpha==0.0 or beta==0.0, because 0.0*NaN != 0.0 00853 using Teuchos::ArrayRCP; 00854 using Teuchos::arcp_const_cast; 00855 #ifdef HAVE_TPETRA_DEBUG 00856 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error, 00857 "Tpetra::MultiVector::update(): MultiVectors do not have compatible Maps:" << std::endl 00858 << "this->getMap(): " << std::endl << this->getMap() 00859 << "A.getMap(): " << std::endl << *A.getMap() << std::endl); 00860 #else 00861 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error, 00862 "Tpetra::MultiVector::update(): MultiVectors do not have the same local length."); 00863 #endif 00864 const size_t myLen = getLocalLength(); 00865 const size_t numVecs = getNumVectors(); 00866 TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error, 00867 "Tpetra::MultiVector::update(): MultiVectors must have the same number of vectors."); 00868 if (isConstantStride() && A.isConstantStride()) { 00869 MVT::GESUM(lclMV_,alpha,(const KMV&)A.lclMV_,beta); 00870 } 00871 else { 00872 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)); 00873 ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_), 00874 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)); 00875 for (size_t j=0; j < numVecs; ++j) { 00876 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00877 avj = A.getSubArrayRCP(avptr,j); 00878 MVT::initializeValues(a,myLen, 1, avj, myLen); 00879 MVT::initializeValues(v,myLen, 1, vj, myLen); 00880 MVT::GESUM(v,alpha,(const KMV &)a,beta); 00881 } 00882 } 00883 } 00884 00885 00886 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00887 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::update( 00888 const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 00889 const Scalar &beta, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, 00890 const Scalar &gamma) { 00891 using Teuchos::ArrayRCP; 00892 using Teuchos::arcp_const_cast; 00893 // this = alpha*A + beta*B + gamma*this 00894 // must support case where &this == &A or &this == &B 00895 // can't short circuit on alpha==0.0 or beta==0.0 or gamma==0.0, because 0.0*NaN != 0.0 00896 #ifdef HAVE_TPETRA_DEBUG 00897 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()) || !this->getMap()->isCompatible(*B.getMap()), 00898 std::runtime_error, 00899 "Tpetra::MultiVector::update(): MultiVectors do not have compatible Maps:" << std::endl 00900 << "this->getMap(): " << std::endl << *this->getMap() 00901 << "A.getMap(): " << std::endl << *A.getMap() << std::endl 00902 << "B.getMap(): " << std::endl << *B.getMap() << std::endl); 00903 #else 00904 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength() || getLocalLength() != B.getLocalLength(), std::runtime_error, 00905 "Tpetra::MultiVector::update(): MultiVectors do not have the same local length."); 00906 #endif 00907 TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors() || B.getNumVectors() != this->getNumVectors(), std::runtime_error, 00908 "Tpetra::MultiVector::update(): MultiVectors must have the same number of vectors."); 00909 const size_t myLen = getLocalLength(); 00910 const size_t numVecs = getNumVectors(); 00911 if (isConstantStride() && A.isConstantStride() && B.isConstantStride()) { 00912 MVT::GESUM(lclMV_,alpha,(const KMV&)A.lclMV_,beta,(const KMV&)B.lclMV_,gamma); 00913 } 00914 else { 00915 KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)), b(MVT::getNode(lclMV_)); 00916 ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_), 00917 avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)), 00918 bvptr = arcp_const_cast<Scalar>(MVT::getValues(B.lclMV_)); 00919 for (size_t j=0; j < numVecs; ++j) { 00920 ArrayRCP<Scalar> vj = getSubArrayRCP( vptr,j), 00921 avj = A.getSubArrayRCP(avptr,j), 00922 bvj = B.getSubArrayRCP(bvptr,j); 00923 MVT::initializeValues(b,myLen, 1, bvj, myLen); 00924 MVT::initializeValues(a,myLen, 1, avj, myLen); 00925 MVT::initializeValues(v,myLen, 1, vj, myLen); 00926 MVT::GESUM(v,alpha,(const KMV&)a,beta,(const KMV&)b,gamma); 00927 } 00928 } 00929 } 00930 00931 00932 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00933 Teuchos::ArrayRCP<const Scalar> 00934 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getData(size_t j) const { 00935 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00936 KOKKOS_NODE_TRACE("MultiVector::getData()") 00937 return node->template viewBuffer<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_),j) ); 00938 } 00939 00940 00941 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00942 Teuchos::ArrayRCP<Scalar> 00943 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDataNonConst(size_t j) { 00944 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00945 KOKKOS_NODE_TRACE("MultiVector::getDataNonConst()") 00946 return node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite, getLocalLength(), getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j) ); 00947 } 00948 00949 00950 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00951 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & 00952 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::operator=(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source) { 00953 // Check for special case of this=Source, in which case we do nothing 00954 if (this != &source) { 00955 #ifdef HAVE_TPETRA_DEBUG 00956 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*source.getMap()), std::runtime_error, 00957 "Tpetra::MultiVector::operator=(): MultiVectors do not have compatible Maps:" << std::endl 00958 << "this->getMap(): " << std::endl << *this->getMap() 00959 << "source.getMap(): " << std::endl << *source.getMap() << std::endl); 00960 #else 00961 TEST_FOR_EXCEPTION( getLocalLength() != source.getLocalLength(), std::runtime_error, 00962 "Tpetra::MultiVector::operator=(): MultiVectors do not have the same local length."); 00963 #endif 00964 TEST_FOR_EXCEPTION(source.getNumVectors() != getNumVectors(), std::runtime_error, 00965 "Tpetra::MultiVector::operator=(): MultiVectors must have the same number of vectors."); 00966 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00967 const size_t numVecs = getNumVectors(); 00968 if (isConstantStride() && source.isConstantStride() && getLocalLength()==getStride() && source.getLocalLength()==source.getStride()) { 00969 // we're both packed, we can copy in one call 00970 KOKKOS_NODE_TRACE("MultiVector::operator=()") 00971 node->template copyBuffers<Scalar>(getLocalLength()*numVecs, MVT::getValues(source.lclMV_), MVT::getValuesNonConst(lclMV_) ); 00972 } 00973 else { 00974 for (size_t j=0; j < numVecs; ++j) { 00975 KOKKOS_NODE_TRACE("MultiVector::operator=()") 00976 node->template copyBuffers<Scalar>(getLocalLength(), source.getSubArrayRCP(MVT::getValues(source.lclMV_),j), 00977 getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j) ); 00978 } 00979 } 00980 } 00981 return(*this); 00982 } 00983 00984 00985 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 00986 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 00987 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subCopy(const Teuchos::ArrayView<const size_t> &cols) const { 00988 TEST_FOR_EXCEPTION(cols.size() < 1, std::runtime_error, 00989 "Tpetra::MultiVector::subCopy(cols): cols must contain at least one column."); 00990 size_t numCopyVecs = cols.size(); 00991 const bool zeroData = false; 00992 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 00993 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mv; 00994 // mv is allocated with constant stride 00995 mv = Teuchos::rcp( new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),numCopyVecs,zeroData) ); 00996 // copy data from *this into mv 00997 for (size_t j=0; j<numCopyVecs; ++j) { 00998 KOKKOS_NODE_TRACE("MultiVector::subCopy()") 00999 node->template copyBuffers<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_), cols[j]), 01000 MVT::getValuesNonConst(mv->lclMV_,j) ); 01001 } 01002 return mv; 01003 } 01004 01005 01006 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01007 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01008 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subCopy(const Teuchos::Range1D &colRng) const { 01009 TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error, 01010 "Tpetra::MultiVector::subCopy(Range1D): range must include at least one vector."); 01011 size_t numCopyVecs = colRng.size(); 01012 const bool zeroData = false; 01013 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01014 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mv; 01015 // mv is allocated with constant stride 01016 mv = Teuchos::rcp( new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),numCopyVecs,zeroData) ); 01017 // copy data from *this into mv 01018 for (size_t js=colRng.lbound(), jd=0; jd<numCopyVecs; ++jd, ++js) { 01019 KOKKOS_NODE_TRACE("MultiVector::subCopy()") 01020 node->template copyBuffers<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_), js), 01021 MVT::getValuesNonConst(mv->lclMV_,jd) ); 01022 } 01023 return mv; 01024 } 01025 01026 01027 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01028 Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01029 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::offsetView(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap, size_t offset) const 01030 { 01031 typedef const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> CMV; 01032 TEST_FOR_EXCEPTION( subMap->getNodeNumElements() + offset > this->getLocalLength(), std::runtime_error, 01033 "Tpetra::MultiVector::offsetView(subMap,offset): sizes are not sane.\noffset == " << offset << "\nsubMap: " << subMap->description() << "\nthis->rowMap: " << this->getMap()->description()); 01034 const size_t numVecs = this->getNumVectors(), 01035 myStride = MVT::getStride(lclMV_), 01036 newLen = subMap->getNodeNumElements(); 01037 using Teuchos::ArrayRCP; 01038 // this is const, so the lclMV_ is const, so that we can only get const buffers 01039 // we will cast away the const; this is okay, because 01040 // a) the constructor doesn't modify the data, and 01041 // b) we are encapsulating in a const MV before returning 01042 ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_); 01043 ArrayRCP<Scalar> ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf); 01044 Teuchos::RCP<CMV> constViewMV; 01045 if (isConstantStride()) { 01046 // view goes from first entry of first vector to last entry of last vector 01047 ArrayRCP<Scalar> subdata = ncbuf.persistingView( offset, myStride * (numVecs-1) + newLen ); 01048 constViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,numVecs) ); 01049 } 01050 else { 01051 // use same which index, but with an offset start pointer 01052 size_t maxSubVecIndex = *std::max_element(whichVectors_.begin(), whichVectors_.end()); 01053 ArrayRCP<Scalar> subdata = ncbuf.persistingView( offset, myStride * maxSubVecIndex + newLen ); 01054 constViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,whichVectors_) ); 01055 } 01056 return constViewMV; 01057 } 01058 01059 01060 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01061 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01062 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::offsetViewNonConst(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap, size_t offset) 01063 { 01064 typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV; 01065 TEST_FOR_EXCEPTION( subMap->getNodeNumElements() + offset > this->getLocalLength(), std::runtime_error, 01066 "Tpetra::MultiVector::offsetView(subMap,offset): sizes are not sane.\noffset == " << offset << "\nsubMap: " << subMap->description() << "\nthis->rowMap: " << this->getMap()->description()); 01067 const size_t numVecs = this->getNumVectors(), 01068 myStride = MVT::getStride(lclMV_), 01069 newLen = subMap->getNodeNumElements(); 01070 using Teuchos::ArrayRCP; 01071 // this is const, so the lclMV_ is const, so that we can only get const buffers 01072 // we will cast away the const; this is okay, because 01073 // a) the constructor doesn't modify the data, and 01074 // b) we are encapsulating in a const MV before returning 01075 ArrayRCP<Scalar> buf = MVT::getValuesNonConst(lclMV_); 01076 Teuchos::RCP<MV> subViewMV; 01077 if (isConstantStride()) { 01078 // view goes from first entry of first vector to last entry of last vector 01079 ArrayRCP<Scalar> subdata = buf.persistingView( offset, myStride * (numVecs-1) + newLen ); 01080 subViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,numVecs) ); 01081 } 01082 else { 01083 // use same which index, but with an offset start pointer 01084 size_t maxSubVecIndex = *std::max_element(whichVectors_.begin(), whichVectors_.end()); 01085 ArrayRCP<Scalar> subdata = buf.persistingView( offset, myStride * maxSubVecIndex + newLen ); 01086 subViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,whichVectors_) ); 01087 } 01088 return subViewMV; 01089 } 01090 01091 01092 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01093 Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01094 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subView(const Teuchos::ArrayView<const size_t> &cols) const { 01095 using Teuchos::ArrayRCP; 01096 typedef const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> CMV; 01097 TEST_FOR_EXCEPTION(cols.size() == 0, std::runtime_error, 01098 "Tpetra::MultiVector::subView(ArrayView): range must include at least one vector."); 01099 // this is const, so the lclMV_ is const, so that we can only get const buffers 01100 // we will cast away the const; this is okay, because 01101 // a) the constructor doesn't modify the data, and 01102 // b) we are encapsulating in a const MV before returning 01103 const size_t myStride = MVT::getStride(lclMV_), 01104 myLen = MVT::getNumRows(lclMV_), 01105 numViewCols = cols.size(); 01106 // use the smallest view possible of the buffer: from the first element of the minInd vector to the last element of the maxInd vector 01107 // this minimizes overlap between views, and keeps view of the minimum amount necessary, in order to allow node to achieve maximum efficiency. 01108 // adjust the indices appropriately; shift so that smallest index is 0 01109 ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_); 01110 ArrayRCP<Scalar> ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf); 01111 Teuchos::Array<size_t> newCols(numViewCols); 01112 size_t minInd = Teuchos::OrdinalTraits<size_t>::max(), 01113 maxInd = Teuchos::OrdinalTraits<size_t>::zero(); 01114 if (isConstantStride()) { 01115 for (size_t j=0; j < numViewCols; ++j) { 01116 newCols[j] = cols[j]; 01117 if (newCols[j] < minInd) minInd = newCols[j]; 01118 if (maxInd < newCols[j]) maxInd = newCols[j]; 01119 } 01120 } 01121 else { 01122 for (size_t j=0; j < numViewCols; ++j) { 01123 newCols[j] = whichVectors_[cols[j]]; 01124 if (newCols[j] < minInd) minInd = newCols[j]; 01125 if (maxInd < newCols[j]) maxInd = newCols[j]; 01126 } 01127 } 01128 ArrayRCP<Scalar> minbuf = ncbuf.persistingView(minInd * myStride, myStride * (maxInd - minInd) + myLen); 01129 for (size_t j=0; j < numViewCols; ++j) { 01130 newCols[j] -= minInd; 01131 } 01132 Teuchos::RCP<CMV> constViewMV; 01133 constViewMV = Teuchos::rcp<CMV>( 01134 new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), minbuf, myStride, newCols()) 01135 ); 01136 return constViewMV; 01137 } 01138 01139 01140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01141 Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01142 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subView(const Teuchos::Range1D &colRng) const { 01143 TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error, 01144 "Tpetra::MultiVector::subView(Range1D): range must include at least one vector."); 01145 size_t numViewVecs = colRng.size(); 01146 using Teuchos::ArrayRCP; 01147 // this is const, so the lclMV_ is const, so that we can only get const buffers 01148 // we will cast away the const; this is okay, because 01149 // a) the constructor doesn't modify the data, and 01150 // b) we are encapsulating in a const MV before returning 01151 ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_); 01152 ArrayRCP<Scalar> ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf); 01153 // resulting MultiVector is constant stride only if *this is 01154 if (isConstantStride()) { 01155 // view goes from first entry of first vector to last entry of last vector 01156 ArrayRCP<Scalar> subdata = ncbuf.persistingView( MVT::getStride(lclMV_) * colRng.lbound(), 01157 MVT::getStride(lclMV_) * (numViewVecs-1) + getLocalLength() ); 01158 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01159 subdata,MVT::getStride(lclMV_),numViewVecs) ); 01160 } 01161 // otherwise, use a subset of this whichVectors_ to construct new multivector 01162 Teuchos::Array<size_t> whchvecs( whichVectors_.begin()+colRng.lbound(), whichVectors_.begin()+colRng.ubound()+1 ); 01163 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01164 ncbuf,MVT::getStride(lclMV_),whchvecs) ); 01165 } 01166 01167 01168 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01169 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01170 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subViewNonConst(const Teuchos::ArrayView<const size_t> &cols) { 01171 TEST_FOR_EXCEPTION(cols.size() == 0, std::runtime_error, 01172 "Tpetra::MultiVector::subViewNonConst(ArrayView): range must include at least one vector."); 01173 if (isConstantStride()) { 01174 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01175 MVT::getValuesNonConst(lclMV_),MVT::getStride(lclMV_), 01176 cols) ); 01177 } 01178 // else, lookup current whichVectors_ using cols 01179 Teuchos::Array<size_t> newcols(cols.size()); 01180 for (size_t j=0; j < Teuchos::as<size_t>(cols.size()); ++j) { 01181 newcols[j] = whichVectors_[cols[j]]; 01182 } 01183 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01184 MVT::getValuesNonConst(lclMV_),MVT::getStride(lclMV_), 01185 newcols()) ); 01186 } 01187 01188 01189 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01190 Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01191 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subViewNonConst(const Teuchos::Range1D &colRng) { 01192 TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error, 01193 "Tpetra::MultiVector::subViewNonConst(Range1D): range must include at least one vector."); 01194 size_t numViewVecs = colRng.size(); 01195 using Teuchos::ArrayRCP; 01196 // resulting MultiVector is constant stride only if *this is 01197 if (isConstantStride()) { 01198 // view goes from first entry of first vector to last entry of last vector 01199 const size_t stride = MVT::getStride(lclMV_); 01200 ArrayRCP<Scalar> data = MVT::getValuesNonConst(lclMV_); 01201 ArrayRCP<Scalar> subdata = data.persistingView( stride * colRng.lbound(), 01202 stride * (numViewVecs-1) + getLocalLength() ); 01203 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01204 subdata,stride,numViewVecs) ); 01205 } 01206 // otherwise, use a subset of this whichVectors_ to construct new multivector 01207 Teuchos::Array<size_t> whchvecs( whichVectors_.begin()+colRng.lbound(), whichVectors_.begin()+colRng.ubound()+1 ); 01208 const size_t stride = MVT::getStride(lclMV_); 01209 ArrayRCP<Scalar> data = MVT::getValuesNonConst(lclMV_); 01210 return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), 01211 data,stride,whchvecs) ); 01212 } 01213 01214 01215 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01216 Teuchos::RCP<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01217 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getVector(size_t j) const { 01218 #ifdef HAVE_TPETRA_DEBUG 01219 TEST_FOR_EXCEPTION(j < 0 || j >= this->getNumVectors(), std::runtime_error, 01220 "Tpetra::MultiVector::getVector(j): index j (== " << j << ") exceeds valid column range for this multivector."); 01221 #endif 01222 // this is const, so lclMV_ is const, so we get const buff 01223 // it is safe to cast away the const because we will wrap it in a const Vector below 01224 Teuchos::ArrayRCP<Scalar> ncbuff; 01225 if (getLocalLength() > 0) { 01226 Teuchos::ArrayRCP<const Scalar> cbuff = getSubArrayRCP(MVT::getValues(lclMV_),j); 01227 ncbuff = Teuchos::arcp_const_cast<Scalar>(cbuff); 01228 } 01229 return Teuchos::rcp<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >( 01230 new Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),ncbuff) 01231 ); 01232 } 01233 01234 01235 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01236 Teuchos::RCP<Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01237 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getVectorNonConst(size_t j) { 01238 #ifdef HAVE_TPETRA_DEBUG 01239 TEST_FOR_EXCEPTION(j < 0 || j >= this->getNumVectors(), std::runtime_error, 01240 "Tpetra::MultiVector::getVectorNonConst(j): index j (== " << j << ") exceeds valid column range for this multivector."); 01241 #endif 01242 Teuchos::ArrayRCP<Scalar> ncbuff; 01243 if (getLocalLength() > 0) { 01244 ncbuff = getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j); 01245 } 01246 return Teuchos::rcp(new Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),ncbuff)); 01247 } 01248 01249 01250 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01251 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dCopy(Teuchos::ArrayView<Scalar> A, size_t LDA) const 01252 { 01253 using Teuchos::ArrayRCP; 01254 TEST_FOR_EXCEPTION(LDA < getLocalLength(), std::runtime_error, 01255 "Tpetra::MultiVector::get1dCopy(A,LDA): specified stride is not large enough for local vector length."); 01256 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(A.size()) < LDA*(getNumVectors()-1)+getLocalLength(), std::runtime_error, 01257 "Tpetra::MultiVector::get1dCopy(A,LDA): specified stride/storage is not large enough for the number of vectors."); 01258 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01259 const size_t myStride = MVT::getStride(lclMV_), 01260 numCols = getNumVectors(), 01261 myLen = getLocalLength(); 01262 if (myLen > 0) { 01263 ArrayRCP<const Scalar> mydata = MVT::getValues(lclMV_); 01264 KOKKOS_NODE_TRACE("MultiVector::getVectorNonConst()") 01265 ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,mydata); 01266 typename Teuchos::ArrayView<Scalar>::iterator Aptr = A.begin(); 01267 for (size_t j=0; j<numCols; j++) { 01268 ArrayRCP<const Scalar> myviewj = getSubArrayRCP(myview,j); 01269 std::copy(myviewj,myviewj+myLen,Aptr); 01270 Aptr += LDA; 01271 } 01272 myview = Teuchos::null; 01273 mydata = Teuchos::null; 01274 } 01275 } 01276 01277 01278 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01279 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dCopy(Teuchos::ArrayView<const Teuchos::ArrayView<Scalar> > ArrayOfPtrs) const 01280 { 01281 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs.size()) != getNumVectors(), std::runtime_error, 01282 "Tpetra::MultiVector::get2dCopy(ArrayOfPtrs): Array of pointers must contain as many pointers as the MultiVector has rows."); 01283 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01284 const size_t numCols = getNumVectors(), 01285 myLen = getLocalLength(); 01286 if (myLen > 0) { 01287 Teuchos::ArrayRCP<const Scalar> mybuff = MVT::getValues(lclMV_); 01288 KOKKOS_NODE_TRACE("MultiVector::get2dCopy()") 01289 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(mybuff.size(), mybuff); 01290 for (size_t j=0; j<numCols; ++j) { 01291 #ifdef HAVE_TPETRA_DEBUG 01292 TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != getLocalLength(), std::runtime_error, 01293 "Tpetra::MultiVector::get2dCopy(ArrayOfPtrs): The ArrayView provided in ArrayOfPtrs[" << j << "] was not large enough to contain the local entries."); 01294 #endif 01295 Teuchos::ArrayRCP<const Scalar> myviewj = getSubArrayRCP(myview,j); 01296 std::copy(myviewj,myviewj+myLen,ArrayOfPtrs[j].begin()); 01297 } 01298 myview = Teuchos::null; 01299 mybuff = Teuchos::null; 01300 } 01301 } 01302 01303 01304 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01305 Teuchos::ArrayRCP<const Scalar> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dView() const 01306 { 01307 TEST_FOR_EXCEPTION(!isConstantStride(), std::runtime_error, 01308 "Tpetra::MultiVector::get1dView() requires that this MultiVector have constant stride."); 01309 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01310 KOKKOS_NODE_TRACE("MultiVector::get1dView()") 01311 return node->template viewBuffer<Scalar>( getStride()*(getNumVectors()-1)+getLocalLength(), MVT::getValues(lclMV_) ); 01312 } 01313 01314 01315 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01316 Teuchos::ArrayRCP<Scalar> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dViewNonConst() 01317 { 01318 TEST_FOR_EXCEPTION(!isConstantStride(), std::runtime_error, 01319 "Tpetra::MultiVector::get1dViewNonConst(): requires that this MultiVector have constant stride."); 01320 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01321 KOKKOS_NODE_TRACE("MultiVector::get1dViewNonConst()") 01322 return node->template viewBufferNonConst<Scalar>( Kokkos::ReadWrite, getStride()*(getNumVectors()-1)+getLocalLength(), MVT::getValuesNonConst(lclMV_) ); 01323 } 01324 01325 01326 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01327 Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > 01328 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dViewNonConst() 01329 { 01330 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01331 Teuchos::ArrayRCP< Teuchos::ArrayRCP<Scalar> > views = Teuchos::arcp<Teuchos::ArrayRCP<Scalar> >(getNumVectors()); 01332 if (isConstantStride()) { 01333 const size_t myStride = MVT::getStride(lclMV_), 01334 numCols = getNumVectors(), 01335 myLen = getLocalLength(); 01336 if (myLen > 0) { 01337 KOKKOS_NODE_TRACE("MultiVector::get2dViewNonConst()") 01338 Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite,myStride*(numCols-1)+myLen,MVT::getValuesNonConst(lclMV_)); 01339 for (size_t j=0; j<numCols; ++j) { 01340 views[j] = myview.persistingView(0,myLen); 01341 myview += myStride; 01342 } 01343 } 01344 } 01345 else { 01346 const size_t myStride = MVT::getStride(lclMV_), 01347 numCols = MVT::getNumCols(lclMV_), 01348 myCols = getNumVectors(), 01349 myLen = MVT::getNumRows(lclMV_); 01350 if (myLen > 0) { 01351 KOKKOS_NODE_TRACE("MultiVector::get2dViewNonConst()") 01352 Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite,myStride*(numCols-1)+myLen,MVT::getValuesNonConst(lclMV_)); 01353 for (size_t j=0; j<myCols; ++j) { 01354 views[j] = myview.persistingView(whichVectors_[j]*myStride,myLen); 01355 } 01356 } 01357 } 01358 return views; 01359 } 01360 01361 01362 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01363 Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > 01364 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dView() const 01365 { 01366 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01367 Teuchos::ArrayRCP< Teuchos::ArrayRCP<const Scalar> > views = Teuchos::arcp<Teuchos::ArrayRCP<const Scalar> >(getNumVectors()); 01368 if (isConstantStride()) { 01369 const size_t myStride = MVT::getStride(lclMV_), 01370 numCols = getNumVectors(), 01371 myLen = getLocalLength(); 01372 if (myLen > 0) { 01373 KOKKOS_NODE_TRACE("MultiVector::get2dView()") 01374 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,MVT::getValues(lclMV_)); 01375 for (size_t j=0; j<numCols; ++j) { 01376 views[j] = myview.persistingView(0,myLen); 01377 myview += myStride; 01378 } 01379 } 01380 } 01381 else { 01382 const size_t myStride = MVT::getStride(lclMV_), 01383 numCols = MVT::getNumCols(lclMV_), 01384 myCols = getNumVectors(), 01385 myLen = MVT::getNumRows(lclMV_); 01386 if (myLen > 0) { 01387 KOKKOS_NODE_TRACE("MultiVector::get2dView()") 01388 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,MVT::getValues(lclMV_)); 01389 for (size_t j=0; j<myCols; ++j) { 01390 views[j] = myview.persistingView(whichVectors_[j]*myStride,myLen); 01391 } 01392 } 01393 } 01394 return views; 01395 } 01396 01397 01398 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01399 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::multiply( 01400 Teuchos::ETransp transA, Teuchos::ETransp transB, 01401 const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, 01402 const Scalar &beta) { 01403 using Teuchos::NO_TRANS; // enums 01404 using Teuchos::TRANS; 01405 using Teuchos::CONJ_TRANS; 01406 using Teuchos::null; 01407 using Teuchos::ScalarTraits; // traits 01408 using Teuchos::as; 01409 using Teuchos::RCP; // data structures 01410 using Teuchos::Array; 01411 using Teuchos::ArrayView; 01412 using Teuchos::rcp; // initializers for data structures 01413 01414 // This routine performs a variety of matrix-matrix multiply operations, interpreting 01415 // the MultiVector (this-aka C , A and B) as 2D matrices. Variations are due to 01416 // the fact that A, B and C can be local replicated or global distributed 01417 // MultiVectors and that we may or may not operate with the transpose of 01418 // A and B. Possible cases are: 01419 // Num 01420 // OPERATIONS cases Notes 01421 // 1) C(local) = A^X(local) * B^X(local) 4 (X=Trans or Not, No comm needed) 01422 // 2) C(local) = A^T(distr) * B (distr) 1 (2D dot product, replicate C) 01423 // 3) C(distr) = A (distr) * B^X(local) 2 (2D vector update, no comm needed) 01424 // 01425 // The following operations are not meaningful for 1D distributions: 01426 // 01427 // u1) C(local) = A^T(distr) * B^T(distr) 1 01428 // u2) C(local) = A (distr) * B^X(distr) 2 01429 // u3) C(distr) = A^X(local) * B^X(local) 4 01430 // u4) C(distr) = A^X(local) * B^X(distr) 4 01431 // u5) C(distr) = A^T(distr) * B^X(local) 2 01432 // u6) C(local) = A^X(distr) * B^X(local) 4 01433 // u7) C(distr) = A^X(distr) * B^X(local) 4 01434 // u8) C(local) = A^X(local) * B^X(distr) 4 01435 // 01436 // Total of 32 case (2^5). 01437 01438 std::string errPrefix("Tpetra::MultiVector::multiply(transOpA,transOpB,alpha,A,B,beta): "); 01439 01440 TEST_FOR_EXCEPTION( ScalarTraits<Scalar>::isComplex && (transA == TRANS || transB == TRANS), std::invalid_argument, 01441 errPrefix << "non-conjugate transpose not supported for complex types."); 01442 transA = (transA == NO_TRANS ? NO_TRANS : CONJ_TRANS); 01443 transB = (transB == NO_TRANS ? NO_TRANS : CONJ_TRANS); 01444 01445 // Compute effective dimensions, w.r.t. transpose operations on 01446 size_t A_nrows = (transA==CONJ_TRANS) ? A.getNumVectors() : A.getLocalLength(); 01447 size_t A_ncols = (transA==CONJ_TRANS) ? A.getLocalLength() : A.getNumVectors(); 01448 size_t B_nrows = (transB==CONJ_TRANS) ? B.getNumVectors() : B.getLocalLength(); 01449 size_t B_ncols = (transB==CONJ_TRANS) ? B.getLocalLength() : B.getNumVectors(); 01450 01451 Scalar beta_local = beta; // local copy of beta; might be reassigned below 01452 01453 TEST_FOR_EXCEPTION( getLocalLength() != A_nrows || getNumVectors() != B_ncols || A_ncols != B_nrows, std::runtime_error, 01454 errPrefix << "dimension of *this, op(A) and op(B) must be consistent."); 01455 01456 bool A_is_local = !A.isDistributed(); 01457 bool B_is_local = !B.isDistributed(); 01458 bool C_is_local = !this->isDistributed(); 01459 bool Case1 = ( C_is_local && A_is_local && B_is_local); // Case 1: C(local) = A^X(local) * B^X(local) 01460 bool Case2 = ( C_is_local && !A_is_local && !B_is_local && transA==CONJ_TRANS && transB==NO_TRANS); // Case 2: C(local) = A^T(distr) * B (distr) 01461 bool Case3 = (!C_is_local && !A_is_local && B_is_local && transA==NO_TRANS ); // Case 3: C(distr) = A (distr) * B^X(local) 01462 01463 // Test that we are considering a meaningful cases 01464 TEST_FOR_EXCEPTION( !Case1 && !Case2 && !Case3, std::runtime_error, 01465 errPrefix << "multiplication of op(A) and op(B) into *this is not a supported use case."); 01466 01467 if (beta != ScalarTraits<Scalar>::zero() && Case2) 01468 { 01469 // if Case2, then C is local and contributions must be summed across all nodes 01470 // however, if beta != 0, then accumulate beta*C into the sum 01471 // when summing across all nodes, we only want to accumulate this once, so 01472 // set beta == 0 on all nodes except node 0 01473 int MyPID = this->getMap()->getComm()->getRank(); 01474 if (MyPID!=0) beta_local = ScalarTraits<Scalar>::zero(); 01475 } 01476 01477 // Check if A, B, C have constant stride, if not then make temp copy (strided) 01478 RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Atmp, Btmp; 01479 RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Ctmp; 01480 if (isConstantStride() == false) Ctmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(*this)); 01481 else Ctmp = rcp(this,false); 01482 01483 if (A.isConstantStride() == false) Atmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A)); 01484 else Atmp = rcp(&A,false); 01485 01486 if (B.isConstantStride() == false) Btmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(B)); 01487 else Btmp = rcp(&B,false); 01488 01489 #ifdef HAVE_TEUCHOS_DEBUG 01490 TEST_FOR_EXCEPTION(!Ctmp->isConstantStride() || !Btmp->isConstantStride() || !Atmp->isConstantStride(), std::logic_error, 01491 errPrefix << "failed making temporary strided copies of input multivectors."); 01492 #endif 01493 01494 KMV &C_mv = Ctmp->lclMV_; 01495 { 01496 // get local multivectors 01497 const KMV &A_mv = Atmp->lclMV_; 01498 const KMV &B_mv = Btmp->lclMV_; 01499 // do the multiply (GEMM) 01500 MVT::GEMM(C_mv,transA,transB,alpha,A_mv,B_mv,beta_local); 01501 } 01502 01503 // Dispose of (possibly) extra copies of A, B 01504 Atmp = null; 01505 Btmp = null; 01506 01507 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01508 // If *this was not strided, copy the data from the strided version and then delete it 01509 if (isConstantStride() == false) { 01510 // *this is not strided, we must put data from Ctmp into *this 01511 TEST_FOR_EXCEPT(&C_mv != &lclMV_); 01512 const size_t numVecs = MVT::getNumCols(lclMV_); 01513 for (size_t j=0; j < numVecs; ++j) { 01514 KOKKOS_NODE_TRACE("MultiVector::multiply()") 01515 node->template copyBuffers<Scalar>(getLocalLength(),MVT::getValues(C_mv,j),MVT::getValuesNonConst(lclMV_,whichVectors_[j])); 01516 } 01517 } 01518 01519 // If Case 2 then sum up *this and distribute it to all processors. 01520 if (Case2) { 01521 this->reduce(); 01522 } 01523 } 01524 01525 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01526 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::elementWiseMultiply(Scalar scalarAB, const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, Scalar scalarThis) 01527 { 01528 using Teuchos::ArrayRCP; 01529 using Teuchos::arcp_const_cast; 01530 #ifdef HAVE_TPETRA_DEBUG 01531 TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()) || !this->getMap()->isCompatible(*B.getMap()), std::runtime_error, 01532 "Tpetra::MultiVector::elementWiseMultiply(): MultiVectors do not have compatible Maps:" << std::endl 01533 << "this->getMap(): " << std::endl << *this->getMap() 01534 << "A.getMap(): " << std::endl << *A.getMap() << std::endl 01535 << "B.getMap(): " << std::endl << *B.getMap() << std::endl); 01536 #else 01537 TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength() || getLocalLength() != B.getLocalLength(), std::runtime_error, 01538 "Tpetra::MultiVector::reciprocal(): MultiVectors do not have the same local length."); 01539 #endif 01540 TEST_FOR_EXCEPTION(B.getNumVectors() != this->getNumVectors(), std::runtime_error, 01541 "Tpetra::MultiVector::reciprocal(): MultiVectors 'this' and B must have the same number of vectors."); 01542 TEST_FOR_EXCEPTION(A.getNumVectors() != 1, std::runtime_error, 01543 "Tpetra::MultiVector::reciprocal(): MultiVectors A must have just 1 vector."); 01544 try { 01545 MVT::ElemMult(lclMV_,scalarThis,scalarAB,(const KMV &)A.lclMV_,(const KMV &)B.lclMV_); 01546 } 01547 catch (std::runtime_error &e) { 01548 TEST_FOR_EXCEPTION(true,std::runtime_error, 01549 "Tpetra::MultiVector::elementWiseMultiply(): caught exception from Kokkos:" << std::endl 01550 << e.what() << std::endl); 01551 } 01552 } 01553 01554 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01555 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::reduce() { 01556 using Teuchos::ArrayRCP; 01557 using Teuchos::ArrayView; 01558 // this should be called only for "local" MultiVectors (!isDistributed()) 01559 TEST_FOR_EXCEPTION(this->isDistributed() == true, std::runtime_error, 01560 "Tpetra::MultiVector::reduce() should only be called for non-distributed MultiVectors."); 01561 Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm(); 01562 if (comm->getSize() == 1) return; 01563 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01564 // sum the data across all multivectors 01565 // need to have separate packed buffers for send and receive 01566 // if we're packed, we'll just set the receive buffer as our data, the send as a copy 01567 // if we're not packed, we'll use allocated buffers for both. 01568 ArrayView<Scalar> target; 01569 const size_t myStride = MVT::getStride(lclMV_), 01570 numCols = MVT::getNumCols(lclMV_), 01571 myLen = MVT::getNumRows(lclMV_); 01572 Teuchos::Array<Scalar> sourceBuffer(numCols*myLen), tmparr(0); 01573 bool packed = isConstantStride() && (myStride == myLen); 01574 KOKKOS_NODE_TRACE("MultiVector::reduce()") 01575 ArrayRCP<Scalar> bufView = node->template viewBufferNonConst<Scalar>( 01576 Kokkos::ReadWrite,myStride*(numCols-1)+myLen, 01577 MVT::getValuesNonConst(lclMV_) ); 01578 if (packed) { 01579 // copy data from view to sourceBuffer, reduce into view below 01580 target = bufView(0,myLen*numCols); 01581 std::copy(target.begin(),target.end(),sourceBuffer.begin()); 01582 } 01583 else { 01584 // copy data into sourceBuffer, reduce from sourceBuffer into tmparr below, copy back to view after that 01585 tmparr.resize(myLen*numCols); 01586 Scalar *sptr = sourceBuffer.getRawPtr(); 01587 ArrayRCP<const Scalar> vptr = bufView; 01588 for (size_t j=0; j<numCols; ++j) 01589 { 01590 std::copy(vptr,vptr+myLen,sptr); 01591 sptr += myLen; 01592 vptr += myStride; 01593 } 01594 target = tmparr(); 01595 } 01596 // reduce 01597 Teuchos::reduceAll<int,Scalar>(*comm,Teuchos::REDUCE_SUM,Teuchos::as<int>(numCols*myLen),sourceBuffer.getRawPtr(),target.getRawPtr()); 01598 if (!packed) { 01599 // copy tmparr back into view 01600 const Scalar *sptr = tmparr.getRawPtr(); 01601 ArrayRCP<Scalar> vptr = bufView; 01602 for (size_t j=0; j<numCols; ++j) 01603 { 01604 std::copy(sptr,sptr+myLen,vptr); 01605 sptr += myLen; 01606 vptr += myStride; 01607 } 01608 } 01609 bufView = Teuchos::null; 01610 } 01611 01612 01613 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01614 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(LocalOrdinal MyRow, size_t VectorIndex, const Scalar &ScalarValue) 01615 { 01616 #ifdef HAVE_TPETRA_DEBUG 01617 TEST_FOR_EXCEPTION(MyRow < this->getMap()->getMinLocalIndex() || MyRow > this->getMap()->getMaxLocalIndex(), std::runtime_error, 01618 "Tpetra::MultiVector::replaceLocalValue(): row index is invalid."); 01619 TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error, 01620 "Tpetra::MultiVector::replaceLocalValue(): vector index is invalid."); 01621 #endif 01622 getDataNonConst(VectorIndex)[MyRow] = ScalarValue; 01623 } 01624 01625 01626 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01627 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(LocalOrdinal MyRow, size_t VectorIndex, const Scalar &ScalarValue) 01628 { 01629 #ifdef HAVE_TPETRA_DEBUG 01630 TEST_FOR_EXCEPTION(MyRow < this->getMap()->getMinLocalIndex() || MyRow > this->getMap()->getMaxLocalIndex(), std::runtime_error, 01631 "Tpetra::MultiVector::sumIntoLocalValue(): row index is invalid."); 01632 TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error, 01633 "Tpetra::MultiVector::sumIntoLocalValue(): vector index is invalid."); 01634 #endif 01635 getDataNonConst(VectorIndex)[MyRow] += ScalarValue; 01636 } 01637 01638 01639 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01640 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(GlobalOrdinal GlobalRow, size_t VectorIndex, const Scalar &ScalarValue) 01641 { 01642 LocalOrdinal MyRow = this->getMap()->getLocalElement(GlobalRow); 01643 #ifdef HAVE_TPETRA_DEBUG 01644 TEST_FOR_EXCEPTION(MyRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error, 01645 "Tpetra::MultiVector::replaceGlobalValue(): row index is not present on this processor."); 01646 TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error, 01647 "Tpetra::MultiVector::replaceGlobalValue(): vector index is invalid."); 01648 #endif 01649 getDataNonConst(VectorIndex)[MyRow] = ScalarValue; 01650 } 01651 01652 01653 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01654 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(GlobalOrdinal GlobalRow, size_t VectorIndex, const Scalar &ScalarValue) 01655 { 01656 LocalOrdinal MyRow = this->getMap()->getLocalElement(GlobalRow); 01657 #ifdef HAVE_TEUCHOS_DEBUG 01658 TEST_FOR_EXCEPTION(MyRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error, 01659 "Tpetra::MultiVector::sumIntoGlobalValue(): row index is not present on this processor."); 01660 TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error, 01661 "Tpetra::MultiVector::sumIntoGlobalValue(): vector index is invalid."); 01662 #endif 01663 getDataNonConst(VectorIndex)[MyRow] += ScalarValue; 01664 } 01665 01666 01667 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01668 template <class T> 01669 Teuchos::ArrayRCP<T> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getSubArrayRCP(Teuchos::ArrayRCP<T> arr, size_t j) const { 01670 const size_t stride = MVT::getStride(lclMV_), 01671 myLen = getLocalLength(); 01672 Teuchos::ArrayRCP<T> ret; 01673 if (isConstantStride()) { 01674 ret = arr.persistingView(j*stride,myLen); 01675 } 01676 else { 01677 ret = arr.persistingView(whichVectors_[j]*stride,myLen); 01678 } 01679 return ret; 01680 } 01681 01682 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01683 const Kokkos::MultiVector<Scalar,Node> & 01684 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalMV() const { 01685 return lclMV_; 01686 } 01687 01688 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01689 Kokkos::MultiVector<Scalar,Node> & 01690 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalMVNonConst() { 01691 return lclMV_; 01692 } 01693 01694 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01695 std::string MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const 01696 { 01697 std::ostringstream oss; 01698 oss << Teuchos::Describable::description(); 01699 oss << "{length = "<<getGlobalLength() 01700 << ", getNumVectors = "<<getNumVectors() 01701 << ", isConstantStride() = "<<isConstantStride() 01702 << "}"; 01703 return oss.str(); 01704 } 01705 01706 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01707 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const { 01708 using std::endl; 01709 using std::setw; 01710 using Teuchos::VERB_DEFAULT; 01711 using Teuchos::VERB_NONE; 01712 using Teuchos::VERB_LOW; 01713 using Teuchos::VERB_MEDIUM; 01714 using Teuchos::VERB_HIGH; 01715 using Teuchos::VERB_EXTREME; 01716 Teuchos::EVerbosityLevel vl = verbLevel; 01717 if (vl == VERB_DEFAULT) vl = VERB_LOW; 01718 Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm(); 01719 const int myImageID = comm->getRank(), 01720 numImages = comm->getSize(); 01721 size_t width = 1; 01722 for (size_t dec=10; dec<getGlobalLength(); dec *= 10) { 01723 ++width; 01724 } 01725 Teuchos::OSTab tab(out); 01726 if (vl != VERB_NONE) { 01727 // VERB_LOW and higher prints description() 01728 if (myImageID == 0) out << this->description() << std::endl; 01729 for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) { 01730 if (myImageID == imageCtr) { 01731 if (vl != VERB_LOW) { 01732 // VERB_MEDIUM and higher prints getLocalLength() 01733 out << "node " << setw(width) << myImageID << ": local length=" << getLocalLength(); 01734 if (vl != VERB_MEDIUM) { 01735 // VERB_HIGH and higher prints isConstantStride() and getStride() 01736 if (isConstantStride()) out << ", constant stride=" << getStride() << endl; 01737 else out << ", non-constant stride" << endl; 01738 if (vl == VERB_EXTREME && getLocalLength() > 0) { 01739 Teuchos::RCP<Node> node = MVT::getNode(lclMV_); 01740 if (isConstantStride()) { 01741 KOKKOS_NODE_TRACE("MultiVector::describe()") 01742 Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>( 01743 getLocalLength()+getStride()*(getNumVectors()-1), 01744 MVT::getValues(lclMV_) ); 01745 // VERB_EXTREME prints values 01746 for (size_t i=0; i<getLocalLength(); ++i) { 01747 out << setw(width) << this->getMap()->getGlobalElement(i) << ": "; 01748 for (size_t j=0; j<getNumVectors(); ++j) { 01749 out << myview[j*getStride()] << " "; 01750 } 01751 ++myview; 01752 out << endl; 01753 } 01754 myview = Teuchos::null; 01755 } 01756 else { 01757 const size_t stride = MVT::getStride(lclMV_), 01758 rows = MVT::getNumRows(lclMV_), 01759 cols = MVT::getNumCols(lclMV_); 01760 KOKKOS_NODE_TRACE("MultiVector::describe()") 01761 Teuchos::ArrayRCP<const Scalar> myview = 01762 node->template viewBuffer<Scalar>( rows + stride * (cols - 1), MVT::getValues(lclMV_) ); 01763 // VERB_EXTREME prints values 01764 for (size_t i=0; i<getLocalLength(); ++i) { 01765 out << setw(width) << this->getMap()->getGlobalElement(i) << ": "; 01766 for (size_t j=0; j<getNumVectors(); ++j) { 01767 out << myview[whichVectors_[j]*stride + i] << " "; 01768 } 01769 out << endl; 01770 } 01771 myview = Teuchos::null; 01772 } 01773 } 01774 } 01775 else { 01776 out << endl; 01777 } 01778 } 01779 } 01780 } 01781 } 01782 } 01783 01784 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01785 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createViews() const { 01786 Teuchos::RCP<Node> node = this->getMap()->getNode(); 01787 if (cview_ == Teuchos::null && getLocalLength() > 0) { 01788 Teuchos::ArrayRCP<const Scalar> buff = MVT::getValues(lclMV_); 01789 KOKKOS_NODE_TRACE("MultiVector::createViews()") 01790 cview_ = node->template viewBuffer<Scalar>(buff.size(),buff); 01791 } 01792 } 01793 01794 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01795 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createViewsNonConst(Kokkos::ReadWriteOption rwo) { 01796 Teuchos::RCP<Node> node = this->getMap()->getNode(); 01797 if (ncview_ == Teuchos::null && getLocalLength() > 0) { 01798 Teuchos::ArrayRCP<Scalar> buff = MVT::getValuesNonConst(lclMV_); 01799 KOKKOS_NODE_TRACE("MultiVector::createViewsNonConst()") 01800 ncview_ = node->template viewBufferNonConst<Scalar>(rwo,buff.size(),buff); 01801 } 01802 } 01803 01804 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01805 void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::releaseViews() const { 01806 cview_ = Teuchos::null; 01807 ncview_ = Teuchos::null; 01808 } 01809 01810 01811 } // namespace Tpetra 01812 01813 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01814 Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01815 Tpetra::createMultiVector(const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > &map, size_t numVectors) { 01816 const bool DO_INIT_TO_ZERO = true; 01817 return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,numVectors,DO_INIT_TO_ZERO) ); 01818 } 01819 01820 // 01821 // Explicit instantiation macro 01822 // 01823 // Must be expanded from within the Tpetra namespace! 01824 // 01825 01826 #define TPETRA_MULTIVECTOR_INSTANT(SCALAR,LO,GO,NODE) \ 01827 \ 01828 template class MultiVector< SCALAR , LO , GO , NODE >; \ 01829 \ 01830 template Teuchos::RCP< MultiVector<SCALAR,LO,GO,NODE> > \ 01831 createMultiVector<SCALAR,LO,GO,NODE>(const Teuchos::RCP<const Map<LO,GO,NODE> > &map, size_t numVectors); 01832 01833 01834 #endif // TPETRA_MULTIVECTOR_DEF_HPP
1.7.4
