/*
 * RealSpaceMatrix.hpp
 *
 *  Created on: Jun 25, 2010
 *      Author: crueger
 */

#ifndef REALSPACEMATRIX_HPP_
#define REALSPACEMATRIX_HPP_

#include <iosfwd>
#include "defs.hpp"

/**
 * Simple class to store matrices with a few extra functions
 */

class Vector;
class MatrixContent;

/** 3x3 Matrix class.
 * This class has some specific features for 3D space such as rotations or
 * hard-coded determinants.
 */
class RealSpaceMatrix
{
  friend Vector operator*(const RealSpaceMatrix&,const Vector&);
public:
  RealSpaceMatrix();

  /**
   * construct a matrix from a 3x3 double array that contains all Elements.
   *
   * The elements are laid out in the following way:
   * array -> matrix
   * 0 -> (0,0)
   * 1 -> (1,0)
   * 2 -> (2,0)
   * 3 -> (0,1)
   * 4 -> (1,1)
   * 5 -> (2,1)
   * 6 -> (0,2)
   * 7 -> (1,2)
   * 8 -> (2,2)
   *
   */
  RealSpaceMatrix(const double*);
  RealSpaceMatrix(const RealSpaceMatrix&);
  RealSpaceMatrix(const MatrixContent&);
  virtual ~RealSpaceMatrix();

  /**
   * Set the matrix to a unit matrix.
   */
  void setIdentity();

  /**
   * Set all matrix entries to zero.
   */
  void setZero();

  /**
   * Sets all matrix corresponding to a rotation matrix,
   * first around the x-, then the y-, finally the z-axis
   * with given angles.
   */
  void setRotation(const double x, const double y, const double z);

  /**
   * Access the matrix at index (i,j)
   */
  double &at(size_t i, size_t j);
  /**
   * Access the matrix at index (i,j)
   */
  const double at(size_t i, size_t j) const;

  /**
   * Set the matrix at index (i,j).
   *
   * Slightly faster than at(i,j)=x
   */
  void set(size_t i, size_t j, const double value);

  /**
   * get the ith row of the matrix as a vector
   */
  Vector &row(size_t);
  const Vector &row(size_t) const;

  /**
   * get the ith column of the matrix as a vector
   */
  Vector &column(size_t);
  const Vector &column(size_t) const;

  /**
   * get the diagonal of the matrix as a vector
   */
  Vector &diagonal();
  const Vector &diagonal() const;

  /**
   * Calculate the determinant of the matrix
   */
  double determinant() const;

  /**
   * Calculate the inverse of the matrix.
   *
   * Rather costly, so use precomputation as often as possible.
   */
  RealSpaceMatrix invert() const;

  /**
   * Diagonalizes a matrix and sets its rows to the resulting eigenvalues.
   * The eigenvalues are returned as a vector.
   *
   * Rather costly, so use precomputation as often as possible.
   */
  Vector transformToEigenbasis();

  /**
   * Calculate the transpose of the matrix.
   */
  RealSpaceMatrix transpose() const;
  RealSpaceMatrix &transpose();

  // operators
  RealSpaceMatrix &operator=(const RealSpaceMatrix&);

  const RealSpaceMatrix &operator+=(const RealSpaceMatrix&);
  const RealSpaceMatrix &operator-=(const RealSpaceMatrix&);
  const RealSpaceMatrix &operator*=(const RealSpaceMatrix&);

  const RealSpaceMatrix &operator*=(const double);

  const RealSpaceMatrix operator+(const RealSpaceMatrix&) const;
  const RealSpaceMatrix operator-(const RealSpaceMatrix&) const;
  const RealSpaceMatrix operator*(const RealSpaceMatrix&) const;

  bool operator==(const RealSpaceMatrix&) const;

private:
  RealSpaceMatrix(MatrixContent*);
  void createViews();
  MatrixContent *content;
  // we keep around some Vector views of the matrix, to return references
  Vector* rows_ptr[NDIM];
  Vector* columns_ptr[NDIM];
  Vector* diagonal_ptr;
};

const RealSpaceMatrix operator*(const double,const RealSpaceMatrix&);
const RealSpaceMatrix operator*(const RealSpaceMatrix&,const double);

/**
 * Takes a symmetric matrix that stores the lower diagonal and produces a
 * full matrix.
 *
 * The array is laid out as follows:
 *
 * array -> matrix
 * 0 -> (0,0)
 * 1 -> (1,0);(0,1)
 * 2 -> (1,1)
 * 3 -> (2,0);(0,2)
 * 4 -> (2,1);(1,2)
 * 5 -> (2,2)
 */
RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const cell_size);

std::ostream &operator<<(std::ostream&,const RealSpaceMatrix&);
Vector operator*(const RealSpaceMatrix&,const Vector&);
Vector& operator*=(Vector&,const RealSpaceMatrix&);

#endif /* REALSPACEMATRIX_HPP_ */