--- /dev/null
+/*---------------------------- 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 <base/exceptions.h>
+#include <cstdio>
+
+
+
+/**
+ * Vector of data.
+ * Memory for Components is supplied explicitly <p>
+ * ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! ) <p>
+ * - all necessary methods for Vectors are supplied <p>
+ * - operators available are `=` , `*` and `( )` <p>
+ * CONVENTIONS for used `equations` : <p>
+ * - THIS vector is always named `U` <p>
+ * - 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 <typename Number>
+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 , <p>
+ * all components are copied from V
+ */
+ Vector (const Vector<Number>& 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<Number>& operator= (const Number s);
+
+ /**
+ * U = V . Copy all components
+ */
+ Vector<Number>& operator= (const Vector<Number>& V);
+
+ /**
+ * U = U * V . Scalar Produkt
+ */
+ Number operator* (const Vector<Number>& 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<Number>& 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<Number> & operator += (const Vector<Number> &V);
+
+ /**
+ * Fast equivalent to #U.add(-1, V)#.
+ */
+ Vector<Number> & operator -= (const Vector<Number> &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<Number>& V);
+
+ /**
+ * U+=a*V.
+ * Simple addition of a scaled vector.
+ */
+ void add (const Number a, const Vector<Number>& V);
+
+ /**
+ * U+=a*V+b*W.
+ * Multiple addition of scaled vectors.
+ */
+ void add (const Number a, const Vector<Number>& V,
+ const Number b, const Vector<Number>& W);
+
+ /**
+ * U=s*U+V.
+ * Scaling and simple vector addition.
+ */
+ void sadd (const Number s, const Vector<Number>& V);
+
+ /**
+ * U=s*U+a*V.
+ * Scaling and simple addition.
+ */
+ void sadd (const Number s, const Number a, const Vector<Number>& V);
+
+ /**
+ * U=s*U+a*V+b*W.
+ * Scaling and multiple addition.
+ */
+ void sadd (const Number s, const Number a,
+ const Vector<Number>& V, const Number b, const Vector<Number>& W);
+
+ /**
+ * U=s*U+a*V+b*W+c*X.
+ * Scaling and multiple addition.
+ */
+ void sadd (const Number s, const Number a,
+ const Vector<Number>& V, const Number b, const Vector<Number>& W,
+ const Number c, const Vector<Number>& 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<Number>& V);
+
+ /**
+ * U=a*V+b*W.
+ * Replacing by sum.
+ */
+ void equ (const Number a, const Vector<Number>& V,
+ const Number b, const Vector<Number>& 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<Number> &a, const Vector<Number> &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 <typename Number>
+inline
+unsigned int Vector<Number>::size () const
+{
+ return dim;
+}
+
+
+
+template <typename Number>
+inline
+Vector<Number>::iterator Vector<Number>::begin () {
+ return &val[0];
+};
+
+
+
+template <typename Number>
+inline
+Vector<Number>::const_iterator Vector<Number>::begin () const {
+ return &val[0];
+};
+
+
+
+template <typename Number>
+inline
+Vector<Number>::iterator Vector<Number>::end () {
+ return &val[dim];
+};
+
+
+
+template <typename Number>
+inline
+Vector<Number>::const_iterator Vector<Number>::end () const {
+ return &val[dim];
+};
+
+
+
+template <typename Number>
+inline
+Number Vector<Number>::operator() (const unsigned int i) const
+{
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ return val[i];
+}
+
+
+
+template <typename Number>
+inline
+Number& Vector<Number>::operator() (const unsigned int i)
+{
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ return val[i];
+}
+
+
+
+
+
+/*---------------------------- vector.h ---------------------------*/
+/* end of #ifndef __vector_H */
+#endif
+/*---------------------------- vector.h ---------------------------*/
--- /dev/null
+// $Id$
+
+// 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, 1998, 1999
+
+#include <lac/vector.h>
+#include <cmath>
+#include <algorithm>
+
+
+template <typename Number>
+static inline Number sqr (const Number x) {
+ return x*x;
+};
+
+
+
+template <typename Number>
+Vector<Number>::Vector () :
+ dim(0),
+ maxdim(0),
+ val(0)
+{}
+
+
+
+template <typename Number>
+Vector<Number>::Vector (const unsigned int n) :
+ dim(0),
+ maxdim(0),
+ val(0)
+{
+ reinit (n, false);
+}
+
+
+
+template <typename Number>
+Vector<Number>::Vector (const Vector<Number>& 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 <typename Number>
+void Vector<Number>::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 <typename Number>
+void Vector<Number>::reinit (const Vector<Number>& v, const bool fast) {
+ reinit (v.size(), fast);
+};
+
+
+
+
+template <typename Number>
+Vector<Number>::~Vector ()
+{
+ if (val) delete[] val;
+}
+
+
+
+template <typename Number>
+void Vector<Number>::clear () {
+ if (dim>0)
+ fill (begin(), end(), 0.);
+}
+
+
+
+template <typename Number>
+bool Vector<Number>::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 <typename Number>
+Number Vector<Number>::operator * (const Vector<Number>& 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 <typename Number>
+Number Vector<Number>::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 <typename Number>
+Number Vector<Number>::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 <typename Number>
+Number Vector<Number>::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 <typename Number>
+Number Vector<Number>::l2_norm () const
+{
+ return sqrt(norm_sqr());
+};
+
+
+
+template <typename Number>
+Number Vector<Number>::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<fabs(val[4*i])) max0=fabs(val[4*i]);
+ if (max1<fabs(val[4*i+1])) max1=fabs(val[4*i+1]);
+ if (max2<fabs(val[4*i+2])) max2=fabs(val[4*i+2]);
+ if (max3<fabs(val[4*i+3])) max3=fabs(val[4*i+3]);
+ };
+ // add up remaining elements
+ for (unsigned int i=(dim/4)*4; i<dim; ++i)
+ if (max0<fabs(val[i]))
+ max0 = fabs(val[i]);
+
+ return max (max(max0, max1),
+ max(max2, max3));
+};
+
+
+
+
+
+template <typename Number>
+Vector<Number>& Vector<Number>::operator += (const Vector<Number>& v)
+{
+ Assert (dim!=0, ExcEmptyVector());
+
+ add (v);
+ return *this;
+}
+
+
+
+template <typename Number>
+Vector<Number>& Vector<Number>::operator -= (const Vector<Number>& 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 <typename Number>
+void Vector<Number>::add (const Number v)
+{
+ Assert (dim!=0, ExcEmptyVector());
+
+ iterator i_ptr = begin(),
+ i_end = end();
+ while (i_ptr!=i_end)
+ *i_ptr++ += v;
+}
+
+
+
+template <typename Number>
+void Vector<Number>::add (const Vector<Number>& 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 <typename Number>
+void Vector<Number>::add (const Number a, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::add (const Number a, const Vector<Number>& v,
+ const Number b, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::sadd (const Number x, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::sadd (const Number x, const Number a, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::sadd (const Number x, const Number a,
+ const Vector<Number>& v, const Number b, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::sadd (const Number x, const Number a,
+ const Vector<Number>& v, const Number b,
+ const Vector<Number>& w, const Number c, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::scale (const Number factor)
+{
+ Assert (dim!=0, ExcEmptyVector());
+
+ iterator ptr=begin(), eptr=end();
+ while (ptr!=eptr)
+ *ptr++ *= factor;
+}
+
+
+
+template <typename Number>
+void Vector<Number>::equ (const Number a, const Vector<Number>& u,
+ const Number b, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::equ (const Number a, const Vector<Number>& 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 <typename Number>
+void Vector<Number>::ratio (const Vector<Number> &a, const Vector<Number> &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 <typename Number>
+Vector<Number>& Vector<Number>::operator = (const Number s)
+{
+ Assert (dim!=0, ExcEmptyVector());
+ fill (begin(), end(), s);
+ return *this;
+}
+
+
+
+template <typename Number>
+Vector<Number>& Vector<Number>::operator = (const Vector<Number>& v)
+{
+ if (v.dim != dim)
+ reinit (v.dim, true);
+
+ if (dim!=0)
+ copy (v.begin(), v.end(), begin());
+ return *this;
+}
+
+
+
+template <typename Number>
+void Vector<Number>::print (FILE* f, const char* format) const
+{
+ Assert (dim!=0, ExcEmptyVector());
+ if (!format) format = " %5.2f";
+ for (unsigned int j=0;j<size();j++)
+ fprintf(f, format, val[j]);
+ fputc('\n',f);
+}
+
+
+
+template <typename Number>
+void Vector<Number>::print (const char* format) const
+{
+ Assert (dim!=0, ExcEmptyVector());
+ if (!format) format = " %5.2f";
+ for (unsigned int j=0;j<size();j++)
+ printf (format, val[j]);
+ printf ("\n");
+}
+
+
+
+template <typename Number>
+void Vector<Number>::print (ostream &out) const {
+ Assert (dim!=0, ExcEmptyVector());
+ for (unsigned int i=0; i<size(); ++i)
+ out << val[i] << endl;
+
+ AssertThrow (out, ExcIO());
+};
+
+
+