#ifndef __MATRIX_GSL_H
#define __MATRIX_GSL_H

/** @file
 * @brief 行列ライブラリのGSLラッパー
 *
 */

#include <cstring>
#include <ostream>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>

#include "param/complex_gsl.h"
#include "param/matrix.h"
#include "util/util.h"

template <class T, class GSLType>
struct gsl_matrix_func_list_t{
  static GSLType *(*alloc)(size_t, size_t);
  static GSLType *(*calloc)(size_t, size_t);
  static void (*free)(GSLType *); // a null pointer is not allowed. <= マニュアル(8.4.1)より 
  static T *(*ptr)(GSLType *, size_t, size_t);
  static void (*set_zero)(GSLType *);
  static int (*memcpy)(GSLType *, const GSLType *);
};

#define MAKE_FUNC_LIST_INSTANCE(gsl_t, gsl_mat_t) \
template <> \
gsl_mat_t *(*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::alloc)(size_t, size_t) = gsl_mat_t ## _alloc; \
template <> \
gsl_mat_t *(*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::calloc)(size_t, size_t) = gsl_mat_t ## _calloc; \
template <> \
void (*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::free)(gsl_mat_t *) = gsl_mat_t ## _free; \
template <> \
gsl_t *(*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::ptr)(gsl_mat_t *, size_t, size_t) = gsl_mat_t ## _ptr; \
template <> \
void (*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::set_zero)(gsl_mat_t *) = gsl_mat_t ## _set_zero; \
template <> \
int (*gsl_matrix_func_list_t<gsl_t, gsl_mat_t>::memcpy)(gsl_mat_t *, const gsl_mat_t *) = gsl_mat_t ## _memcpy;

MAKE_FUNC_LIST_INSTANCE(double, gsl_matrix);
MAKE_FUNC_LIST_INSTANCE(gsl_complex, gsl_matrix_complex);

#undef MAKE_FUNC_LIST_INSTANCE

template <class T, class GSLType, class GSLMatrixType>
class Array2D_GSL : public Array2D<T> {
  protected:
    typedef T elm_t;
    typedef GSLType gsl_t;
    typedef GSLMatrixType gsl_mat_t;

  private:
    typedef Array2D<elm_t> base_t;
    typedef Array2D_GSL<elm_t, gsl_t, gsl_mat_t> self_t;

  protected:
    static const gsl_matrix_func_list_t<gsl_t, gsl_mat_t> func_list;
#define gsl_func(name) (func_list.name)
    gsl_mat_t *m_gsl_mat;  ///< 確保したメモリ
    int *ref;   ///< 参照カウンタ

    /**
     * Array2D_Denseクラスのコンストラクタ。
     * 指定の行サイズ、指定の列サイズで2次元配列を生成します。
     *
     * @param rows 行数
     * @param columns 列数
     * @throw StorageException メモリが確保できなかった場合
     */
    Array2D_GSL(
        gsl_mat_t *mat) throw(StorageException)
          : base_t(mat->size1, mat->size2),
            m_gsl_mat(mat),
            ref(new int(0)) {
      // どうせbad_alloc例外がでるので、メモリ確保失敗を捕捉しないことにした
      //if(!ref){throw StorageException("Lack of memory!!");}
      (*ref)++;
    }

  public:
    using base_t::rows;
    using base_t::columns;
    
    /**
     * 単純配列化したものを返します。
     *
     * @return (T *) 単純配列
     */
    gsl_t *buffer() const {return m_gsl_mat->data;}

    /**
     * 内部のバッファを返します。
     *
     * @return (GSLType *) 内部のバッファ
     */
    gsl_mat_t *raw() const {return m_gsl_mat;}

    /**
     * 指定した行列成分を返します。
     *
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     * @throw StorageException 参照インデックスが不正の場合
     */
    elm_t &operator()(
        const unsigned int &row,
        const unsigned int &column) throw(StorageException){
      gsl_t *res(gsl_func(ptr)(m_gsl_mat, row, column));
      if(!res){throw StorageException("Index incorrect");}
      return *(elm_t *)res;
    }

