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

#ifndef MATRIX_HPP_
#define MATRIX_HPP_

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

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

class Vector;
class MatrixContent;

class Matrix
{
  friend Vector operator*(const Matrix&,const Vector&);
public:
  Matrix();

  /**
   * 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)
   *
   */
  Matrix(const double*);
  Matrix(const Matrix&);
  virtual ~Matrix();

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

  /**
   * 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.
   */
  Matrix invert() const;

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

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

  Matrix &operator*=(const double);

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

  bool operator==(const Matrix&) const;

private:
  Matrix(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;
};

Matrix operator*(const double,const Matrix&);
Matrix operator*(const Matrix&,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)
 */
Matrix ReturnFullMatrixforSymmetric(const double * const cell_size);

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

#endif /* MATRIX_HPP_ */