source: ThirdParty/LinearAlgebra/src/LinearAlgebra/VectorContent.hpp@ 7a434f

/*
 * VectorContent.hpp
 *
 *  Created on: Jul 2, 2010
 *      Author: crueger
 */

#ifndef VECTORCONTENT_HPP_
#define VECTORCONTENT_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/Assert.hpp"

#include <boost/serialization/array.hpp>
#include <boost/serialization/split_member.hpp>

/**
 * !file
 * The way GSL works does not allow for forward definitions of the structures.
 * Because of this the pointer to the gsl_vector struct is wrapped inside another
 * (dumb) object that allows for forward definitions.
 *
 * DO NOT USE OUTSIDE OF VECTOR UNLESS YOU ABSOLUTELY HAVE TO AND KNOW WHAT YOU ARE DOING.
 */

#include <iosfwd>

#include <gsl/gsl_vector.h>

#include "MatrixVector_ops.hpp"

class Vector;
class MatrixContent;
struct VectorViewContent;
namespace boost {
  namespace serialization {
    class access;
  }
}

/** This class safely stores vector dimension away.
 *
 * This way we simulate const-ness of the dimension while still allowing
 * serialization to modify them.
 *
 */
struct VectorDimension {
  friend struct VectorViewContent;
public:
  VectorDimension(size_t _dimension) : dimension(_dimension) {}

  size_t getDimension() const {return dimension;}

private:
  void setDim(size_t _dimension) {dimension = _dimension;}

  friend class boost::serialization::access;
  // serialization (due to gsl_vector separate versions needed)
  template<class Archive>
  void serialize(Archive & ar, const unsigned int version)
  {
    ar & dimension;
  }

private:
  size_t dimension;      // store for internal purposes
};

/** Dummy structure to create a unique signature.
 *
 */
struct VectorBaseCase{};

class VectorContent : public VectorDimension
{
  friend std::ostream & operator<< (std::ostream& ost, const VectorContent &m);
  friend VectorContent const operator*(const VectorContent& a, const double m);
  friend VectorContent const operator*(const double m, const VectorContent& a);
  friend VectorContent const operator+(const VectorContent& a, const VectorContent& b);
  friend VectorContent const operator-(const VectorContent& a, const VectorContent& b);

  // matrix vector products
  friend VectorContent const operator*(const VectorContent& vec, const MatrixContent& mat);
  friend VectorContent const operator*(const MatrixContent& mat, const VectorContent& vec);

public:
  explicit VectorContent(size_t _dim);
  VectorContent(VectorBaseCase);
  VectorContent(const VectorContent * const src);
  VectorContent(const VectorContent & src);
  VectorContent(gsl_vector * _src, bool _free_content_on_exit = false);
  VectorContent(const size_t _dimension, std::istream &inputstream);
  virtual ~VectorContent();

  static size_t preparseVectorDimensions(std::istream &inputstream);

  // Accessing
  double &at(size_t m);
  const double at(size_t m) const;
  double & operator[](size_t i);
  const double operator[](size_t i) const;
  double *Pointer(size_t m) const;
  const double *const_Pointer(size_t m) const;
  void write(std::ostream &ost) const;

  // Assignment operator
  VectorContent &operator=(const VectorContent& src);

  // Initializing
  void setFromDoubleArray(double * x);
  void setFromVector(Vector &v);
  void setValue(double x);
  void setZero();
  int setBasis(size_t m);

  // Exchanging elements
  int SwapElements(size_t i, size_t j);
  int Reverse();

  // checking state
  bool IsZero() const;
  bool IsOne() const;

  // properties
  //bool IsNormalTo(const VectorContent &normal) const;
  //bool IsEqualTo(const VectorContent &a) const;

  // Norms
  double Norm() const;
  double NormSquared() const;
  void Normalize();
  VectorContent getNormalized() const;

  // properties relative to another VectorContent
  double DistanceSquared(const VectorContent &y) const;
  double ScalarProduct(const VectorContent &y) const;
  double Angle(const VectorContent &y) const;

  // operators
  bool operator==(const VectorContent& b) const;
  const VectorContent& operator+=(const VectorContent& b);
  const VectorContent& operator-=(const VectorContent& b);
  const VectorContent& operator*=(const double m);
  const double operator*(const VectorContent& b) const;

  bool operator!=(const VectorContent &other) const {
    return !(*this == other);
  }

private:
  VectorContent();
  friend class boost::serialization::access;
  // serialization (due to gsl_vector separate versions needed)
  template<class Archive>
  void save(Archive & ar, const unsigned int version) const
  {
    ar & boost::serialization::base_object<VectorDimension>(*this);
    ar & content->size;
    ar & content->stride;
    ar & boost::serialization::make_array<double>(content->data, content->size);
  }
  template<class Archive>
  void load(Archive & ar, const unsigned int version)
  {
    ar & boost::serialization::base_object<VectorDimension>(*this);
    content = gsl_vector_calloc(getDimension());
    ar & content->size;
    ASSERT(getDimension() == content->size,
        "MatrixContent::load() - dimension mismatch: "
        +toString(getDimension())+" != "+toString(content->size)+".");
    ar & content->stride;
    ar & boost::serialization::make_array<double>(content->data, content->size);
    // this is always a copy, hence always free
    free_content_on_exit = true;
  }
  BOOST_SERIALIZATION_SPLIT_MEMBER()

public:
  gsl_vector *content;
private:
  bool free_content_on_exit;
};
//BOOST_CLASS_EXPORT_GUID(VectorContent, "VectorContent")


std::ostream & operator << (std::ostream& ost, const VectorContent &m);
VectorContent const operator*(const VectorContent& a, const double m);
VectorContent const operator*(const double m, const VectorContent& a);
VectorContent const operator+(const VectorContent& a, const VectorContent& b);
VectorContent const operator-(const VectorContent& a, const VectorContent& b);

/** Vector view.
 * Extension of VectorContent to contain not a gsl_vector but only a view on a
 * gsl_vector (or row/column in a gsl_matrix).
 *
 * We need the above VectorBaseCase here:
 * content, i.e. the gsl_vector, must not be allocated, as it is just a view
 * with an internal gsl_vector_view. Hence, VectorBaseCase is just dummy class
 * to give the constructor a unique signature.
 */
struct VectorViewContent : public VectorContent
{
  VectorViewContent(gsl_vector_view _view) :
    VectorContent(VectorBaseCase()),
    view(_view)
  {
    dimension = _view.vector.size;
    content=&view.vector;
  }
  ~VectorViewContent(){
    content=0;
  }
  gsl_vector_view view;
};

#endif /* VECTORCONTENT_HPP_ */