function/vector.hh

Go to the documentation of this file.

/* Function derivation: Vector, stores the values in an array
 */

#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) 
      : Function<F>(dim), data_(this->dim_, 0), v_((F*)data_)  {}

    template<class G>
    Vector(const Space<G>& spc, 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 Function<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>
    Function<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]; }

    int size() {
      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;

    F operator*(const Vector<F>& fnc) const;

    template<class G>
    Vector<F>& assemble(const Space<G>& spc, LinearForm<F,G>& lf);

    Vector<F>& apply(F fnc(const F&));

    operator F*() const { return v_; }
    
    F* data() const { return v_; }

    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_; }

    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, 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, 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>
  Function<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;
  }

} // namespace concepts

#endif // vector_hh