    /**
     * 要素のゼロクリアをします。
     *
     */
    void clear(){
      gsl_func(set_zero)(m_gsl_mat);
    }

    /**
     * Array2D_Denseクラスのコンストラクタ。
     * 指定の行サイズ、指定の列サイズで2次元配列を生成します。
     *
     * @param rows 行数
     * @param columns 列数
     * @throw StorageException メモリが確保できなかった場合
     */
    Array2D_GSL(
        const unsigned int &rows,
        const unsigned int &columns) throw(StorageException)
          : base_t(rows, columns),
            m_gsl_mat(NULL),
            ref(new int(0)) {
      // どうせbad_alloc例外がでるので、refについてメモリ確保失敗を捕捉しないことにした
      if(((rows * columns > 0) && !(m_gsl_mat = gsl_func(alloc)(rows, columns)))
          /*|| !ref*/){
        // サイズゼロの行列を生成しようとすると怒られる
        throw StorageException("Lack of memory!!");
      }
      (*ref)++;
    }

    /**
     * Array2D_Denseクラスのコンストラクタ。
     * 指定の行サイズ、指定の列サイズで2次元配列を生成します。
     * また成分はserilaizedによって指定された値で生成されます
     *
     * @param rows 行数
     * @param columns 列数
     * @param serialized 成分
     * @throw StorageException メモリが確保できなかった場合
     */
    Array2D_GSL(
        const unsigned int &rows,
        const unsigned int &columns,
        const gsl_t *serialized)
      throw(StorageException)
        : base_t(rows, columns),
          m_gsl_mat(NULL), ref(new int(0)) {
      // どうせbad_alloc例外がでるので、refについてメモリ確保失敗を捕捉しないことにした
      if(((rows * columns > 0) && !(m_gsl_mat = gsl_func(alloc)(rows, columns)))
          /*|| !ref*/){
        // サイズゼロの行列を生成しようとすると怒られる
        throw StorageException("Lack of memory!!");
      }  
      (*ref)++;
      memcpy(m_gsl_mat->data, serialized, sizeof(gsl_t) * rows * columns);
    }

    /**
     * コピーコンストラクタ
     *
     * @param array コピー元
     */
    Array2D_GSL(const self_t &array)
        : base_t(array.m_rows, array.m_columns){
      if(m_gsl_mat = array.m_gsl_mat){(*(ref = array.ref))++;}
    }
    /**
     * デストラクタ。
     * 参照カウンタを減算すると共に、もしカウンタが0の場合、
     * 確保したメモリをdeleteによって解放します。
     */
    ~Array2D_GSL(){
      if(ref && ((--(*ref)) <= 0)){
        if(m_gsl_mat){gsl_func(free)(m_gsl_mat);}
        delete ref;
      }
    }

    /**
     * 代入演算子。
     * 高速なシャローコピーを行います。
     *
     * @return (Array2D_Dense<T>) 自分自身
     */
    self_t &operator=(const self_t &array){
      if(this != &array){
        if(ref && ((--(*ref)) <= 0)){
          if(m_gsl_mat){gsl_func(free)(m_gsl_mat);}
          delete ref;
        }
        if(m_gsl_mat = array.m_gsl_mat){
          base_t::m_rows = array.m_rows;
          base_t::m_columns = array.m_columns;
          (*(ref = array.ref))++;
        }
      }
      return *this;
    }
    
    struct AllElements {
      typedef AllElements iterator;
      self_t &array2d;
      unsigned int r, c;
      AllElements(
          self_t &_array2d,
          const unsigned int _row = 0, const unsigned int _column = 0) 
          : array2d(_array2d),
          r(_row), c(_column) {}
      AllElements(
          const AllElements &orig) 
          : array2d(orig.array2d),
          r(orig.r), c(orig.c) {}
      ~AllElements(){}
      AllElements begin() const {return AllElements(array2d, 0, 0);}
      AllElements end() const {return AllElements(array2d, array2d.rows(), 0);}
      AllElements &operator++(){
        if(++c >= array2d.columns()){
          c = 0;
          ++r;
        }
        return *this;
      }
      bool operator!=(const AllElements &another){
        return (r != another.r) || (c != another.c);
      }
      T &operator*(){return array2d(r, c);}
      unsigned int row() {return r;}
      unsigned int column() {return c;}
    };
    
