vector.hh

 
#ifndef vector_hh
#define vector_hh
 
#include <cmath>
#include <string>
 
#include "basics/typedefs.hh"
#include "basis.hh"
#include "space/spaceSet.hh"
#include "toolbox/array.hh"
 
namespace concepts {
 
  // forward declaration
  template<class F, class G>
  class LinearForm;
 
  template<class F>
  class Space;
 
  // **************************************************************** Vector **
 
  template<class F>
  class Vector : public Function<F> {
  public:
    Vector(const Vector<F>& f);
    template<class H>
    Vector(const Vector<H>& f);
 
    template<class G>
    Vector(const Space<G>& spc) 
      : Function<F>(spc), data_(this->dim_, 0), v_((F*)data_) {}
 
    Vector(const uint dim=0)
      : Function<F>(dim), data_(this->dim_, 0), v_((F*)data_)  {}
 
    template<class G>
    Vector(const Space<G>& spc, const LinearForm<F,G>& lf);
 
    template<class G>
    Vector(const Space<G>& spc, const std::string& fname);
 
    template<class G>
    Vector(const Space<G>& spc, F* data) :
      Function<F>(spc), data_(spc.dim()), v_((F*)data_) {
      std::memcpy(v_, data, this->dim_ * sizeof(v_[0]));
    }
    Vector(uint dim, F* data) :
      Function<F>(dim), data_(dim), v_((F*)data_) {
      std::memcpy(v_, data, this->dim_ * sizeof(v_[0]));
    }
 
    template<class H>
    Vector(const Vector<H>& fncX, F fnc(const H&));
    template<class H>
    Vector(const Vector<H>& fncX, const F& fnc(const H&));
 
    Vector(const Vector<typename Realtype<F>::type>& V_R,
           const Vector<typename Realtype<F>::type>& V_I);
    
    template<class H>
    Vector(const Vector<H>& fnc, const Set<IndexRange>& indices);
 
    virtual ~Vector();
 
    virtual Vector<F>& operator=(const Function<F>& fnc);
 
    Vector<F>& operator=(const Vector<F>& fnc);
 
    /* Assignement operator from different number type
 
       This is extra, because virtual member templates don't exist.
 
       The compiler checks, if the assigment is possible, e.g.
         Cmplx = Real, possible
         Real = Cmplx, not possible
     */
    template<class H>
    Vector<F>& operator=(const Function<H>& fnc);
 
    Vector<F>& operator=(F c);
 
    virtual F& operator()(uint i) {
      conceptsAssert3(i < this->dim_, Assertion(), "i = " << i << ", n = " << this->dim_);
      return v_[i]; }
    virtual F operator()(uint i) const {
      conceptsAssert3(i < this->dim_, Assertion(), "i = " << i << ", n = " << this->dim_);
      return v_[i]; }
 
    uint size() const{
      return n();
    }
 
    virtual Function<F>& operator+=(const Function<F>& fnc);
    virtual Function<F>& operator+=(F c);
    Vector<F> operator+ (const Function<F>& fnc) const;
    Vector<F> operator+ (F c) const;
 
    virtual Function<F>& operator-=(const Function<F>& fnc);
    virtual Function<F>& operator-=(F c);
    Vector<F> operator- (const Function<F>& fnc) const;
    Vector<F> operator- (F c) const;
 
    virtual Function<F>& operator*=(F sc);
    Vector<F> operator* (F c) const;
    virtual Function<F>& operator/=(F sc);
    Vector<F> operator/ (F c) const;
 
    //    template <class T>
    //    friend std::ostream& operator<<(std::ostream & os, const Vector<T> & v);
 
    F operator*(const Vector<F>& fnc) const;
 
    template<class G>
    Vector<F>& assemble(const Space<G>& spc, const LinearForm<F,G>& lf);
 
    Vector<F>& apply(F fnc(const F&));
 
    operator F*() const { return v_; }
    
    F* data() const { return v_; }
 
    virtual void resize(uint n) {
      this->dim_ = n;
      data_.resize(n);
      v_ = (F*)data_;
    }
 
    virtual Function<F>& add(const Function<F>& fnc, F sc);
 
    template<class H>
    Vector<F>& add(const Vector<H>& fnc, F sc, uint offset = 0);
 
    template<class H>
    Vector<F>& add(const Vector<H>& fnc) { return add(fnc, F(1.0)); }
 
    Real l1() const;
 
    inline Real l2() const { return sqrt(l2_2()); }
 
    Real l2_2() const;
 
    Real max() const;
 
    uint n() const { return this->dim_; }
 
    // reverse the vector
    void reverse();
 
    const Vector<F>& write(const std::string& fname) const;
 
    void storeMatlab(const char* filename, const char* name = 0,
                     bool append = false) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
 
    template<class G>
    void fillEntries_(const Space<G>& spc, const LinearForm<F,G>& lf);
    Array<F> data_;
  private:
    F* v_;
    static uint storeMatlabCounter_;
  };
 
  template<class F>
  template<class H>
  Vector<F>::Vector(const Vector<H>& f) : 
    Function<F>(f.dim()), data_(this->dim_), v_((F*)data_) {
    memorycpy(v_, (H*)f, this->dim_);
  }
 
  template<class F>
  template<class H>
  Vector<F>::Vector(const Vector<H>& fncX, F fnc(const H&)) 
    : Function<F>(fncX.dim()), data_(this->dim_, 0), v_((F*)data_)
  {
    F* v = v_;
    H* w = (H*)fncX;
    for(uint i = this->dim_; i; --i) *v++= fnc(*w++);
  }
 
  template<class F>
  template<class G>
  Vector<F>::Vector(const Space<G>& spc, const LinearForm<F,G>& lf) :
    Function<F>(spc), data_(this->dim_), v_((F*)data_) {
    fillEntries_(spc, lf);
  }
 
  template<class F>
  template<class H>
  Vector<F>::Vector(const Vector<H>& fncX, const F& fnc(const H&))
    : Function<F>(fncX.dim()), data_(this->dim_, 0), v_((F*)data_)
  {
    F* v = v_;
    H* w = (H*)fncX;
    for(uint i = this->dim_; i; --i) *v++= fnc(*w++);
  }
 
  template<class F>
  template<class H>
  Vector<F>::Vector(const Vector<H>& fnc, const Set<IndexRange>& indices) :
    Function<F>(indices.dim()), data_(this->dim_), v_((F*)data_)
  {
    F* d = data_;
    for (Set<IndexRange>::index_iterator i = indices.index_begin();
         i != indices.index_end(); )
      *d++ = fnc(*i++);
  }
 
  template<class F>
  template<class H>
  Vector<F>& Vector<F>::operator=(const Function<H>& fnc)
  {
    for(uint i = 0; i < this->dim_; ++i) v_[i] = fnc(i);
    return *this;
  }
 
 
  template<class F>
  template<class H>
  Vector<F>& Vector<F>::add(const Vector<H>& fnc, F sc, uint offset)
  {
    F* d = v_ + offset;
    const H* s = (H*)fnc;
    for (uint i = std::min(this->dim_-offset,fnc.dim()); i; --i) *d++ += sc * *s++;
    return *this;
  }
 
  template<class F>
  void Vector<F>::reverse(){
  F *p1 = (F*)data_, *p2 = (F*)data_ + (data_.size()-1);
    for(uint i = std::floor(double(data_.size())/2.0); i--; )
      std::swap(*p1++, *p2--);
  }
 
} // namespace concepts
 
#endif // vector_hh