// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file MatrixVector.h */

#ifndef AMDIS_MATRIXVECTOR_H
#define AMDIS_MATRIXVECTOR_H

#include <vector>
#include "Global.h"
#include "MemoryManager.h"
#include "Serializable.h"

namespace AMDiS {

  template<typename T> class DOFVector;

  // ============================================================================
  // ===== class Vector =========================================================
  // ============================================================================

  /** \brief
   * Class for efficient vector operations of fixed size numerical vectors.
   */
  template<typename T>
  class Vector : public Serializable
  {
  public:
    MEMORY_MANAGED(Vector<T>);

    /** \brief
     * Constructor.
     */
    Vector(int i = 0) 
      : size(i) 
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
    };

    inline bool used() const
    {
      return (valArray != NULL);
    };

    /** \brief
     * Change the size of the vector to newSize.
     */
    inline void resize(int newSize) {
      if (size != newSize) {
	if (valArray) 
	  delete [] valArray;
	valArray = new T[newSize];
	size = newSize;
      }
    };

    /** \brief
     * Copy constructor.
     */
    Vector(const Vector<T>& rhs) 
      : Serializable(),size(rhs.size)
    {
      valArray = new T[rhs.size];
      *this = rhs; // uses operator=()
    };

    /** \brief
     * Destructor.
     */
    virtual ~Vector() { 
      delete [] valArray; 
    };

    /** \brief
     * Assignement operator
     */
    inline const Vector<T>& operator=(const Vector<T>& rhs) {
      TEST_EXIT_DBG(rhs.size == size)("invalid size\n");
      T *rhsIt, *thisIt;
      for(rhsIt = rhs.begin(), thisIt = this->begin();
	  rhsIt != rhs.end();
	  ++rhsIt, ++thisIt)
	{
	  *thisIt = *rhsIt;
	}
      return *this;
    };

    /** \brief
     * Assignement operator
     */
    inline const Vector<T>& operator=(const T& scal) {
      T *thisIt;
      for(thisIt = this->begin();
	  thisIt != this->end();
	  ++thisIt)
	{
	  *thisIt = scal;
	}
      return *this;
    };

    /** \brief
     * Assignement operator
     */
    inline const Vector<T>& operator=(const T* vec) {
      T *thisIt;
      const T *vecIt;
      for(thisIt = this->begin(), vecIt = &vec[0];
	  thisIt != this->end();
	  ++thisIt, ++vecIt)
	{
	  *thisIt = *vecIt;
	}
      return *this;
    };

    /** \brief
     * Sets all entries to scal.
     */
    inline const Vector<T>& set(const T& scal) {
      return *this = scal;
    };

    /** \brief
     * Sets all entries.
     */
    inline const Vector<T>& setValues(const T* values) {
      T *thisIt;
      const T *valuesIt;
      for(thisIt = this->begin(), valuesIt = values; 
	  thisIt != this->end(); 
	  ++thisIt, ++valuesIt) 
	{
	  *thisIt = *valuesIt;
	}
      return *this;
    };

    /** \brief
     * Sets all entries.
     */
    inline void fill(const T value) {
      for (T *thisIt = this->begin(); thisIt != this->end(); thisIt++) {
	*thisIt = value;
      }
    }

    /** \brief
     * Comparison operator.
     */
    inline bool operator==(const Vector<T>& rhs) const {
      if(size != rhs.size) return false;
      T *rhsIt, *thisIt;
      for(rhsIt = rhs.begin(), thisIt = this->begin();
	  rhsIt != rhs.end();
	  ++rhsIt, ++thisIt)
	{
	  if(*thisIt != *rhsIt) return false;
	}
      return true;
    };

    /** \brief
     * Comparison operator.
     */
    inline bool operator!=(const Vector<T>& rhs) const {
      return !(*this==rhs);
    }

    /** \brief
     * Access to the i-th vector element.
     */
    inline T& operator[](int i) {
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
      return valArray[i];
    };

    /** \brief
     * Access to the i-th vector element for const vectors.
     */
    inline const T& operator[] (int i) const {
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
      return valArray[i];
    };

    /** \brief
     * Returns pointer to the first vector element.
     */
    inline T *begin() const { 
      return valArray; 
    };

    /** \brief
     * Returns pointer after the last vector element.
     */
    inline T *end() const { 
      return valArray + size; 
    };

    /** \brief
     * Returns \ref size.
     */
    virtual int getSize() const { 
      return size; 
    };

    /** \brief
     * Returns \ref valArray as T-array
     */
    inline T *getValArray() { 
      return valArray; 
    };

    void print() const {
      std::cout << this->size << " vector: " << std::endl;
      for (int i = 0; i < size; i++) {
	std::cout << this->valArray[i] << " ";
      }
      std::cout << std::endl;
    };

    // ===== Serializable implementation =====
  
    void serialize(std::ostream &out) {
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
    };

    void deserialize(std::istream &in) {
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
    };

    std::string getTypeName() const { 
      return "Vector"; 
    };

  protected:
    /** \brief
     * Size of the vector.
     */
    int size;