    using base_t::iterate_operation;

    virtual void all_elements(void (*op)(T &t)){
      iterate_operation(AllElements(*this), op);
    }
    
    virtual void all_elements(typename base_t::IterateOperator &op){
      iterate_operation(AllElements(*this), op);
    }
    
    virtual void all_elements(typename base_t::IterateOperator2 &op){
      iterate_operation(AllElements(*this), op);
    }
#undef gsl_func
};

#define MAKE_ARRAY2D_DENSE_INSTANCE(T, GSLType, GSLMatrixType) \
template <> \
class Array2D_Dense<T > \
    : public Array2D_GSL<T, GSLType, GSLMatrixType> { \
  private: \
    typedef Array2D_GSL<T, GSLType, GSLMatrixType> super_t; \
    typedef Array2D<super_t::elm_t> base_t; \
    typedef Array2D_Dense<super_t::elm_t> self_t; \
\
  protected: \
    Array2D_Dense( \
        super_t::gsl_mat_t *mat) throw(StorageException) \
          : super_t(mat) {} \
\
  public: \
    Array2D_Dense( \
        const unsigned int &rows, \
        const unsigned int &columns) throw(StorageException) \
          : super_t(rows, columns){} \
\
    Array2D_Dense( \
        const unsigned int &rows, \
        const unsigned int &columns, \
        const super_t::gsl_t *serialized) \
      throw(StorageException) \
        : super_t(rows, columns, serialized){} \
\
    Array2D_Dense(const self_t &array) \
        : super_t(array){} \
\
    base_t *shallow_copy() const {return new self_t(*this);} \
\
    base_t *copy() const throw(StorageException){ \
      self_t *array( \
          new self_t(rows(), columns())); \
      if(super_t::m_gsl_mat){ \
        super_t::func_list.memcpy(array->m_gsl_mat, super_t::m_gsl_mat); \
      } \
      return array; \
    } \
\
    self_t dense() const {return self_t(*this);} \
\
    self_t *unchained_instance() const { \
      if((super_t::ref) && (*super_t::ref > 1)){ \
        return static_cast<self_t *>(copy()); \
      }else{ \
        return static_cast<self_t *>(shallow_copy()); \
      } \
    } \
\
    ~Array2D_Dense(){} \
\
    static self_t *calloc( \
        const unsigned int &rows, \
        const unsigned int &columns){ \
      return (rows * columns > 0) \
          ? new self_t(super_t::func_list.calloc(rows, columns)) \
          : new self_t(rows, columns); \
    } \
}

MAKE_ARRAY2D_DENSE_INSTANCE(double, double, gsl_matrix);
MAKE_ARRAY2D_DENSE_INSTANCE(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_ARRAY2D_DENSE_INSTANCE

/**
 * 転置行列で密2次元配列
 *
 * @return (Array2D_Dense<T> *)
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template <> \
Array2D_Dense<T > Array2D_Transpose<T >::dense() const { \
  Array2D_Dense<T > array(rows(), columns()); \
  GSLMatrixType ## _transpose_memcpy(array.raw(), getTarget().dense().raw()); \
  return array; \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * Matrixクラスのコンストラクタ
 * 指定の行数、指定の列数で行列を生成
 * 成分はすべてT(0)で初期化
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template <> \
Matrix<T >::Matrix( \
      const unsigned int &rows, \
      const unsigned int &columns \
    ) throw(MatrixException) \
    : m_Storage(Array2D_Dense<T >::calloc(rows, columns)){}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 単位行列の生成
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template <> \
Matrix<T > Matrix<T >::getI(const unsigned int &size){ \
  Array2D_Dense<T > *array2d(new Array2D_Dense<T >(size, size)); \
  GSLMatrixType ## _set_identity(array2d->raw()); \
  return Matrix<T >(array2d); \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 行列の成分全てを指定倍
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template <> \
Matrix<T > Matrix<T >::operator*(const T &scalar) const{ \
  Array2D_Dense<T > *dense(storage()->dense().unchained_instance()); \
  GSLMatrixType ## _scale( \
      dense->raw(), \
      (const GSLType &)scalar); \
  return Matrix<T >(dense); \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 行列の成分全てを除算
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template<> \
Matrix<T > Matrix<T >::operator/(const T &scalar) const{ \
  return (*this) * (1. / scalar); \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 行列を成分ごとに加算
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template<> \
Matrix<T > Matrix<T >::operator+(const Matrix<T > &matrix) const{ \
  Array2D_Dense<T > *dense(storage()->dense().unchained_instance()); \
  GSLMatrixType ## _add( \
      dense->raw(), \
      matrix.storage()->dense().raw()); \
  return Matrix<T >(dense); \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 行列を成分ごとに減算
 */
#define MAKE_SPECIALIZED(T, GSLType, GSLMatrixType) \
template<> \
Matrix<T > Matrix<T >::operator-(const Matrix<T > &matrix) const{ \
  Array2D_Dense<T > *dense(storage()->dense().unchained_instance()); \
  GSLMatrixType ## _sub( \
      dense->raw(), \
      matrix.storage()->dense().raw()); \
  return Matrix<T >(dense); \
}

