00001
00002
00003
00004
00005 #ifndef PLAYA_SINGLECHUNKVECTOR_HPP
00006 #define PLAYA_SINGLECHUNKVECTOR_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaVectorBaseDecl.hpp"
00010
00011 extern "C"
00012 {
00013 void daxpy_(int*, double*, double*, int*, double*, int*);
00014 }
00015
00016
00017 namespace Playa
00018 {
00019
00020
00021
00022
00023
00024
00025 template <class Scalar>
00026 class SingleChunkVector : public VectorBase<Scalar>
00027 {
00028 public:
00029
00030 SingleChunkVector() : rewound_(true) {}
00031
00032 virtual ~SingleChunkVector() {}
00033
00034
00035 virtual ConstDataChunk<Scalar> nextConstChunk() const
00036 {rewound_ = false; return ConstDataChunk<Scalar>(chunkSize(), dataPtr());}
00037
00038
00039 virtual NonConstDataChunk<Scalar> nextChunk()
00040 {rewound_ = false; return NonConstDataChunk<Scalar>(chunkSize(), dataPtr());}
00041
00042
00043 virtual bool hasMoreChunks() const
00044 {return rewound_;}
00045
00046
00047 virtual void rewind() const
00048 {rewound_=true;}
00049
00050
00051
00052
00053 virtual const double& operator[](int localIndex) const
00054 {return dataPtr()[localIndex];}
00055
00056
00057 virtual double& operator[](int localIndex)
00058 {return dataPtr()[localIndex];}
00059
00060
00061
00062 virtual int chunkSize() const = 0 ;
00063
00064 virtual const Scalar* dataPtr() const = 0 ;
00065
00066 virtual Scalar* dataPtr() = 0 ;
00067
00068 virtual void update(const Scalar& alpha, const VectorBase<Scalar>* other,
00069 const Scalar& gamma)
00070 {
00071 const SingleChunkVector<Scalar>* sco
00072 = dynamic_cast<const SingleChunkVector<Scalar>* >(other);
00073 TEUCHOS_TEST_FOR_EXCEPT(sco==0);
00074
00075 Scalar* const myVals = this->dataPtr();
00076 const Scalar* const yourVals = sco->dataPtr();
00077 int n = chunkSize();
00078 int stride = 1;
00079 if (gamma==1.0)
00080 {
00081 daxpy_(&n, (double*) &alpha, (double*) yourVals, &stride,
00082 myVals, &stride);
00083 }
00084 else
00085 {
00086 for (int i=0; i<n; i++)
00087 {
00088 myVals[i] = gamma*myVals[i] + alpha*yourVals[i];
00089 }
00090 }
00091 }
00092
00093
00094 virtual void update(
00095 const Scalar& alpha, const VectorBase<Scalar>* x,
00096 const Scalar& beta, const VectorBase<Scalar>* y,
00097 const Scalar& gamma)
00098 {
00099 const SingleChunkVector<Scalar>* scx
00100 = dynamic_cast<const SingleChunkVector<Scalar>* >(x);
00101 TEUCHOS_TEST_FOR_EXCEPT(scx==0);
00102 const SingleChunkVector<Scalar>* scy
00103 = dynamic_cast<const SingleChunkVector<Scalar>* >(y);
00104 TEUCHOS_TEST_FOR_EXCEPT(scy==0);
00105
00106 Scalar* const myVals = this->dataPtr();
00107 const Scalar* const xVals = scx->dataPtr();
00108 const Scalar* const yVals = scy->dataPtr();
00109 int n = chunkSize();
00110 if (gamma==1.0)
00111 {
00112 for (int i=0; i<n; i++)
00113 {
00114 myVals[i] += alpha*xVals[i] + beta*yVals[i];
00115 }
00116 }
00117 else
00118 {
00119 for (int i=0; i<n; i++)
00120 {
00121 myVals[i] = gamma*myVals[i] + alpha*xVals[i] + beta*yVals[i];
00122 }
00123 }
00124 }
00125
00126
00127 virtual void update(
00128 const Scalar& alpha, const VectorBase<Scalar>* x,
00129 const Scalar& beta, const VectorBase<Scalar>* y,
00130 const Scalar& gamma, const VectorBase<Scalar>* z,
00131 const Scalar& delta)
00132 {
00133 const SingleChunkVector<Scalar>* scx
00134 = dynamic_cast<const SingleChunkVector<Scalar>* >(x);
00135 TEUCHOS_TEST_FOR_EXCEPT(scx==0);
00136 const SingleChunkVector<Scalar>* scy
00137 = dynamic_cast<const SingleChunkVector<Scalar>* >(y);
00138 TEUCHOS_TEST_FOR_EXCEPT(scy==0);
00139 const SingleChunkVector<Scalar>* scz
00140 = dynamic_cast<const SingleChunkVector<Scalar>* >(z);
00141 TEUCHOS_TEST_FOR_EXCEPT(scz==0);
00142
00143 Scalar* const myVals = this->dataPtr();
00144 const Scalar* const xVals = scx->dataPtr();
00145 const Scalar* const yVals = scy->dataPtr();
00146 const Scalar* const zVals = scz->dataPtr();
00147
00148 int n = chunkSize();
00149 if (delta==1.0)
00150 {
00151 for (int i=0; i<n; i++)
00152 {
00153 myVals[i] += alpha*xVals[i] + beta*yVals[i] + gamma*zVals[i];
00154 }
00155 }
00156 else
00157 {
00158 for (int i=0; i<n; i++)
00159 {
00160 myVals[i] = delta*myVals[i] + alpha*xVals[i] + beta*yVals[i]
00161 + gamma*zVals[i];
00162 }
00163 }
00164 }
00165
00166
00167
00168
00169 virtual Scalar dot(const VectorBase<Scalar>* other) const
00170 {
00171 const SingleChunkVector<Scalar>* const sco
00172 = dynamic_cast<const SingleChunkVector<Scalar>* >(other);
00173 TEUCHOS_TEST_FOR_EXCEPT(sco==0);
00174
00175 const Scalar* const yourVals = sco->dataPtr();
00176 const Scalar* const myVals = this->dataPtr();
00177
00178 Scalar rtn = 0.0;
00179 int n = this->chunkSize();
00180 for (int i=0; i<n; i++) rtn += myVals[i]*yourVals[i];
00181 return rtn;
00182 }
00183
00184
00185 virtual Scalar norm2() const
00186 {
00187 const Scalar* const myVals = this->dataPtr();
00188
00189 Scalar rtn = 0.0;
00190 int n = this->chunkSize();
00191 for (int i=0; i<n; i++) rtn += myVals[i]*myVals[i];
00192 return ::sqrt(rtn);
00193 }
00194
00195 private:
00196 mutable bool rewound_;
00197 };
00198 }
00199
00200
00201 #endif