From 248560811becefa05c34b685c7cefbf2f7a81070 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 23 Feb 1999 13:54:57 +0000 Subject: [PATCH] Add templated data vectors. They will replace our old dVector eventually. git-svn-id: https://svn.dealii.org/trunk@878 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/vector.h | 494 ++++++++++++++++++ deal.II/lac/include/lac/vector.templates.h | 556 +++++++++++++++++++++ deal.II/lac/source/vector.cc | 7 + 3 files changed, 1057 insertions(+) create mode 100644 deal.II/lac/include/lac/vector.h create mode 100644 deal.II/lac/include/lac/vector.templates.h create mode 100644 deal.II/lac/source/vector.cc diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h new file mode 100644 index 0000000000..7ae125e0f5 --- /dev/null +++ b/deal.II/lac/include/lac/vector.h @@ -0,0 +1,494 @@ +/*---------------------------- vector.h ---------------------------*/ +/* $Id$ */ +#ifndef __vector_H +#define __vector_H +/*---------------------------- vector.h ---------------------------*/ + +// This file was once part of the DEAL Library +// DEAL is Copyright(1995) by +// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier +// Revised and extended by Wolfgang Bangerth + +#include +#include + + + +/** + * Vector of data. + * Memory for Components is supplied explicitly

+ * ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! )

+ * - all necessary methods for Vectors are supplied

+ * - operators available are `=` , `*` and `( )`

+ * CONVENTIONS for used `equations` :

+ * - THIS vector is always named `U`