MAKE_SPECIALIZED(double, double, gsl_matrix);
MAKE_SPECIALIZED(Complex<double>, gsl_complex, gsl_matrix_complex);

#undef MAKE_SPECIALIZED

/*
 * 行列を乗算
 */
template<>
Matrix<double> Matrix<double>::operator*(const Matrix<double> &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), matrix.columns()));
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1,
      storage()->dense().raw(), matrix.storage()->dense().raw(), 0,
      array2d->raw());
  return Matrix<double>(array2d);
}

/*
 * 行列を乗算(非転置 * 転置)
 */
template<>
Matrix<double> Matrix<double>::operator*(
    const TransposedMatrix<double> &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), matrix.columns()));
  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1,
      storage()->dense().raw(), 
      matrix.untranspose().storage()->dense().raw(), 0,
      array2d->raw());
  return Matrix<double>(array2d);
}

/*
 * 行列を乗算(転置 * 非転置)
 */
template<>
Matrix<double> TransposedMatrix<double>::operator*(
    const Matrix<double> &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), matrix.columns()));
  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1,
      untranspose().storage()->dense().raw(), 
      matrix.storage()->dense().raw(), 0,
      array2d->raw());
  return Matrix<double>::make_instance(array2d);
}

/*
 * 行列を乗算(転置 * 転置)
 */
template<>
Matrix<double> TransposedMatrix<double>::operator*(
    const TransposedMatrix<double> &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), matrix.columns()));
  gsl_blas_dgemm(CblasTrans, CblasTrans, 1,
      untranspose().storage()->dense().raw(), 
      matrix.untranspose().storage()->dense().raw(), 0,
      array2d->raw());
  return Matrix<double>::make_instance(array2d);
}

/*
 * 行列を乗算
 */
template<>
Matrix<Complex<double> > Matrix<Complex<double> >::operator*(
    const Matrix<Complex<double> > &matrix)
      const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), matrix.columns()));
  gsl_complex unity, zero;
  GSL_REAL(unity) = 1;  GSL_IMAG(unity) = 0;
  GSL_REAL(zero) = 0;   GSL_IMAG(zero) = 0;
  gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, unity,
      storage()->dense().raw(), matrix.storage()->dense().raw(), zero,
      array2d->raw());
  return Matrix<Complex<double> >(array2d);
}

/*
 * 行列を乗算(非転置 * 転置)
 */
template<>
Matrix<Complex<double> > Matrix<Complex<double> >::operator*(
    const TransposedMatrix<Complex<double> > &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), matrix.columns()));
  gsl_complex unity, zero;
  GSL_REAL(unity) = 1;  GSL_IMAG(unity) = 0;
  GSL_REAL(zero) = 0;   GSL_IMAG(zero) = 0;
  gsl_blas_zgemm(CblasNoTrans, CblasTrans, unity,
      storage()->dense().raw(), 
      matrix.untranspose().storage()->dense().raw(), zero,
      array2d->raw());
  return Matrix<Complex<double> >(array2d);
}

