00001
00002
00003
00004
00005 #ifndef PLAYA_VECTORSPACEIMPL_HPP
00006 #define PLAYA_VECTORSPACEIMPL_HPP
00007
00008 #include "PlayaVectorSpaceDecl.hpp"
00009 #include "PlayaBlockVectorSpaceDecl.hpp"
00010 #include "PlayaVectorDecl.hpp"
00011 #include "PlayaExceptions.hpp"
00012 #include "Teuchos_Time.hpp"
00013 #include "Teuchos_TimeMonitor.hpp"
00014
00015 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00016 #include "PlayaBlockVectorSpaceImpl.hpp"
00017 #endif
00018
00019 using namespace Teuchos;
00020
00021 namespace Playa
00022 {
00023
00024
00025 static inline Time& createVecTimer()
00026 {
00027 static RCP<Time> rtn
00028 = TimeMonitor::getNewTimer("vector allocation");
00029 return *rtn;
00030 }
00031
00032
00033
00034 template <class Scalar> inline
00035 bool VectorSpace<Scalar>::operator==(const VectorSpace<Scalar>& other) const
00036 {
00037 return isCompatible(other);
00038 }
00039
00040
00041
00042 template <class Scalar> inline
00043 bool VectorSpace<Scalar>::operator!=(const VectorSpace<Scalar>& other) const
00044 {
00045 return !(operator==(other));
00046 }
00047
00048
00049
00050
00051 template <class Scalar> inline
00052 Vector<Scalar> VectorSpace<Scalar>::createMember() const
00053 {
00054 Vector<Scalar> rtn = this->ptr()->createMember(*this);
00055 rtn.setToConstant(0.0);
00056 return rtn;
00057 }
00058
00059
00060
00061
00062 template <class Scalar> inline
00063 int VectorSpace<Scalar>::baseGlobalNaturalIndex() const
00064 {
00065 return this->ptr()->baseGlobalNaturalIndex();
00066 }
00067
00068
00069 template <class Scalar> inline
00070 int VectorSpace<Scalar>::numLocalElements() const
00071 {
00072 return this->ptr()->numLocalElements();
00073 }
00074
00075
00076
00077
00078
00079 template <class Scalar> inline
00080 bool VectorSpace<Scalar>::isCompatible(const VectorSpace<Scalar>& vecSpc) const
00081 {
00082 TEUCHOS_TEST_FOR_EXCEPTION(vecSpc.ptr().get() == 0, std::runtime_error,
00083 "null argument in VectorSpace<Scalar>::isCompatible()");
00084 return this->ptr().get()->isCompatible(vecSpc.ptr().get());
00085 }
00086
00087
00088
00089
00090
00091
00092 template <class Scalar> inline
00093 bool VectorSpace<Scalar>::contains(const Vector<Scalar> &vec) const
00094 {
00095 return (operator==(vec.space()));
00096 }
00097
00098
00099
00100 template <class Scalar>
00101 int VectorSpace<Scalar>::numBlocks() const
00102 {
00103 return this->ptr()->numBlocks();
00104 }
00105
00106
00107
00108 template <class Scalar>
00109 bool VectorSpace<Scalar>::isBlockSpace() const
00110 {
00111 return dynamic_cast<const BlockVectorSpaceBase<Scalar>*>(this->ptr().get())!=0;
00112 }
00113
00114
00115
00116
00117 template <class Scalar> inline
00118 const VectorSpace<Scalar>& VectorSpace<Scalar>::getBlock(const int i) const
00119 {
00120 const BlockVectorSpaceBase<Scalar>* bvs =
00121 dynamic_cast<const BlockVectorSpaceBase<Scalar>* > (this->ptr().get());
00122 TEUCHOS_TEST_FOR_EXCEPTION(bvs == 0 && numBlocks()!=1, std::runtime_error,
00123 "getBlock called for a space that "
00124 "is not a BlockVectorSpace" << std::endl);
00125 if (bvs != 0)
00126 {
00127 return bvs->getBlock(i);
00128 }
00129 return *this;
00130 }
00131
00132
00133 template <class Scalar> inline
00134 const VectorSpace<Scalar>& VectorSpace<Scalar>
00135 ::getBlock(const BlockIterator<Scalar>& b) const
00136 {
00137
00138 TEUCHOS_TEST_FOR_EXCEPTION(b.atEnd(), RuntimeError,
00139 "Attempt to use a block iterator that's run off end");
00140
00141 return this->getBlock(b.blockIndex());
00142 }
00143
00144
00145
00146 template <class Scalar> inline
00147 const VectorSpace<Scalar>& VectorSpace<Scalar>
00148 ::getBlock(const std::deque<int>& b) const
00149 {
00150
00151 if (b.size()==0) return *this;
00152
00153 if (b.size()==1)
00154 {
00155 return this->getBlock(b.front());
00156 }
00157
00158 int b0 = b.front();
00159 std::deque<int> tmp = b;
00160 tmp.pop_front();
00161 return this->getBlock(b0).getBlock(tmp);
00162 }
00163
00164
00165
00166
00167
00168 template <class Scalar> inline
00169 BlockIterator<Scalar> VectorSpace<Scalar>::beginBlock() const
00170 {
00171 return BlockIterator<Scalar>(*this, false);
00172 }
00173
00174 template <class Scalar> inline
00175 BlockIterator<Scalar> VectorSpace<Scalar>::endBlock() const
00176 {
00177 return BlockIterator<Scalar>(*this, true);
00178 }
00179
00180
00181 template <class Scalar> inline
00182 int VectorSpace<Scalar>::mapToGNI(
00183 const BlockIterator<Scalar>& b,
00184 int indexWithinBlock) const
00185 {
00186 int rtn = this->baseGlobalNaturalIndex();
00187 if (this->numBlocks() > 1)
00188 {
00189 for (BlockIterator<Scalar> a=this->beginBlock(); a < b; a++)
00190 {
00191 rtn += getBlock(a).numLocalElements();
00192 }
00193 }
00194 return rtn + indexWithinBlock;
00195 }
00196
00197 template <class Scalar> inline
00198 bool VectorSpace<Scalar>::containsGNI(int gni) const
00199 {
00200 return (gni >= this->baseGlobalNaturalIndex()
00201 && gni < this->baseGlobalNaturalIndex()+this->numLocalElements());
00202 }
00203
00204 template <class Scalar> inline
00205 void VectorSpace<Scalar>::getBlockAndOffsetFromGNI(int gni,
00206 BlockIterator<Scalar>& block, int& indexWithinBlock) const
00207 {
00208 int low = this->baseGlobalNaturalIndex();
00209
00210 if (this->numBlocks() == 1)
00211 {
00212 indexWithinBlock = gni - low;
00213 }
00214 else
00215 {
00216 int sizeSum = low;
00217 for (BlockIterator<Scalar> b=this->beginBlock();
00218 b != this->endBlock(); b++)
00219 {
00220 if (this->getBlock(b).containsGNI(gni))
00221 {
00222 block = b;
00223 indexWithinBlock = gni - sizeSum;
00224 }
00225 sizeSum += this->getBlock(b).numLocalElements();
00226 }
00227 }
00228 }
00229
00230 }
00231
00232
00233
00234 #endif