+ * - vectors are always uppercase , scalars are lowercase + * + * @author Roland Becker, Guido Kanschat, Franz-Theo Suttmeier, revised and extended by Wolfgang Bangerth, documented by Klaus Mampel and Wolfgang Bangerth + */ +template +class Vector { + friend class dFMatrix; + + protected: + + /** + * Dimension. Actual number of components + * contained in the vector. + * Get this number by calling #size()#. + */ + unsigned int dim; + + /** + * Amount of memory actually reserved for + * this vector. This number may be greater + * than #dim# if a #reinit# was called with + * less memory requirements than the vector + * needed last time. At present #reinit# + * does not free memory when the number of + * needed elements is reduced. + */ + unsigned int maxdim; + + /** + * Pointer to the array of components. + */ + Number *val; + + public: + + /** + * Declare iterator types just like those + * for the C++ standard library: + * + * Data type stored by this container. + */ + typedef Number value_type; + + /** + * Declare standard types used in all + * containers. + */ + typedef value_type* pointer; + typedef const value_type* const_pointer; + typedef value_type* iterator; + typedef const value_type* const_iterator; + typedef value_type& reference; + typedef const value_type& const_reference; + typedef size_t size_type; + + + /** + * @name 1: Basic Object-handling + */ + //@{ + /** + * Dummy-Constructor. Dimension=0 + */ + Vector (); + + /** + * Copy-Constructor. Dimension set to that of V ,

+ * all components are copied from V + */ + Vector (const Vector& V); + + /** + * Constructor. Set dimension to #n# and + * initialize all elements with zero. + */ + Vector (const unsigned int n); + + /** + * Destructor. Clears memory + */ + ~Vector (); + + /** + * Set all entries to zero. Equivalent to + * #v = 0#, but more obvious and faster. + * Note that this function does not change + * the size of the vector, unlike the + * STL's #vector<>::clear# function. + */ + void clear (); + + /** + * U(0-N) = s . Fill all components + */ + Vector& operator= (const Number s); + + /** + * U = V . Copy all components + */ + Vector& operator= (const Vector& V); + + /** + * U = U * V . Scalar Produkt + */ + Number operator* (const Vector& V) const; + + /** + * Return square of the l2-norm. + */ + Number norm_sqr () const; + + /** + * Return the mean value of the elements of + * this vector. + */ + Number mean_value () const; + + /** + * Return the l1-norm of the vector, i.e. + * the sum of the absolute values. + */ + Number l1_norm () const; + + /** + * Return the l2-norm of the vector, i.e. + * the square root of the sum of the + * squares of the elements. + */ + Number l2_norm () const; + + /** + * Return the maximum absolute value of the + * elements of this vector. + */ + Number linfty_norm () const; + + + /** + * Change the dimension of the vector to + * #N#. The reserved memory for this vector + * remains unchanged if possible, to make + * things faster, but this may waste some + * memory, so take this in the back of your + * head. + * However, if #N==0# all memory is freed, + * i.e. if you want to resize the vector + * and release the memory not needed, you + * have to first call #reinit(0)# and then + * #reinit(N)#. This cited behaviour is + * analogous to that of the STL containers. + * + * On #fast==false#, the vector is filled by + * zeros. + */ + void reinit (const unsigned int N, const bool fast=false); + + /** + * Change the dimension to that of the + * vector #V#. The same applies as for + * the other #reinit# function. + * + * The elements of #V# are not copied, i.e. + * this function is the same as calling + * #reinit (V.size(), fast)#. + */ + void reinit (const Vector& V, const bool fast=false); + + /** + * Return dimension of the vector. This + * function was formerly called #n()#, but + * was renamed to get the #Vector# class + * closer to the C++ standard library's + * #vector# container. + */ + unsigned int size () const; + + /** + * Return whether the vector contains only + * elements with value zero. This function + * is mainly for internal consistency + * check and should seldomly be used when + * not in debug mode since it uses quite + * some time. + */ + bool all_zero () const; + + /** + * Make the #Vector# class a bit like the + * #vector<># class of the C++ standard + * library by returning iterators to + * the start and end of the elements of this + * vector. + */ + iterator begin (); + + /** + * Return constant iterator to the start of + * the vectors. + */ + const_iterator begin () const; + + /** + * Return an iterator pointing to the + * element past the end of the array. + */ + iterator end (); + + /** + * Return a constant iterator pointing to + * the element past the end of the array. + */ + const_iterator end () const; + //@} + + + /** + * @name 2: Data-Access + */ + //@{ + /** + * Access Components. returns U(i) , + * INLINE + */ + Number operator() (const unsigned int i) const; + + /** + * Access Components. returns U(i) , + * INLINE + */ + Number& operator() (const unsigned int i); + //@} + + + /** + * @name 3: Modification of vectors + */ + //@{ + /** + * Fast equivalent to #U.add(1, V)#. + */ + Vector & operator += (const Vector &V); + + /** + * Fast equivalent to #U.add(-1, V)#. + */ + Vector & operator -= (const Vector &V); + + /** + * U(0-DIM)+=s. + * Addition of #s# to all components. Note + * that #s# is a scalar and not a vector. + */ + void add (const Number s); + + /** + * U+=V. + * Simple vector addition, equal to the + * #operator +=#. + */ + void add (const Vector& V); + + /** + * U+=a*V. + * Simple addition of a scaled vector. + */ + void add (const Number a, const Vector& V); + + /** + * U+=a*V+b*W. + * Multiple addition of scaled vectors. + */ + void add (const Number a, const Vector& V, + const Number b, const Vector& W); + + /** + * U=s*U+V. + * Scaling and simple vector addition. + */ + void sadd (const Number s, const Vector& V); + + /** + * U=s*U+a*V. + * Scaling and simple addition. + */ + void sadd (const Number s, const Number a, const Vector& V); + + /** + * U=s*U+a*V+b*W. + * Scaling and multiple addition. + */ + void sadd (const Number s, const Number a, + const Vector& V, const Number b, const Vector& W); + + /** + * U=s*U+a*V+b*W+c*X. + * Scaling and multiple addition. + */ + void sadd (const Number s, const Number a, + const Vector& V, const Number b, const Vector& W, + const Number c, const Vector& X); + + /** + * Scale each element of the vector by the + * given factor. This function was + * previously called #equ(Number)#, which + * in my eyes is an extremely unintuitive + * naming and was thus replaced. + */ + void scale (const Number factor); + + /** + * U=a*V. Replacing + */ + void equ (const Number a, const Vector& V); + + /** + * U=a*V+b*W. + * Replacing by sum. + */ + void equ (const Number a, const Vector& V, + const Number b, const Vector& W); + + /** + * Compute the elementwise ratio of the + * two given vectors, that is let + * #this[i] = a[i]/b[i]#. This is useful + * for example if you want to compute + * the cellwise ratio of true to estimated + * error. + * + * This vector is appropriately scaled to + * hold the result. + * + * If any of the #b[i]# is zero, the result + * is undefined. No attempt is made to + * catch such situations. + */ + void ratio (const Vector &a, const Vector &b); + //@} + + + /** + * @name 5: Mixed stuff + */ + //@{ + /** + * Output of vector in user-defined format. + */ + void print (FILE* fp, const char* format = 0) const; + + /** + * Output of vector in user-defined format. + */ + void print (const char* format = 0) const; + + /** + * Print to given stream, one element per line. + */ + void print (ostream &) const; + //@} + + /** + * Exception + */ + DeclException2 (ExcDimensionsDontMatch, + int, int, + << "The dimensions " << arg1 << " and " << arg2 + << " do not match here."); + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The given index " << arg1 + << " should be less than " << arg2 << "."); + /** + * Exception + */ + DeclException1 (ExcInvalidNumber, + int, + << "The provided number is invalid here: " << arg1); + /** + * Exception + */ + DeclException0 (ExcOutOfMemory); + /** + * Exception + */ + DeclException0 (ExcEmptyVector); + /** + * Exception + */ + DeclException0 (ExcIO); +}; + + + + + + +/*----------------------- Inline functions ----------------------------------*/ + + +template +inline +unsigned int Vector::size () const +{ + return dim; +} + + + +template +inline +Vector::iterator Vector::begin () { + return &val[0]; +}; + + + +template +inline +Vector::const_iterator Vector::begin () const { + return &val[0]; +}; + + + +template +inline +Vector::iterator Vector::end () { + return &val[dim]; +}; + + + +template +inline +Vector::const_iterator Vector::end () const { + return &val[dim]; +}; + + + +template +inline +Number Vector::operator() (const unsigned int i) const +{ + Assert (i +inline +Number& Vector::operator() (const unsigned int i) +{ + Assert (i +#include +#include + + +template +static inline Number sqr (const Number x) { + return x*x; +}; + + + +template +Vector::Vector () : + dim(0), + maxdim(0), + val(0) +{} + + + +template +Vector::Vector (const unsigned int n) : + dim(0), + maxdim(0), + val(0) +{ + reinit (n, false); +} + + + +template +Vector::Vector (const Vector& v) : + dim(v.size()), + maxdim(v.size()), + val(0) +{ + if (dim) + { + val = new Number[maxdim]; + Assert (val != 0, ExcOutOfMemory()); + copy (v.begin(), v.end(), begin()); + } +} + + + +template +void Vector::reinit (const unsigned int n, const bool fast) { + if (n==0) + { + if (val) delete[] val; + val = 0; + maxdim = dim = 0; + return; + }; + + if (n>maxdim) + { + if (val) delete[] val; + val = new Number[n]; + Assert (val != 0, ExcOutOfMemory()); + maxdim = n; + }; + dim = n; + if (fast == false) + clear (); +} + + + +template +void Vector::reinit (const Vector& v, const bool fast) { + reinit (v.size(), fast); +}; + + + + +template +Vector::~Vector () +{ + if (val) delete[] val; +} + + + +template +void Vector::clear () { + if (dim>0) + fill (begin(), end(), 0.); +} + + + +template +bool Vector::all_zero () const { + Assert (dim!=0, ExcEmptyVector()); + + const_iterator p = begin(), + e = end(); + while (p!=e) + if (*p++ != 0.0) + return false; + return true; +}; + + + +template +Number Vector::operator * (const Vector& v) const +{ + Assert (dim!=0, ExcEmptyVector()); + + if (&v == this) + return norm_sqr(); + + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + + Number sum0 = 0, + sum1 = 0, + sum2 = 0, + sum3 = 0; + + // use modern processors better by + // allowing pipelined commands to be + // executed in parallel + const_iterator ptr = begin(), + vptr = v.begin(), + eptr = ptr + (dim/4)*4; + while (ptr!=eptr) + { + sum0 += (*ptr++ * *vptr++); + sum1 += (*ptr++ * *vptr++); + sum2 += (*ptr++ * *vptr++); + sum3 += (*ptr++ * *vptr++); + }; + // add up remaining elements + while (ptr != end()) + sum0 += *ptr++ * *vptr++; + + return sum0+sum1+sum2+sum3; +} + + + +template +Number Vector::norm_sqr () const +{ + Assert (dim!=0, ExcEmptyVector()); + + Number sum0 = 0, + sum1 = 0, + sum2 = 0, + sum3 = 0; + + // use modern processors better by + // allowing pipelined commands to be + // executed in parallel + const_iterator ptr = begin(), + eptr = ptr + (dim/4)*4; + while (ptr!=eptr) + { + sum0 += sqr(*ptr++); + sum1 += sqr(*ptr++); + sum2 += sqr(*ptr++); + sum3 += sqr(*ptr++); + }; + // add up remaining elements + while (ptr != end()) + sum0 += sqr(*ptr++); + + return sum0+sum1+sum2+sum3; +}; + + + +template +Number Vector::mean_value () const +{ + Assert (dim!=0, ExcEmptyVector()); + + Number sum0 = 0, + sum1 = 0, + sum2 = 0, + sum3 = 0; + + // use modern processors better by + // allowing pipelined commands to be + // executed in parallel + const_iterator ptr = begin(), + eptr = ptr + (dim/4)*4; + while (ptr!=eptr) + { + sum0 += *ptr++; + sum1 += *ptr++; + sum2 += *ptr++; + sum3 += *ptr++; + }; + // add up remaining elements + while (ptr != end()) + sum0 += *ptr++; + + return (sum0+sum1+sum2+sum3)/size(); +}; + + + +template +Number Vector::l1_norm () const +{ + Assert (dim!=0, ExcEmptyVector()); + + Number sum0 = 0, + sum1 = 0, + sum2 = 0, + sum3 = 0; + + // use modern processors better by + // allowing pipelined commands to be + // executed in parallel + const_iterator ptr = begin(), + eptr = ptr + (dim/4)*4; + while (ptr!=eptr) + { + sum0 += fabs(*ptr++); + sum1 += fabs(*ptr++); + sum2 += fabs(*ptr++); + sum3 += fabs(*ptr++); + }; + // add up remaining elements + while (ptr != end()) + sum0 += fabs(*ptr++); + + return sum0+sum1+sum2+sum3; +}; + + + +template +Number Vector::l2_norm () const +{ + return sqrt(norm_sqr()); +}; + + + +template +Number Vector::linfty_norm () const { + Assert (dim!=0, ExcEmptyVector()); + + Number max0=0., + max1=0., + max2=0., + max3=0.; + for (unsigned int i=0; i<(dim/4); ++i) + { + if (max0 +Vector& Vector::operator += (const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + + add (v); + return *this; +} + + + +template +Vector& Vector::operator -= (const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(); + while (i_ptr!=i_end) + *i_ptr++ -= *v_ptr++; + + return *this; +} + + + +template +void Vector::add (const Number v) +{ + Assert (dim!=0, ExcEmptyVector()); + + iterator i_ptr = begin(), + i_end = end(); + while (i_ptr!=i_end) + *i_ptr++ += v; +} + + + +template +void Vector::add (const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(); + while (i_ptr!=i_end) + *i_ptr++ += *v_ptr++; +} + + + +template +void Vector::add (const Number a, const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(); + while (i_ptr!=i_end) + *i_ptr++ += a * *v_ptr++; +} + + + +template +void Vector::add (const Number a, const Vector& v, + const Number b, const Vector& w) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(), + w_ptr = w.begin(); + while (i_ptr!=i_end) + *i_ptr++ += a * *v_ptr++ + b * *w_ptr++; +} + + + +template +void Vector::sadd (const Number x, const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(); + for (; i_ptr!=i_end; ++i_ptr) + *i_ptr = x * *i_ptr + *v_ptr++; +} + + + +template +void Vector::sadd (const Number x, const Number a, const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(); + for (; i_ptr!=i_end; ++i_ptr) + *i_ptr = x * *i_ptr + a * *v_ptr++; +} + + + +template +void Vector::sadd (const Number x, const Number a, + const Vector& v, const Number b, const Vector& w) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(), + w_ptr = w.begin(); + for (; i_ptr!=i_end; ++i_ptr) + *i_ptr = x * *i_ptr + a * *v_ptr++ + b * *w_ptr++; +} + + + +template +void Vector::sadd (const Number x, const Number a, + const Vector& v, const Number b, + const Vector& w, const Number c, const Vector& y) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim)); + Assert (dim == y.dim, ExcDimensionsDontMatch(dim, y.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator v_ptr = v.begin(), + w_ptr = w.begin(), + y_ptr = y.begin(); + + for (; i_ptr!=i_end; ++i_ptr) + *i_ptr = (x * *i_ptr) + (a * *v_ptr++) + (b * *w_ptr++) + (c * *y_ptr++); +} + + + +template +void Vector::scale (const Number factor) +{ + Assert (dim!=0, ExcEmptyVector()); + + iterator ptr=begin(), eptr=end(); + while (ptr!=eptr) + *ptr++ *= factor; +} + + + +template +void Vector::equ (const Number a, const Vector& u, + const Number b, const Vector& v) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim)); + Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator u_ptr = u.begin(), + v_ptr = v.begin(); + while (i_ptr!=i_end) + *i_ptr++ = a * *u_ptr++ + b * *v_ptr++; +} + + + +template +void Vector::equ (const Number a, const Vector& u) +{ + Assert (dim!=0, ExcEmptyVector()); + Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim)); + iterator i_ptr = begin(), + i_end = end(); + const_iterator u_ptr = u.begin(); + while (i_ptr!=i_end) + *i_ptr++ = a * *u_ptr++; +} + + + +template +void Vector::ratio (const Vector &a, const Vector &b) { + Assert (dim!=0, ExcEmptyVector()); + Assert (a.dim == b.dim, ExcDimensionsDontMatch (a.dim, b.dim)); + + // no need to reinit with zeros, since + // we overwrite them anyway + reinit (a.size(), true); + iterator i_ptr = begin(), + i_end = end(); + const_iterator a_ptr = a.begin(), + b_ptr = b.begin(); + while (i_ptr!=i_end) + *i_ptr++ = *a_ptr++ / *b_ptr++; +}; + + + +template +Vector& Vector::operator = (const Number s) +{ + Assert (dim!=0, ExcEmptyVector()); + fill (begin(), end(), s); + return *this; +} + + + +template +Vector& Vector::operator = (const Vector& v) +{ + if (v.dim != dim) + reinit (v.dim, true); + + if (dim!=0) + copy (v.begin(), v.end(), begin()); + return *this; +} + + + +template +void Vector::print (FILE* f, const char* format) const +{ + Assert (dim!=0, ExcEmptyVector()); + if (!format) format = " %5.2f"; + for (unsigned int j=0;j +void Vector::print (const char* format) const +{ + Assert (dim!=0, ExcEmptyVector()); + if (!format) format = " %5.2f"; + for (unsigned int j=0;j +void Vector::print (ostream &out) const { + Assert (dim!=0, ExcEmptyVector()); + for (unsigned int i=0; i + +// explicit instantiations +template class Vector; +template class Vector; -- 2.39.5