/*
 * 行列を乗算(転置 * 非転置)
 */
template<>
Matrix<Complex<double> > TransposedMatrix<Complex<double> >::operator*(
    const Matrix<Complex<double> > &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), matrix.columns()));
  gsl_complex unity, zero;
  GSL_REAL(unity) = 1;  GSL_IMAG(unity) = 0;
  GSL_REAL(zero) = 0;   GSL_IMAG(zero) = 0;
  gsl_blas_zgemm(CblasTrans, CblasNoTrans, unity,
      untranspose().storage()->dense().raw(), 
      matrix.storage()->dense().raw(), zero,
      array2d->raw());
  return Matrix<Complex<double> >::make_instance(array2d);
}

/*
 * 行列を乗算(転置 * 転置)
 */
template<>
Matrix<Complex<double> > TransposedMatrix<Complex<double> >::operator*(
    const TransposedMatrix<Complex<double> > &matrix)
    const throw(MatrixException){
  if(columns() != matrix.rows()){
    throw MatrixException("Operation void!!");
  }
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), matrix.columns()));
  gsl_complex unity, zero;
  GSL_REAL(unity) = 1;  GSL_IMAG(unity) = 0;
  GSL_REAL(zero) = 0;   GSL_IMAG(zero) = 0;
  gsl_blas_zgemm(CblasTrans, CblasTrans, unity,
      untranspose().storage()->dense().raw(), 
      matrix.untranspose().storage()->dense().raw(), zero,
      array2d->raw());
  return Matrix<Complex<double> >::make_instance(array2d);
}

/**
 * 単項演算子-
 */
template<>
Matrix<double> Matrix<double>::operator-() const{
  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());
  gsl_matrix_scale(
      dense->raw(), 
      -1);
  return Matrix<double>(dense);
}

/**
 * 単項演算子-
 */
template<>
Matrix<Complex<double> > Matrix<Complex<double> >::operator-() const{
  Array2D_Dense<Complex<double> > *dense(storage()->dense().unchained_instance());
  gsl_complex neg_unity;
  GSL_REAL(neg_unity) = -1;  GSL_IMAG(neg_unity) = 0;
  gsl_matrix_complex_scale(
      dense->raw(), 
      neg_unity);
  return Matrix<Complex<double> >(dense);
}

/**
 * 行列式
 */
template<>
double Matrix<double>::determinant(bool do_check)
    const throw(MatrixException){
  if(do_check && !isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());
  int sign;
  gsl_linalg_LU_decomp(dense->raw(), p, &sign);
  double det(gsl_linalg_LU_det(dense->raw(), sign));
  gsl_permutation_free(p);
  delete dense;
  return det;
}

/**
 * 行列式
 */
template<>
Complex<double> Matrix<Complex<double> >::determinant(bool do_check)
    const throw(MatrixException){
  if(do_check && !isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<Complex<double> > *dense(storage()->dense().unchained_instance());
  int sign;
  gsl_linalg_complex_LU_decomp(dense->raw(), p, &sign);
  gsl_complex det(gsl_linalg_complex_LU_det(dense->raw(), sign));
  gsl_permutation_free(p);
  delete dense;
  return (Complex<double>)det;
}

/**
 * 逆行列
 */
template <>
Matrix<double> Matrix<double>::inverse() const throw(MatrixException){
  if(!isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());
  int sign;
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), rows()));
  gsl_linalg_LU_decomp(dense->raw(), p, &sign);
  gsl_linalg_LU_invert(dense->raw(), p, array2d->raw());
  gsl_permutation_free(p);
  delete dense;
  return Matrix<double>(array2d);
}

/**
 * 逆行列
 */
template <>
Matrix<Complex<double> > Matrix<Complex<double> >::inverse()
    const throw(MatrixException){
  if(!isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<Complex<double> > *dense(storage()->dense().unchained_instance());
  int sign;
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), rows()));
  gsl_linalg_complex_LU_decomp(dense->raw(), p, &sign);
  gsl_linalg_complex_LU_invert(dense->raw(), p, array2d->raw());
  gsl_permutation_free(p);
  delete dense;
  return Matrix<Complex<double> >(array2d);
}

