vector.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_VECTOR
#define MFEM_VECTOR

// Data type vector

#include "../general/array.hpp"
#include <cmath>
#include <iostream>
#if defined(_MSC_VER) && (_MSC_VER < 1800)
#include <float.h>
#define isfinite _finite
#endif

namespace mfem
{

inline int CheckFinite(const double *v, const int n);

class Vector
{
protected:

   int size, allocsize;
   double * data;

public:

   Vector () { allocsize = size = 0; data = 0; }

   Vector(const Vector &);

   explicit Vector (int s);

   Vector (double *_data, int _size)
   { data = _data; size = _size; allocsize = -size; }

   void Load (std::istream ** in, int np, int * dim);

   void Load(std::istream &in, int Size);

   void Load(std::istream &in) { int s; in >> s; Load (in, s); }

   void SetSize(int s);

   void SetData(double *d) { data = d; }

   void SetDataAndSize(double *d, int s)
   { data = d; size = s; allocsize = -s; }

   void NewDataAndSize(double *d, int s)
   {
      if (allocsize > 0) { delete [] data; }
      SetDataAndSize(d, s);
   }

   void MakeDataOwner() { allocsize = abs(allocsize); }

   void Destroy();

   inline int Size() const { return size; }

   // double *GetData() { return data; }

   inline double *GetData() const { return data; }

   inline operator double *() { return data; }

   inline operator const double *() const { return data; }

   inline bool OwnsData() const { return (allocsize > 0); }

   inline void StealData(double **p)
   { *p = data; data = 0; size = allocsize = 0; }

   inline double *StealData() { double *p; StealData(&p); return p; }

   double & Elem (int i);

   const double & Elem (int i) const;

   inline double & operator() (int i);

   inline const double & operator() (int i) const;

   double operator*(const double *) const;

   double operator*(const Vector &v) const;

   Vector & operator=(const double *v);

   Vector & operator=(const Vector &v);

   Vector & operator=(double value);

   Vector & operator*=(double c);

   Vector & operator/=(double c);

   Vector & operator-=(double c);

   Vector & operator-=(const Vector &v);

   Vector & operator+=(const Vector &v);

   Vector & Add(const double a, const Vector &Va);

   Vector & Set(const double a, const Vector &x);

   void SetVector (const Vector &v, int offset);

   void Neg();

   inline void Swap(Vector &other);

   friend void add(const Vector &v1, const Vector &v2, Vector &v);

   friend void add(const Vector &v1, double alpha, const Vector &v2, Vector &v);

   friend void add(const double a, const Vector &x, const Vector &y, Vector &z);

   friend void add (const double a, const Vector &x,
                    const double b, const Vector &y, Vector &z);

   friend void subtract(const Vector &v1, const Vector &v2, Vector &v);

   friend void subtract(const double a, const Vector &x,
                        const Vector &y, Vector &z);

   void median(const Vector &lo, const Vector &hi);

   void GetSubVector(const Array<int> &dofs, Vector &elemvect) const;
   void GetSubVector(const Array<int> &dofs, double *elem_data) const;

   void SetSubVector(const Array<int> &dofs, const Vector &elemvect);
   void SetSubVector(const Array<int> &dofs, double *elem_data);

   void AddElementVector(const Array<int> & dofs, const Vector & elemvect);
   void AddElementVector(const Array<int> & dofs, double *elem_data);
   void AddElementVector(const Array<int> & dofs, const double a,
                         const Vector & elemvect);

   void SetSubVectorComplement(const Array<int> &dofs, const double val);

   void Print(std::ostream & out = std::cout, int width = 8) const;

   void Print_HYPRE(std::ostream &out) const;

   void Randomize(int seed = 0);
   double Norml2() const;
   double Normlinf() const;
   double Norml1() const;
   double Normlp(double p) const;
   double Max() const;
   double Min() const;
   double Sum() const;
   double DistanceTo (const double *p) const;

   int CheckFinite() const { return mfem::CheckFinite(data, size); }

   virtual ~Vector ();
};

// Inline methods

inline int CheckFinite(const double *v, const int n)
{
   // isfinite didn't appear in a standard until C99, and later C++11
   // It wasn't standard in C89 or C++98.  PGI as of 14.7 still defines
   // it as a macro, which sort of screws up everybody else.
   int bad = 0;
   for (int i = 0; i < n; i++)
   {
#ifdef isfinite
      if (!isfinite(v[i]))
#else
      if (!std::isfinite(v[i]))
#endif
      {
         bad++;
      }
   }
   return bad;
}

inline Vector::Vector (int s)
{
   if (s > 0)
   {
      allocsize = size = s;
      data = new double[s];
   }
   else
   {
      allocsize = size = 0;
      data = NULL;
   }
}

inline void Vector::SetSize(int s)
{
   if (s == size)
   {
      return;
   }
   if (s <= abs(allocsize))
   {
      size = s;
      return;
   }
   if (allocsize > 0)
   {
      delete [] data;
   }
   allocsize = size = s;
   data = new double[s];
}

inline void Vector::Destroy()
{
   if (allocsize > 0)
   {
      delete [] data;
   }
   allocsize = size = 0;
   data = NULL;
}

inline double & Vector::operator() (int i)
{
   MFEM_ASSERT(data && i >= 0 && i < size,
               "index [" << i << "] is out of range [0," << size << ")");

   return data[i];
}

inline const double & Vector::operator() (int i) const
{
   MFEM_ASSERT(data && i >= 0 && i < size,
               "index [" << i << "] is out of range [0," << size << ")");

   return data[i];
}

inline void Vector::Swap(Vector &other)
{
   mfem::Swap(size, other.size);
   mfem::Swap(allocsize, other.allocsize);
   mfem::Swap(data, other.data);
}

template<> inline void Swap<Vector>(Vector &a, Vector &b)
{
   a.Swap(b);
}

inline Vector::~Vector()
{
   if (allocsize > 0)
   {
      delete [] data;
   }
}

inline double Distance(const double *x, const double *y, const int n)
{
   using namespace std;
   double d = 0.0;

   for (int i = 0; i < n; i++)
   {
      d += (x[i]-y[i])*(x[i]-y[i]);
   }

   return sqrt(d);
}

}

#endif