00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef TSFVECTORIMPL_HPP
00030 #define TSFVECTORIMPL_HPP
00031
00032
00033 #include "TSFVectorDecl.hpp"
00034 #include "Thyra_SUNDIALS_Ops.hpp"
00035 #include "TSFIndexableVector.hpp"
00036 #include "Thyra_DefaultProductVectorSpace.hpp"
00037 #include "Thyra_DefaultProductVector.hpp"
00038 #include "Thyra_DefaultSpmdVector.hpp"
00039 #include "SundancePrintable.hpp"
00040 #include "SundanceOut.hpp"
00041 #include "SundanceTabs.hpp"
00042
00043 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00044 #include "TSFVectorSpaceImpl.hpp"
00045 #include "TSFSequentialIteratorImpl.hpp"
00046 #endif
00047
00048
00049
00050
00051 namespace TSFExtended
00052 {
00053 using Sundance::Out;
00054 using Sundance::Tabs;
00055 using Sundance::Printable;
00056
00057
00058 template <class Scalar>
00059 void Vector<Scalar>::setBlock(int i, const Vector<Scalar>& v)
00060 {
00061 Thyra::DefaultProductVector<Scalar>* pv =
00062 dynamic_cast<Thyra::DefaultProductVector<Scalar>* >(this->ptr().get());
00063 TEST_FOR_EXCEPTION(pv == 0, std::runtime_error,
00064 "vector is not a product vector");
00065 Thyra::assign(pv->getNonconstVectorBlock(i).ptr(), *(v.ptr().ptr()));
00066 }
00067
00068
00069
00070
00071 template <class Scalar>
00072 Vector<Scalar> Vector<Scalar>::getBlock(int i) const
00073 {
00074 const Thyra::DefaultProductVector<Scalar>* pv =
00075 dynamic_cast <const Thyra::DefaultProductVector<Scalar>* >(this->ptr().get());
00076 if (pv==0)
00077 {
00078 TEST_FOR_EXCEPTION(i != 0, std::runtime_error,
00079 "Nonzero block index " << i << " into a vector that is not "
00080 "a product vector");
00081 return *this;
00082 }
00083 Teuchos::RCP<const Thyra::VectorBase<Scalar> > b = pv->getVectorBlock(i);
00084 Teuchos::RCP<Thyra::VectorBase<Scalar> > bb
00085 = rcp_const_cast<Thyra::VectorBase<Scalar> >(b);
00086 return bb;
00087 }
00088
00089
00090
00091 template <class Scalar>
00092 void Vector<Scalar>::print(std::ostream& os) const
00093 {
00094
00095 const Printable* p =
00096 dynamic_cast<const Printable* >(this->ptr().get());
00097 if (p != 0)
00098 {
00099 p->print(os);
00100 return;
00101 }
00102 const Thyra::ProductMultiVectorBase<Scalar>* pv =
00103 dynamic_cast <const Thyra::ProductMultiVectorBase<Scalar>* >(this->ptr().get());
00104 if (pv != 0)
00105 {
00106 os << "ProductVector[" << std::endl;
00107 for (int i=0; i<this->space().numBlocks(); i++)
00108 {
00109 os << "block=" << i << std::endl;
00110 this->getBlock(i).print(os);
00111 os << std::endl;
00112 }
00113 os << "]" << std::endl;
00114 return;
00115 }
00116 else
00117 {
00118 for (SequentialIterator<Scalar> i=this->space().begin(); i!=this->space().end(); i++)
00119 {
00120 os << i.globalIndex() << '\t' << (*this)[i] << std::endl;
00121 }
00122
00123 }
00124 }
00125
00126
00127
00128
00129
00130 template <class Scalar> inline
00131 const AccessibleVector<Scalar>* Vector<Scalar>::castToAccessible() const
00132 {
00133 const AccessibleVector<Scalar>* av
00134 = dynamic_cast<const AccessibleVector<Scalar>*>(this->ptr().get());
00135 TEST_FOR_EXCEPTION(av==0, std::runtime_error,
00136 "Attempted to cast non-accessible vector "
00137 << this->description() << " to an AccessibleVector");
00138 return av;
00139 }
00140
00141
00142 template <class Scalar> inline
00143 LoadableVector<Scalar>* Vector<Scalar>::castToLoadable()
00144 {
00145 LoadableVector<Scalar>* lv
00146 = dynamic_cast<LoadableVector<Scalar>*>(this->ptr().get());
00147 TEST_FOR_EXCEPTION(lv==0, std::runtime_error,
00148 "Attempted to cast non-loadable vector "
00149 << this->description() << " to a LoadableVector");
00150 return lv;
00151 }
00152
00153
00154
00155 template <class Scalar> inline
00156 const RawDataAccessibleVector<Scalar>* Vector<Scalar>::castToRawDataAccessible() const
00157 {
00158 const RawDataAccessibleVector<Scalar>* av
00159 = dynamic_cast<const RawDataAccessibleVector<Scalar>*>(this->ptr().get());
00160 TEST_FOR_EXCEPTION(av==0, std::runtime_error,
00161 "Attempted to cast non-accessible vector "
00162 << this->description()
00163 << " to an RawDataAccessibleVector");
00164 return av;
00165 }
00166
00167
00168 template <class Scalar> inline
00169 RawDataAccessibleVector<Scalar>* Vector<Scalar>::castToRawDataAccessible()
00170 {
00171 RawDataAccessibleVector<Scalar>* av
00172 = dynamic_cast<RawDataAccessibleVector<Scalar>*>(this->ptr().get());
00173 TEST_FOR_EXCEPTION(av==0, std::runtime_error,
00174 "Attempted to cast non-accessible vector "
00175 << this->description()
00176 << " to an RawDataAccessibleVector");
00177 return av;
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187 template <class Scalar> inline
00188 Vector<Scalar>& Vector<Scalar>::scale(const Scalar& alpha)
00189 {
00190 {
00191 TimeMonitor t(*opTimer());
00192 Thyra::Vt_S(this->ptr().ptr(), alpha);
00193 }
00194 return *this;
00195 }
00196
00197
00198
00199
00200
00201 template <class Scalar> inline
00202 Vector<Scalar>& Vector<Scalar>::update(const Scalar& alpha,
00203 const Vector<Scalar>& x)
00204 {
00205 Ptr<Thyra::VectorBase<Scalar> > p = this->ptr().ptr();
00206 p.assert_not_null();
00207 const Thyra::VectorBase<Scalar>* px = x.ptr().get();
00208 TEST_FOR_EXCEPT(px==0);
00209 {
00210 TimeMonitor t(*opTimer());
00211 Thyra::Vp_StV(p, alpha, *px);
00212 }
00213
00214 return *this;
00215 }
00216
00217
00218
00219
00220
00221
00222 template <class Scalar> inline
00223 Vector<Scalar>& Vector<Scalar>::acceptCopyOf(const Vector<Scalar>& x)
00224 {
00225 TimeMonitor t(*opTimer());
00226 if (this->ptr().get()==0)
00227 {
00228 Vector<Scalar> me = x.space().createMember();
00229 this->ptr() = me.ptr();
00230 }
00231
00232 Thyra::assign(this->ptr().ptr(), *(x.ptr()));
00233
00234 return *this;
00235 }
00236
00237 template <class Scalar> inline
00238 Vector<Scalar> Vector<Scalar>::copy() const
00239 {
00240 Vector<Scalar> rtn = space().createMember();
00241 {
00242 TimeMonitor t(*opTimer());
00243 rtn.acceptCopyOf(*this);
00244 }
00245 return rtn;
00246 }
00247
00248
00249
00250
00251 template <class Scalar> inline
00252 Vector<Scalar> Vector<Scalar>::dotStar(const Vector<Scalar>& other) const
00253 {
00254 Vector<Scalar> rtn = space().createMember();
00255 {
00256 TimeMonitor t(*opTimer());
00257
00258 Thyra::ele_wise_prod(1.0, *(this->ptr()), *(other.ptr()), rtn.ptr().ptr());
00259 }
00260 return rtn;
00261 }
00262
00263
00264
00265
00266 template <class Scalar> inline
00267 void Vector<Scalar>::dotStarInto(
00268 const Vector<Scalar>& a, const Vector<Scalar>& b) const
00269 {
00270 TimeMonitor t(*opTimer());
00271
00272 Thyra::ele_wise_prod(1.0, *(a.ptr()), *(b.ptr()), this->ptr().ptr());
00273 }
00274
00275
00276
00277
00278
00279
00280 template <class Scalar> inline
00281 Vector<Scalar> Vector<Scalar>::dotSlash(const Vector<Scalar>& other) const
00282 {
00283 Vector<Scalar> rtn = space().createMember();
00284 {
00285 TimeMonitor t(*opTimer());
00286 Thyra::ele_wise_divide(1.0, *(this->ptr)(), *(other.ptr()), rtn.ptr().ptr());
00287 }
00288 return rtn;
00289 }
00290
00291
00292
00293
00294
00295
00296
00297 template <class Scalar> inline
00298 Vector<Scalar> Vector<Scalar>::abs() const
00299 {
00300 Vector<Scalar> rtn = space().createMember();
00301 {
00302 TimeMonitor t(*opTimer());
00303 rtn.acceptCopyOf(*this);
00304 rtn.abs();
00305 }
00306 return rtn;
00307 }
00308
00309
00310
00311
00312
00313
00314 template <class Scalar> inline
00315 Vector<Scalar> Vector<Scalar>::reciprocal() const
00316 {
00317 Vector<Scalar> rtn = space().createMember();
00318 {
00319 TimeMonitor t(*opTimer());
00320 rtn.acceptCopyOf(*this);
00321 rtn.reciprocal();
00322 }
00323 return rtn;
00324 }
00325
00326
00327
00328
00329 template <class Scalar> inline
00330 Vector<Scalar>& Vector<Scalar>::abs()
00331 {
00332 {
00333 TimeMonitor t(*opTimer());
00334 Thyra::abs(this->ptr().ptr(), *(this->ptr()));
00335 }
00336 return *this;
00337 }
00338
00339
00340
00341
00342
00343
00344 template <class Scalar> inline
00345 Vector<Scalar>& Vector<Scalar>::reciprocal()
00346 {
00347 TimeMonitor t(*opTimer());
00348 Thyra::reciprocal(this->ptr().ptr(), *(this->ptr()));
00349
00350 return *this;
00351 }
00352
00353
00354
00355
00356 template <class Scalar> inline
00357 Vector<Scalar>& Vector<Scalar>::update(const Scalar& alpha,
00358 const Vector<Scalar>& x,
00359 const Scalar& gamma)
00360 {
00361 TimeMonitor t(*opTimer());
00362
00363 ArrayView<const Scalar> a(&alpha, 1);
00364 Teuchos::Ptr<const Thyra::VectorBase<Scalar> > px = x.ptr().ptr();
00365 ArrayView<const Teuchos::Ptr<const Thyra::VectorBase<Scalar> > > apx(&px,1);
00366
00367 Thyra::linear_combination(a, apx, gamma, this->ptr().ptr());
00368
00369 return *this;
00370 }
00371
00372
00373
00374
00375
00376 template <class Scalar> inline
00377 Vector<Scalar>& Vector<Scalar>::update(const Scalar& alpha,
00378 const Vector<Scalar>& x,
00379 const Scalar& beta,
00380 const Vector<Scalar>& y,
00381 const Scalar& gamma)
00382 {
00383 TimeMonitor t(*opTimer());
00384
00385 Scalar a[2];
00386 a[0] = alpha;
00387 a[1] = beta;
00388 ArrayView<const Scalar> av(a,2);
00389
00390 Ptr<const Thyra::VectorBase<Scalar> > vecs[2];
00391 vecs[0] = x.ptr().ptr().getConst();
00392 vecs[1] = y.ptr().ptr().getConst();
00393 ArrayView<const Ptr<const Thyra::VectorBase<Scalar> > > vv(vecs,2);
00394
00395 Thyra::linear_combination(av, vv, gamma, this->ptr().ptr());
00396
00397 return *this;
00398 }
00399
00400
00401
00402
00403
00404 template <class Scalar> inline
00405 Scalar Vector<Scalar>::dot(const Vector<Scalar>& other) const
00406 {
00407 TimeMonitor t(*opTimer());
00408
00409 return Thyra::dot(*(this->ptr)(), *(other.ptr()));
00410 }
00411
00412
00413
00414 template <class Scalar> inline
00415 Scalar Vector<Scalar>::operator*(const Vector<Scalar>& other) const
00416 {
00417 return dot(other);
00418 }
00419
00420
00421
00422
00423
00424 template <class Scalar> inline
00425 Scalar Vector<Scalar>::norm1() const
00426 {
00427 TimeMonitor t(*opTimer());
00428
00429 return Thyra::norm_1(*(this->ptr)());
00430 }
00431
00432
00433
00434
00435
00436 template <class Scalar> inline
00437 Scalar Vector<Scalar>::norm2() const
00438 {
00439 TimeMonitor t(*opTimer());
00440 return Thyra::norm_2(*(this->ptr)());
00441 }
00442
00443
00444
00445
00446
00447 template <class Scalar> inline
00448 Scalar Vector<Scalar>::norm2(const Vector<Scalar>& weights) const
00449 {
00450 TimeMonitor t(*opTimer());
00451
00452 return Thyra::norm_2(*(weights.ptr()), *(this->ptr)());
00453 }
00454
00455
00456
00457
00458
00459
00460 template <class Scalar> inline
00461 Scalar Vector<Scalar>::normInf() const
00462 {
00463 TimeMonitor t(*opTimer());
00464
00465 return Thyra::norm_inf(*(this->ptr)());
00466 }
00467
00468
00469
00470
00471
00472 template <class Scalar> inline
00473 void Vector<Scalar>::zero()
00474 {
00475 TimeMonitor t(*opTimer());
00476
00477 Thyra::assign(this->ptr().ptr(), 0.0);
00478 }
00479
00480
00481
00482
00483
00484 template <class Scalar> inline
00485 void Vector<Scalar>::setToConstant(const Scalar& alpha)
00486 {
00487 TimeMonitor t(*opTimer());
00488
00489 Thyra::assign(this->ptr().ptr(), alpha);
00490 }
00491
00492
00493
00494
00495
00496
00497 template <class Scalar> inline
00498 Scalar Vector<Scalar>::max()const
00499 {
00500 TimeMonitor t(*opTimer());
00501 return Thyra::max(*(this->ptr)());
00502 }
00503
00504
00505
00506 template <class Scalar> inline
00507 Scalar Vector<Scalar>::max(int& index)const
00508 {
00509 TimeMonitor t(*opTimer());
00510 Scalar maxEl;
00511 Ptr<Scalar> maxElP(&maxEl);
00512 OrdType loc_index = -1;
00513 Ptr<OrdType> pind(&loc_index);
00514 Thyra::max(*(this->ptr()), maxElP, pind);
00515 index = loc_index;
00516 return maxEl;
00517 }
00518
00519
00520
00521 template <class Scalar> inline
00522 Scalar Vector<Scalar>::max(const Scalar& bound, int& index)const
00523 {
00524 TimeMonitor t(*opTimer());
00525 Scalar maxEl;
00526 Ptr<Scalar> maxElP(&maxEl);
00527 OrdType loc_index = -1;
00528 Ptr<OrdType> pind(&loc_index);
00529 Thyra::maxLessThanBound(*(this->ptr()), bound, maxElP, pind);
00530 index = loc_index;
00531 return maxEl;
00532
00533 }
00534
00535
00536
00537 template <class Scalar> inline
00538 Scalar Vector<Scalar>::min()const
00539 {
00540 TimeMonitor t(*opTimer());
00541 return Thyra::min(*(this->ptr()));
00542 }
00543
00544
00545
00546 template <class Scalar> inline
00547 Scalar Vector<Scalar>::min(int& index)const
00548 {
00549 TimeMonitor t(*opTimer());
00550 Scalar minEl;
00551 Ptr<Scalar> minElP(&minEl);
00552 OrdType loc_index = -1;
00553 Ptr<OrdType> pind(&loc_index);
00554 Thyra::min(*(this->ptr()), minElP, pind);
00555 index = loc_index;
00556 return minEl;
00557 }
00558
00559
00560
00561 template <class Scalar> inline
00562 Scalar Vector<Scalar>::min(const Scalar& bound, int& index)const
00563 {
00564 TimeMonitor t(*opTimer());
00565 Scalar minEl;
00566 Ptr<Scalar> minElP(&minEl);
00567 OrdType loc_index = -1;
00568 Ptr<OrdType> pind(&loc_index);
00569
00570 Thyra::minGreaterThanBound(*(this->ptr)(), bound, minElP, pind);
00571 index = loc_index;
00572 return minEl;
00573 }
00574
00575
00576
00577
00578
00579 template <class Scalar> inline
00580 Scalar Vector<Scalar>::getElement(OrdType globalIndex) const
00581 {
00582 Thyra::ProductVectorBase<Scalar>* p
00583 = dynamic_cast<Thyra::ProductVectorBase<Scalar>*>(&*this->ptr());
00584
00585 if (p)
00586 {
00587 const Thyra::ProductVectorSpaceBase<Scalar>* pvs
00588 = dynamic_cast<const Thyra::ProductVectorSpaceBase<Scalar>*>(space().ptr().get());
00589 TEST_FOR_EXCEPT(pvs == 0);
00590
00591 const Thyra::DefaultProductVectorSpace<Scalar>* dpvs
00592 = dynamic_cast<const Thyra::DefaultProductVectorSpace<Scalar>*>(pvs);
00593 if (dpvs)
00594 {
00595 int blockIndex=-1;
00596 OrdType globOffsetInBlock=-1;
00597 dpvs->getVecSpcPoss(globalIndex, &blockIndex, &globOffsetInBlock);
00598 RCP<Thyra::VectorBase<Scalar> > vec_i
00599 = p->getNonconstVectorBlock(blockIndex);
00600 Vector<Scalar> vv(vec_i);
00601 return vv.getElement(globOffsetInBlock);
00602 }
00603 else
00604 {
00605 int k = 0;
00606 for (int i = 0; i < pvs->numBlocks(); i++)
00607 {
00608 RCP<Thyra::VectorBase<Scalar> > vec_i
00609 = p->getNonconstVectorBlock(i);
00610 int len = vec_i->space()->dim();
00611 if (globalIndex < k + len )
00612 {
00613 Vector<Scalar> vv(vec_i);
00614 int globalIndexWithinBlock = globalIndex - k;
00615 return vv.getElement(globalIndexWithinBlock);
00616 break;
00617 }
00618 k += len;
00619 }
00620 }
00621 TEST_FOR_EXCEPT(true);
00622 return 0.0;
00623 }
00624 else
00625 {
00626 Thyra::DefaultSpmdVector<Scalar>* dsv
00627 = dynamic_cast<Thyra::DefaultSpmdVector<Scalar>*>(this->ptr().get());
00628
00629 if (dsv)
00630 {
00631 OrdType stride = dsv->getStride();
00632 OrdType low = dsv->spmdSpace()->localOffset();
00633 OrdType subdim = dsv->spmdSpace()->localSubDim();
00634 TEST_FOR_EXCEPTION( globalIndex < low || globalIndex >= low+subdim,
00635 std::runtime_error,
00636 "Bounds violation: " << globalIndex << "is out of range [low"
00637 << ", " << low+subdim << "]");
00638 return dsv->getPtr()[stride*(globalIndex - low)];
00639 }
00640 else
00641 {
00642 return castToAccessible()->getElement(globalIndex);
00643 }
00644 }
00645 }
00646
00647
00648 template <class Scalar> inline
00649 void Vector<Scalar>::setElement(OrdType globalIndex, const Scalar& value)
00650 {
00651 Thyra::ProductVectorBase<Scalar>* p
00652 = dynamic_cast<Thyra::ProductVectorBase<Scalar>*>(&*this->ptr());
00653
00654 if (p)
00655 {
00656 const Thyra::ProductVectorSpaceBase<Scalar>* pvs
00657 = dynamic_cast<const Thyra::ProductVectorSpaceBase<Scalar>*>(space().ptr().get());
00658 TEST_FOR_EXCEPT(pvs == 0);
00659
00660 const Thyra::DefaultProductVectorSpace<Scalar>* dpvs
00661 = dynamic_cast<const Thyra::DefaultProductVectorSpace<Scalar>*>(pvs);
00662 if (dpvs)
00663 {
00664 int blockIndex=-1;
00665 OrdType globOffsetInBlock=-1;
00666 dpvs->getVecSpcPoss(globalIndex, &blockIndex, &globOffsetInBlock);
00667 RCP<Thyra::VectorBase<Scalar> > vec_i
00668 = p->getNonconstVectorBlock(blockIndex);
00669 Vector<Scalar> vv(vec_i);
00670 vv.setElement(globOffsetInBlock, value);
00671 }
00672 else
00673 {
00674 int k = 0;
00675 for (int i = 0; i < pvs->numBlocks(); i++)
00676 {
00677 RCP<Thyra::VectorBase<Scalar> > vec_i
00678 = p->getNonconstVectorBlock(i);
00679 int len = vec_i->space()->dim();
00680 if (globalIndex < k + len )
00681 {
00682 Vector<Scalar> vv(vec_i);
00683 int globalIndexWithinBlock = globalIndex - k;
00684 vv.setElement(globalIndexWithinBlock, value);
00685 break;
00686 }
00687 k += len;
00688 }
00689 }
00690
00691 }
00692 else
00693 {
00694 Thyra::DefaultSpmdVector<Scalar>* dsv
00695 = dynamic_cast<Thyra::DefaultSpmdVector<Scalar>*>(this->ptr().get());
00696
00697 if (dsv)
00698 {
00699 OrdType stride = dsv->getStride();
00700 OrdType low = dsv->spmdSpace()->localOffset();
00701 OrdType subdim = dsv->spmdSpace()->localSubDim();
00702 TEST_FOR_EXCEPTION( globalIndex < low || globalIndex >= low+subdim,
00703 std::runtime_error,
00704 "Bounds violation: " << globalIndex << "is out of range [low"
00705 << ", " << low+subdim << "]");
00706 dsv->getPtr()[stride*(globalIndex - low)] = value;
00707 }
00708 else
00709 {
00710 castToLoadable()->setElement(globalIndex, value);
00711 }
00712 }
00713 }
00714
00715
00716
00717 template <class Scalar> inline
00718 Scalar& Vector<Scalar>::operator[](const SequentialIterator<Scalar>& iter)
00719 {
00720 const OrdType& blockIndex = iter.blockIndex();
00721 const OrdType& indexInBlock = iter.indexInBlock();
00722
00723 return localElement(blockIndex, indexInBlock);
00724 }
00725
00726
00727
00728
00729 template <class Scalar> inline
00730 const Scalar& Vector<Scalar>::operator[](const SequentialIterator<Scalar>& iter) const
00731 {
00732 const OrdType& blockIndex = iter.blockIndex();
00733 const OrdType& indexInBlock = iter.indexInBlock();
00734
00735 return localElement(blockIndex, indexInBlock);
00736 }
00737
00738
00739
00740 template <class Scalar> inline
00741 const Scalar& Vector<Scalar>::localElement(const OrdType& blockIndex, const OrdType& indexInBlock) const
00742 {
00743 const Thyra::ProductVectorBase<Scalar>* p
00744 = dynamic_cast<const Thyra::ProductVectorBase<Scalar>*>(&*this->ptr());
00745 RCP<const Thyra::VectorBase<Scalar> > vec;
00746 if (p)
00747 {
00748 vec = p->getVectorBlock(blockIndex);
00749 }
00750 else
00751 {
00752 vec = this->ptr();
00753 }
00754 const Thyra::DefaultSpmdVector<Scalar>* dsv
00755 = dynamic_cast<const Thyra::DefaultSpmdVector<Scalar>*>(this->ptr().get());
00756 if (dsv != 0)
00757 {
00758 return dsv->getPtr()[indexInBlock*dsv->getStride()];
00759 }
00760 return castToRawDataAccessible()->dataPtr()[indexInBlock];
00761 }
00762
00763
00764
00765
00766 template <class Scalar> inline
00767 Scalar& Vector<Scalar>::localElement(const OrdType& blockIndex, const OrdType& indexInBlock)
00768 {
00769 Thyra::ProductVectorBase<Scalar>* p
00770 = dynamic_cast<Thyra::ProductVectorBase<Scalar>*>(&*this->ptr());
00771 RCP<Thyra::VectorBase<Scalar> > vec;
00772 if (p)
00773 {
00774 vec = p->getNonconstVectorBlock(blockIndex);
00775 }
00776 else
00777 {
00778 vec = this->ptr();
00779 }
00780
00781 Thyra::DefaultSpmdVector<Scalar>* dsv
00782 = dynamic_cast<Thyra::DefaultSpmdVector<Scalar>*>(this->ptr().get());
00783 if (dsv != 0)
00784 {
00785 return dsv->getPtr()[indexInBlock*dsv->getStride()];
00786 }
00787 TSFExtended::RawDataAccessibleVector<Scalar>* d
00788 = this->castToRawDataAccessible();
00789 return d->dataPtr()[indexInBlock];
00790 }
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800 template <class Scalar> inline
00801 void Vector<Scalar>::addToElement(OrdType globalIndex, const Scalar& value)
00802 {
00803 Thyra::ProductVectorBase<Scalar>* p
00804 = dynamic_cast<Thyra::ProductVectorBase<Scalar>*>(&*this->ptr());
00805
00806 if (p)
00807 {
00808 const Thyra::ProductVectorSpaceBase<Scalar>* pvs
00809 = dynamic_cast<const Thyra::ProductVectorSpaceBase<Scalar>*>(space().ptr().get());
00810 TEST_FOR_EXCEPT(pvs == 0);
00811
00812 const Thyra::DefaultProductVectorSpace<Scalar>* dpvs
00813 = dynamic_cast<const Thyra::DefaultProductVectorSpace<Scalar>*>(pvs);
00814 if (dpvs)
00815 {
00816 int blockIndex=-1;
00817 OrdType globOffsetInBlock=-1;
00818 dpvs->getVecSpcPoss(globalIndex, &blockIndex, &globOffsetInBlock);
00819 RCP<Thyra::VectorBase<Scalar> > vec_i
00820 = p->getNonconstVectorBlock(blockIndex);
00821 Vector<Scalar> vv(vec_i);
00822 vv.addToElement(globOffsetInBlock, value);
00823 }
00824 else
00825 {
00826 int k = 0;
00827 for (int i = 0; i < pvs->numBlocks(); i++)
00828 {
00829 RCP<Thyra::VectorBase<Scalar> > vec_i
00830 = p->getNonconstVectorBlock(i);
00831 int len = vec_i->space()->dim();
00832 if (globalIndex < k + len )
00833 {
00834 Vector<Scalar> vv(vec_i);
00835 int globalIndexWithinBlock = globalIndex - k;
00836 vv.addToElement(globalIndexWithinBlock, value);
00837 break;
00838 }
00839 k += len;
00840 }
00841 }
00842
00843 }
00844 else
00845 {
00846 Thyra::DefaultSpmdVector<Scalar>* dsv
00847 = dynamic_cast<Thyra::DefaultSpmdVector<Scalar>*>(this->ptr().get());
00848
00849 LoadableVector<Scalar>* loadable = dynamic_cast<LoadableVector<Scalar>*>(this->ptr().get());
00850 if (loadable)
00851 {
00852 loadable->addToElement(globalIndex, value);
00853 }
00854 else if (dsv)
00855 {
00856 OrdType stride = dsv->getStride();
00857 OrdType low = dsv->spmdSpace()->localOffset();
00858 OrdType subdim = dsv->spmdSpace()->localSubDim();
00859 TEST_FOR_EXCEPTION( globalIndex < low || globalIndex >= low+subdim,
00860 std::runtime_error,
00861 "Bounds violation: " << globalIndex << "is out of range [low"
00862 << ", " << low+subdim << "]");
00863 dsv->getPtr()[stride*(globalIndex - low)] += value;
00864 }
00865 else
00866 {
00867 TEST_FOR_EXCEPT(true);
00868 }
00869 }
00870 }
00871
00872
00873
00874 template <class Scalar> inline
00875 void Vector<Scalar>::boundscheck(OrdType i, int dim) const
00876 {
00877 TEST_FOR_EXCEPTION( i < 0 || i >= dim, std::runtime_error,
00878 "Bounds violation: " << i << "is out of range [0"
00879 << ", " << dim << "]");
00880 }
00881
00882
00883 }
00884
00885 #endif