/*
 * LU分解
 */
template<>
Matrix<double> Matrix<double>::decomposeLU(bool do_check)
    const throw(MatrixException){
  if(do_check && !isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());
  gsl_matrix *lu_gsl(dense->raw());
  int sign;
  gsl_linalg_LU_decomp(lu_gsl, p, &sign);
  Array2D_Dense<double> *array2d(
      new Array2D_Dense<double>(rows(), rows() * 2));
  gsl_matrix *res(array2d->raw());
#define l(i, j) *gsl_matrix_ptr(res, i, j)
#define u(i, j) *gsl_matrix_ptr(res, i, j + rows())
#define lu(i, j) gsl_matrix_get(lu_gsl, i, j)
  for(unsigned i(0); i < rows(); i++){
    l(i, i) = 1;
    u(i, i) = lu(i, i);
    for(unsigned j(i + 1); j < rows(); j++){
      l(i, j) = 0;
      l(j, i) = lu(j, i);
      u(i, j) = lu(i, j);
      u(j, i) = 0;
    }
  }
#undef l
#undef u
#undef lu
  gsl_permutation_free(p);
  delete dense;
  return Matrix<double>(array2d);
}

/*
 * LU分解
 */
template<>
Matrix<Complex<double> > Matrix<Complex<double> >::decomposeLU(bool do_check)
    const throw(MatrixException){
  if(do_check && !isSquare()){throw MatrixException("Operation void!!");}
  gsl_permutation *p(gsl_permutation_calloc(rows()));
  Array2D_Dense<Complex<double> > *dense(storage()->dense().unchained_instance());
  gsl_matrix_complex *lu_gsl(dense->raw());
  int sign;
  gsl_linalg_complex_LU_decomp(lu_gsl, p, &sign);
  Array2D_Dense<Complex<double> > *array2d(
      new Array2D_Dense<Complex<double> >(rows(), rows() * 2));
  gsl_matrix_complex *res(array2d->raw());
  gsl_complex unity, zero;
  GSL_REAL(unity) = 1;  GSL_IMAG(unity) = 0;
  GSL_REAL(zero) = 0;   GSL_IMAG(zero) = 0;
#define l(i, j) *gsl_matrix_complex_ptr(res, i, j)
#define u(i, j) *gsl_matrix_complex_ptr(res, i, j + rows())
#define lu(i, j) gsl_matrix_complex_get(lu_gsl, i, j)
  for(unsigned i(0); i < rows(); i++){
    l(i, i) = unity;
    u(i, i) = lu(i, i);
    for(unsigned j(i + 1); j < rows(); j++){
      l(i, j) = zero;
      l(j, i) = lu(j, i);
      u(i, j) = lu(i, j);
      u(j, i) = zero;
    }
  }
#undef l
#undef u
#undef lu
  gsl_permutation_free(p);
  delete dense;
  return Matrix<Complex<double> >(array2d);
}

/*
 * ヘッセンベルク行列
 */
template<>
Matrix<double> Matrix<double>::hessenberg(Matrix<double> *transform)
    const throw(MatrixException) {
  if(!isSquare()){throw MatrixException("Operation void!!");}
  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());

  gsl_vector *tau(gsl_vector_alloc(rows()));
  gsl_linalg_hessenberg_decomp(
      dense->raw(), 
      tau);

  if(transform){
    gsl_linalg_hessenberg_unpack_accum(
        dense->raw(),
        tau, transform->storage()->dense().raw());
  }

  gsl_linalg_hessenberg_set_zero(
      dense->raw());
  gsl_vector_free(tau);

  return Matrix<double>(dense);
}

/*
 * 固有値、固有ベクトル
 */
