]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add templated data vectors. They will replace our old dVector eventually.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Feb 1999 13:54:57 +0000 (13:54 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Feb 1999 13:54:57 +0000 (13:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@878 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/vector.h [new file with mode: 0644]
deal.II/lac/include/lac/vector.templates.h [new file with mode: 0644]
deal.II/lac/source/vector.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h
new file mode 100644 (file)
index 0000000..7ae125e
--- /dev/null
@@ -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 <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     ---------------------------*/
diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h
new file mode 100644 (file)
index 0000000..b05f8cf
--- /dev/null
@@ -0,0 +1,556 @@
+// $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());
+};
+
+
+
diff --git a/deal.II/lac/source/vector.cc b/deal.II/lac/source/vector.cc
new file mode 100644 (file)
index 0000000..f825a49
--- /dev/null
@@ -0,0 +1,7 @@
+// $Id$
+
+#include <lac/vector.templates.h>
+
+// explicit instantiations
+template class Vector<double>;
+template class Vector<float>;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.