Sundance 2.1

00001 #ifndef TSFDENSESERIALVECTOR_H
00002 #define TSFDENSESERIALVECTOR_H
00003 
00004 #include "TSFConfigDefs.hpp"
00005 #include "Teuchos_TestForException.hpp"
00006 #include "Teuchos_Array.hpp"
00007 #include "Teuchos_BLAS.hpp"
00008 
00009 namespace TSFExtended
00010 {
00011   using std::string;
00012 
00013   /**\ingroup DenseSerial
00014    * Serial vector with math operations
00015    */
00016 
00017   class DenseSerialVector : public Teuchos::Array<double>
00018   {
00019   public:
00020     /** {\bf Constructors, Destructors, and Assignment Operator} */
00021     //@{
00022     /** Empty ctor */
00023     DenseSerialVector() : Teuchos::Array<double>() {;}
00024     /** Create a vector of length n */
00025     inline DenseSerialVector(int n) : Teuchos::Array<double>(n) {;}
00026     /** Create a vector of length n, and fill it with value */
00027     inline DenseSerialVector(int n, const double& value);
00028     /** Create a vector of length n, assuming responsibility for a C array */
00029     inline DenseSerialVector(int n, const double* cArray);
00030     //@}
00031 
00032 
00033     /** {\bf some reasonably efficient math operations} */
00034     //@{
00035     /** change sign */
00036     void negate();
00037 
00038     /** vector addition with result returned through a reference argument */
00039     inline void add(const DenseSerialVector& other,
00040                     DenseSerialVector& result) const ;
00041     /** self-modifying vector addition */
00042     void add(const DenseSerialVector& other) ;
00043 
00044     /** vector subtraction with result returned through a reference argument */
00045     inline void subtract(const DenseSerialVector& other,
00046                          DenseSerialVector& result) const ;
00047     /** self-modifying vector subtraction */
00048     void subtract(const DenseSerialVector& other) ;
00049 
00050     /** daxpy (z = a*x + y) with result returned through a reference argument */
00051     inline void daxpy(const DenseSerialVector& other, const double& a,
00052                       DenseSerialVector& result) const ;
00053     /** self-modifying daxpy */
00054     void daxpy(const DenseSerialVector& other, const double& a) ;
00055 
00056     /** element-by-element multiplication
00057      * with result returned through a reference argument */
00058     inline void eMult(const DenseSerialVector& other,
00059                       DenseSerialVector& result) const ;
00060 
00061     /** self-modifying element-by-element multiplication */
00062     void eMult(const DenseSerialVector& other) ;
00063 
00064 
00065     /** absolute value of each element */
00066     void abs()  ;
00067 
00068     /** return the value of the max element  */
00069     double max() const ;
00070 
00071     /** return the value of the min element  */
00072     double min() const ;
00073 
00074     /** compute the matlab style ".*" operation, i.e.,
00075      *      this[i] = y[i] * z[i]  */
00076     void dotStar(const DenseSerialVector& y,
00077                  const DenseSerialVector& z);
00078 
00079 
00080     /** compute the matlab style "./" operation, i.e.,
00081      *      this[i] = y[i] / z[i]  */
00082     void dotSlash(const DenseSerialVector& y,
00083                   const DenseSerialVector& z);
00084 
00085 
00086     /** multiplication by a scalar
00087      * with result returned through a reference argument */
00088     inline void scalarMult(const double& scalar,
00089                            DenseSerialVector& result) const ;
00090     /** self-modifying multiplication by a scalar */
00091     void scalarMult(const double& scalar);
00092 
00093     /** exponentiation by a scalar
00094      * with result returned through a reference argument */
00095     inline void scalarPow(const double& scalar,
00096                           DenseSerialVector& result ) const ;
00097     /** Self-modifying exponentiation by a scalar */
00098     void scalarPow(const double& scalar);
00099 
00100     /** dot product */
00101     double dot(const DenseSerialVector& other) const ;
00102 
00103     /** dot product with self */
00104     double norm2Squared() const ;
00105 
00106     /** 2-norm */
00107     double norm2() const {return sqrt(norm2Squared());}
00108 
00109     /** sum elements */
00110     inline double sumElements() const ;
00111 
00112     /** return maximum element value */
00113     double maxNorm() const ;
00114 
00115     /** set all elements to zero */
00116     void zero() {setScalar(0.0);}
00117 
00118     /** set all elements to the given value */
00119     void setScalar(const double& a);
00120     //@}
00121 
00122     /** {\bf overloaded math operators; for performance reasons,
00123      * avoid these in performance-critical
00124      * code} */
00125     //@{
00126     /** unary minus */
00127     inline DenseSerialVector operator-() const ;
00128     /** reflexive addition */
00129     inline DenseSerialVector& operator+=(const DenseSerialVector& other);
00130     /** reflexive subtraction */
00131     inline DenseSerialVector& operator-=(const DenseSerialVector& other);
00132     /** reflexive scalar mult */
00133     inline DenseSerialVector& operator*=(const double& scalar);
00134     /** reflexive scalar division */
00135     inline DenseSerialVector& operator/=(const double& scalar);
00136 
00137     /** addition */
00138     inline DenseSerialVector operator+(const DenseSerialVector& other) const ;
00139     /** subtraction */
00140     inline DenseSerialVector operator-(const DenseSerialVector& other) const ;
00141     /** dot product */
00142     inline double operator*(const DenseSerialVector& other) const ;
00143     /** scalar mult */
00144     inline DenseSerialVector operator*(const double& scalar) const ;
00145     /** scalar division */
00146     inline DenseSerialVector operator/(const double& scalar) const ;
00147     //@}
00148 
00149     /** write a brief description to string */
00150     string summary() const ;
00151 
00152     /** a BLAS object */
00153     inline static const Teuchos::BLAS<int, double>& blasObject()
00154     {static Teuchos::BLAS<int, double> rtn; return rtn;}
00155 
00156   private:
00157     double* x() {return &(operator[](0));}
00158     const double* x() const {return &(operator[](0));}
00159     void checkLength(const DenseSerialVector& other, 
00160                      const string& funcName) const ;
00161 
00162   };
00163 
00164 
00165   inline DenseSerialVector::DenseSerialVector(int n, const double& value)
00166     : Teuchos::Array<double>(n, value)
00167   {;}
00168 
00169   inline DenseSerialVector::DenseSerialVector(int n, const double* cArray)
00170     :  Teuchos::Array<double>(n)
00171   {
00172     memcpy(x(), (const void*) cArray, (size_t) (n*sizeof(double)));
00173   }
00174 
00175   inline void DenseSerialVector::checkLength(const DenseSerialVector& other, 
00176                                              const string& funcName) const
00177   {
00178     TEST_FOR_EXCEPTION(size() != other.size(), runtime_error,
00179                        "mismatch between operands " 
00180                        << summary() << " and " << other.summary()
00181                        << " in method " << funcName);
00182   }
00183 
00184   inline void DenseSerialVector::add(const DenseSerialVector& other,
00185                                      DenseSerialVector& result) const
00186   {
00187     result = *this;
00188     result.add(other);
00189   }
00190 
00191   inline void DenseSerialVector::subtract(const DenseSerialVector& other,
00192                                           DenseSerialVector& result) const
00193   {
00194     result = *this;
00195     result.subtract(other);
00196   }
00197 
00198   inline void DenseSerialVector::daxpy(const DenseSerialVector& other,
00199                                        const double& a,
00200                                        DenseSerialVector& result) const
00201   {
00202     result = *this;
00203     result.daxpy(other, a);
00204   }
00205 
00206   inline void DenseSerialVector::eMult(const DenseSerialVector& other,
00207                                        DenseSerialVector& result) const
00208   {
00209     result = *this;
00210     result.eMult(other);
00211   }
00212 
00213   inline void DenseSerialVector::scalarMult(const double& a,
00214                                             DenseSerialVector& result) const
00215   {
00216     result = *this;
00217     result.scalarMult(a);
00218   }
00219 
00220 
00221   inline DenseSerialVector DenseSerialVector::operator-() const
00222   {
00223     DenseSerialVector rtn = *this;
00224     rtn.negate();
00225     return rtn;
00226   }
00227 
00228   inline DenseSerialVector& DenseSerialVector::operator+=(const DenseSerialVector& other)
00229   {
00230     add(other);
00231     return *this;
00232   }
00233 
00234   inline DenseSerialVector DenseSerialVector::operator+(const DenseSerialVector& other) const
00235   {
00236     DenseSerialVector rtn = *this;
00237     rtn.add(other);
00238     return rtn;
00239   }
00240 
00241   inline DenseSerialVector& DenseSerialVector::operator-=(const DenseSerialVector& other)
00242   {
00243     subtract(other);
00244     return *this;
00245   }
00246 
00247   inline DenseSerialVector DenseSerialVector::operator-(const DenseSerialVector& other) const
00248   {
00249     DenseSerialVector rtn = *this;
00250     rtn.subtract(other);
00251     return rtn;
00252   }
00253 
00254   inline DenseSerialVector& DenseSerialVector::operator*=(const double& scalar)
00255   {
00256     scalarMult(scalar);
00257     return *this;
00258   }
00259 
00260   inline DenseSerialVector& DenseSerialVector::operator/=(const double& scalar)
00261   {
00262     TEST_FOR_EXCEPTION(scalar==0, std::runtime_error, 
00263                        "DenseSerialVector::operator/= divide by zero");
00264     scalarMult(1.0/scalar);
00265     return *this;
00266   }
00267 
00268   inline DenseSerialVector DenseSerialVector::operator*(const double& scalar) const
00269   {
00270     DenseSerialVector rtn = *this;
00271     rtn.scalarMult(scalar);
00272     return rtn;
00273   }
00274 
00275   inline double DenseSerialVector::operator*(const DenseSerialVector& other) const
00276   {
00277     return dot(other);
00278   }
00279 
00280   inline DenseSerialVector DenseSerialVector::operator/(const double& scalar) const
00281   {
00282     DenseSerialVector rtn = *this;
00283     TEST_FOR_EXCEPTION(scalar==0, std::runtime_error, 
00284                        "DenseSerialVector::operator/ divide by zero");
00285     rtn.scalarMult(1.0/scalar);
00286     return rtn;
00287   }
00288 
00289   inline DenseSerialVector operator*(const double& scalar,
00290                                      const DenseSerialVector& v)
00291   {
00292     return v*scalar;
00293   }
00294 
00295 
00296 }
00297 
00298 
00299 
00300 
00301 #endif
TSFDenseSerialVector.hpp