template <>
Matrix<Complex<double> > Matrix<double>::eigen(
    const double &threshold_abs,
    const double &threshold_rel) const throw(MatrixException){

  if(!isSquare()){throw MatrixException("Operation void!!");}

  Array2D_Dense<double> *dense(storage()->dense().unchained_instance());

  gsl_vector_complex *eigen_val(gsl_vector_complex_alloc(rows()));
  gsl_matrix_complex *eigen_vec(gsl_matrix_complex_alloc(rows(), rows()));
  gsl_eigen_nonsymmv_workspace *w(gsl_eigen_nonsymmv_alloc(rows()));
  gsl_eigen_nonsymmv(dense->raw(), eigen_val, eigen_vec, w);
  delete dense;

  //結果の格納用の行列
  Matrix<Complex<double> > result(rows(), rows() + 1);
  for(unsigned j(0); j < rows(); j++){
    result(j, rows()) = gsl_vector_complex_get(eigen_val, j);
    for(unsigned i(0); i < rows(); i++){
      result(i, j) = gsl_matrix_complex_get(eigen_vec, i, j);
    }
  }

  gsl_eigen_nonsymmv_free(w);
  gsl_vector_complex_free(eigen_val);
  gsl_matrix_complex_free(eigen_vec);

  return result;
}

/*
 * 固有値、固有ベクトル
 */
template <>
Matrix<Complex<double> > Matrix<double>::eigen()
    const throw(MatrixException){
  return eigen(0, 0);
}

#if 0
/**
 * @brief 行列
 *
 * 行列を定義したクラス。
 * 様々な行列の演算を定義しています。
 *
 * なお、内部的に参照カウンタを利用したライトウエイトな実装になっているため、
 * メモリや演算回数が節約されることが機体されます。
 * そのために明示的にcopy()メソッドを使用しないとディープコピーがされないので、
 * 再利用を行う際は注意してください。
 *
 * @see Array2D 内部的に利用する2次元配列を定義したクラス。
 */
template <>
class Matrix<double>{

    /**
     * LU行列であること利用して線型方程式(Ax = y)のxを解きます。
     *
     * @param y 右辺
     * @param do_check LU分解済み行列の定義を満たしているか、確認する
     * @return (Matrix<T2>) x(解)
     */
    template <class T2>
    Matrix<T2> solve_linear_eq_with_LU(
        const Matrix<T2> &y, bool do_check = true)
        const throw(MatrixException) {
      bool not_LU(false);
      if(do_check){
        if(rows() * 2 != columns()){not_LU = true;}
      }
      self_t L(partial(rows(), rows(), 0, 0)),
             U(partial(rows(), rows(), 0, rows()));
      if(do_check){
        for(unsigned i(1); i < rows(); i++){
          for(unsigned j(0); j < i; j++){
            if(U(i, j) != elm_t(0)){not_LU = true;}
            if(L(j, i) != elm_t(0)){not_LU = true;}
          }
        }
      }
      if(not_LU){
        throw MatrixException("Not LU decomposed matrix!!");
      }
      if((y.columns() != 1)
          || (y.rows() != rows())){
        throw MatrixException("Operation void!!");
      }


      // L(Ux) = y で y' = (Ux)をまず解く
      Matrix<T2> y_copy(y.copy());
      Matrix<T2> y_prime(Matrix<T2>::naked(y.rows(), 1));
      for(unsigned i(0); i < rows(); i++){
        y_prime(i, 0) = y_copy(i, 0) / L(i, i);
        for(unsigned j(i + 1); j < rows(); j++){
          y_copy(j, 0) -= L(j, i) * y_prime(i, 0);
        }
      }

      // 続いてUx = y'で xを解く
      Matrix<T2> x(Matrix<T2>::naked(y.rows(), 1));
      for(unsigned i(rows()); i > 0;){
        i--;
        x(i, 0) = y_prime(i, 0) / U(i, i);
        for(unsigned j(i); j > 0;){
          j--;
          y_prime(j, 0) -= U(j, i) * x(i, 0);
        }
      }

      return x;
    }