    /** \brief
     * Internal storage pointer.
     */
    T *valArray;
  };



  // ============================================================================
  // ===== class Matrix =========================================================
  // ============================================================================

  /** \brief
   * Class for efficient matrix operations of fixed size numerical matrices.
   */
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
    MEMORY_MANAGED(Matrix<T>);

    /** \brief
     * Constructor.
     */
    Matrix(int r, int c)
      : Vector<T>(r*c), 
	rows(r),
	cols(c)
    {};

    /** \brief
     * Changes the size of the matrix to newRows x newCols.
     */
    inline void resize(int newRows, int newCols) {
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
    };

    /** \brief
     * Assignement operator.
     */
    inline const Matrix<T>& operator=(const T& scal) {
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
    };

    inline bool operator==(const Matrix<T>& rhs) const {
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
    };

    /** \brief
     * Comparison operator.
     */
    inline bool operator!=(const Matrix<T>& rhs) const {
      return !(*this == rhs);
    }

    /** \brief
     * Acces to i-th matrix row.
     */
    inline T *operator[](int i) {
      return this->valArray + cols * i;
    };

    /** \brief
     * Acces to i-th matrix row for constant matrices.
     */
    inline const T *operator[](int i) const {
      return this->valArray + cols * i;
    };

    /** \brief
     * Returns \ref rows.
     */
    inline int getNumRows() const { 
      return rows; 
    };

    /** \brief
     * Return \ref cols.
     */
    inline int getNumCols() const { 
      return cols; 
    };

    /** \brief
     * Returns \ref rows.
     */
    inline int getSize() const {
      return rows; 
    }; 

    /** \brief
     * Returns pointer after the last vector element.
     */
    inline T *end() const { 
      return this->valArray + (cols * rows); 
    };

    void print() const {
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
      for (int i = 0; i < rows; i++) {
	for (int j = 0; j < cols; j++) {
	  std::cout << this->valArray[i * cols + j] << " ";
	}
	std::cout << std::endl;
      }
    };

  protected:
    /** \brief
     * Number of matrix rows.
     */
    int rows;

    /** \brief
     * Number of matrix columns.
     */
    int cols;
  };

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
	 ++resultIt) {
      *resultIt = 0;
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
	*resultIt += *mIt * *vIt;
      }
    }
    return result;
  };

  /** \brief
   * Computes x(Ay)^T, with x and y vectors, and A a matrix.
   */
  template<typename T>
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
  {
    TEST_EXIT_DBG(a.getNumCols() == x.getSize())("A and x not compatible\n");
    TEST_EXIT_DBG(a.getNumCols() == y.getSize())("A and y not compatible\n");

    T result = 0;
    T tmp = 0;
    T *aIt, *xIt, *yIt;

    for (aIt = a.begin(), xIt = x.begin(); aIt != a.end(); ++xIt) {
      tmp = 0;
      for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt) {
	tmp += *aIt * *yIt;
      }
      result += *xIt * tmp;
    }

    return result;
  }

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) {
    return mv(m, v, v);
  };

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) {
    Vector<T> result(v.getSize());
    return mv(m, v, result);
  };


  /** \brief
   * Scalar product.
   */
  template<typename T> 
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
  {
    double result = 0.0;

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It) {
      result += *v1It * *v2It;
    }

    return result;
  };

  /** \brief
   * Vector addition.
   */
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
    TEST_EXIT_DBG(v1.getSize() == v2.getSize())("invalid size in test v1 == v2\n");
    TEST_EXIT_DBG(v2.getSize() == result.getSize())("invalid size in test v2 == result\n");
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
      {
	*resultIt = *v1It + *v2It;
      }

    return result;
  };

  /** \brief
   * scalar * vector
   */
  template<typename T>
  inline const Vector<T>& mult(const T& scal, 
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");

    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt) {
      *resultIt = scal * *vIt;
    }

    return result;
  };

  /** \brief
   * vector + scalar
   */
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt)
      {
	*resultIt = *vIt + scal;
      }

    return result;
  };

  /** \brief
   * y = a * x + y.
   */
  template<typename T>
  inline const Vector<T>& axpy(const T& a,
			       const Vector<T> &x,
			       Vector<T> &y)
  {
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
	 ++xIt, ++yIt) {
      *yIt += a * *xIt;
    }

    return y;
  };

  template<typename T>
  inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) {
    return mult(scal, v, v);
  };

  template<typename T>
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) {
    Vector<T> result = v;
    result *= scal;
    return result;
  };

  template<typename T>
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) {
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
  };

  template<typename T>
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) {
    return add(v1, v2, v1);
  };

  template<typename T>
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) {
    Vector<T> result = v1;
    result += v2;
    return result;
  };

  template<typename T>
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2){
    return axpy(-1.0, v2, v1);
  };

  template<typename T>
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2){
    Vector<T> result = v1;
    result -= v2;
    return result;
  };

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
    for (vIt = v->begin(); vIt != v->end(); ++vIt) {
      result += *vIt * *vIt;
    }
    return sqrt(result);
  };

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
    z[0] = x[1] * y[2] - x[2] * y[1];
    z[1] = x[2] * y[0] - x[0] * y[2];
    z[2] = x[0] * y[1] - x[1] * y[0];
  };

}

#endif // AMDIS_MATRIXVECTOR_H