    /**
     * UD分解をします。
     * (0, 0)〜(n-1,n-1):  U行列
     * (0, n)〜(n-1,2n-1): D行列
     *
     * @param do_check 対称行列チェックを行うか(デフォルトtrue)
     * @return (self_t) UD分解
     * @throw MatrixException 対称行列ではなくUD分解を計算することができない場合
     */
    self_t decomposeUD(bool do_check = true) const throw(MatrixException){
      if(do_check && !isSymmetric()){throw MatrixException("Operation void");}
      self_t P(copy());
      self_t UD(rows(), columns() * 2);
#define U(i, j) UD(i, j)
#define D(i, j) UD(i, j + columns())
      for(int i = rows() - 1; i >= 0; i--){
        D(i, i) = P(i, i);
        U(i, i) = elm_t(1);
        for(int j = 0; j < i; j++){
          U(j, i) = P(j, i) / D(i, i);
          for(int k = 0; k <= j; k++){
            P(k, j) -= U(k, i) * D(i, i) * U(j, i);
          }
        }
      }
#undef U
#undef D
      return UD;
    }

    /**
     * 2次小行列の固有値を求めます。
     *
     * @param row 2次小行列の左上項の行インデックス
     * @param column 2次小行列の左上項の列インデックス
     * @param upper 結果(固有値1)
     * @param lower 結果(固有値2)
     */
    void eigen22(const int &row, const int &column, Complex<elm_t> &upper, Complex<elm_t> &lower) const throw(MatrixException){
      elm_t a((*const_cast<Matrix *>(this))(row, column)),
        b((*const_cast<Matrix *>(this))(row, column + 1)),
        c((*const_cast<Matrix *>(this))(row + 1, column)),
        d((*const_cast<Matrix *>(this))(row + 1, column + 1));
      elm_t root(pow((a - d), 2) + b * c * 4);
      if(root >= elm_t(0)){
        root = ::sqrt(root);
        upper = Complex<elm_t>((a + d + root) / 2);
        lower = Complex<elm_t>((a + d - root) / 2);
      }else{
        root = ::sqrt(root * -1);
        upper = Complex<elm_t>((a + d) / 2, root / 2);
        lower = Complex<elm_t>((a + d) / 2, root / 2 * -1);
      }
    }

  protected:
    /**
     * 行列の平方根を求めます。
     * 返却値はMatrix型です。
     *
     * 行列Aが
     * @f[
     *    A = V D V^{-1}
     * @f]
     * と固有値(D)と固有ベクトル(V)に分解できれば
     * @f[
     *    A^{1/2} = V D^{1/2} V^{-1}
     * @f]
     * である。
     *
     * @param eigen_mat 固有値、固有ベクトルが入った(n,n+1)の行列
     * @return (Matrix<Complex<elm_t> >) 平方根
     */
    Matrix<Complex<elm_t> > sqrt(const Matrix<Complex<elm_t> > &eigen_mat) const throw(MatrixException) {
      int n(eigen_mat.rows());
      Matrix<Complex<elm_t> > VsD(eigen_mat.partial(n, n, 0, 0));
      Matrix<Complex<elm_t> > nV(VsD.inverse());
      for(int i(0); i < n; i++){
        VsD.partial(n, 1, 0, i) *= ::sqrt(const_cast<Matrix<Complex<elm_t> > &>(eigen_mat)(i, n));
      }

      return VsD * nV;
    }

  public:
    /**
     * 行列の平方根を求めます。
     * 返却値はMatrix型です。
     *
     * @param threshold_abs 固有値、固有ベクトル求める際に収束判定に用いる絶対誤差
     * @param threshold_abs 固有値、固有ベクトル求める際に収束判定に用いる相対誤差
     * @return (Matrix<Complex<elm_t> >) 平方根
     */
    Matrix<Complex<elm_t> > sqrt(
        const elm_t &threshold_abs,
        const elm_t &threshold_rel) const throw(MatrixException){
      return sqrt(eigen(threshold_abs, threshold_rel));
    }

    /**
     * 行列の平方根を求めます。
     * 返却値はMatrix型です。
     *
     * @return (Matrix<Complex<elm_t> >) 平方根
     */
    Matrix<Complex<elm_t> > sqrt() const throw(MatrixException){
      return sqrt(eigen());
    }
};

#endif

#endif /* __MATRIX_GSL_H */