/*****************************************************************

This software is part of Jean-Sébastien Caux's Abacus toolsuite.

Copyright © Jean-Sébastien Caux.

*****************************************************************/


export module matrix;


import std;


import conveniences;


// version using flattened vector

////////////////////
// ↓ Class Matrix //
////////////////////

export template <class T>
class Matrix {

private:
  std::size_t dim_;
  std::vector<T> element_;

public:
  Matrix (std::size_t dim);

  inline std::size_t size() const { return dim_; };

  inline T operator() (IndexU i, IndexU j) const { return element_[i*dim_ + j]; };
  inline T& operator() (IndexU i, IndexU j) { return element_[i*dim_ + j]; };

  void setZero();

  Matrix<T>& operator*= (const T& a);
  Matrix<T>& operator/= (const T& a);

  void ludcmp (std::vector<IndexU>& indx, T& d);
  void lubksb (std::vector<IndexU>& indx, std::vector<T>& b);
  std::complex<long double> lndet_LU_destroy ();
};

template <class T>
Matrix<T>::Matrix (std::size_t dim)
  : dim_(dim)
  , element_(std::vector<T>(dim_*dim_))
{}

template <class T>
void Matrix<T>::setZero ()
{
  std::fill(element_.begin(), element_.end(), T(0));
}


template <class T>
Matrix<T>& Matrix<T>::operator*= (const T& a)
{
  std::transform(element_.begin(), element_.end(), element_.begin(), [a](T el) { return el * a; });
  return *this;
}

template <class T>
Matrix<T>& Matrix<T>::operator/= (const T& a)
{
  T oneovera = T(1)/a;
  std::transform(element_.begin(), element_.end(), element_.begin(), [oneovera](T el) { return el * oneovera; });
  return *this;
}


template <class T>
void Matrix<T>::ludcmp (std::vector<IndexU>& indx, T& d)
{
  IndexU i, j, k;
  IndexU imax { 0 };
  IndexU idim_, jdim_, imaxdim_;
  T big, dum, sum, temp;

  IndexU n { size() };

  std::vector<T> vv(n);
  d = T(1);
  for (i = 0; i < n; i++) {
    big = T(0);
    idim_ = i*dim_;
    for (j = 0; j < n; j++) {
      if ((std::fabs(temp = element_[idim_ + j])) > std::fabs(big)) big = temp;
    }
    if (big == T(0)) throw;
    vv[i] = T(1)/big;
  }

  for (j = 0; j < n; j++) {
    for (i = 0; i < j; i++) {
      idim_ = i*dim_;
      sum = element_[idim_ + j];
      for (k = 0; k < i; k++) sum -= element_[idim_ + k] * element_[k*dim_ + j];
      element_[idim_ + j] = sum;
    }
    big = T(0);
    for (i = j; i < n; i++) {
      idim_ = i*dim_;
      sum = element_[idim_ + j];
      for (k = 0; k < j; k++) sum -= element_[idim_ + k] * element_[k*dim_ + j];
      element_[idim_ + j] = sum;
      if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) {
	big = dum;
	imax = i;
      }
    }
    jdim_ = j*dim_;
    if (j != imax) {
      imaxdim_ = imax*dim_;
      for (k = 0; k < n; k++) {
	dum = element_[imaxdim_ + k];
	element_[imaxdim_ + k] = element_[jdim_ + k];
	element_[jdim_ + k] = dum;
      }
      d = -d;
      vv[imax] = vv[j];
    }
    indx[j] = imax;
    if (j !=n-1) {
      dum = T(1)/(element_[jdim_ + j]);
      for (i = j + 1; i < n; i++) element_[i*dim_ + j] *= dum;
    }
  }
}

template <class T>
void Matrix<T>::lubksb (std::vector<IndexU>& indx, std::vector<T>& b)
{
  // int i, ip, j;
  // int ii { 0 };
  // int idim_;
  // T sum;

  // int n { int(size()) };
  // for (i = 0; i < n; i++) {
  //   ip = indx[i];
  //   sum = b[ip];
  //   b[ip] = b[i];
  //   idim_ = i*dim_;
  //   if (ii != 0)
  //     for (j = ii-1; j < i; j++) sum -= element_[idim_ + j] * b[j];
  //   else if (sum != T(0))
  //     ii = i + 1;
  //   b[i] = sum;
  // }
  // for (i = n - 1; i >= 0; i--) {
  //   sum = b[i];
  //   idim_ = i*dim_;
  //   for (j = i + 1; j < n; j++) sum -= element_[idim_ + j] * b[j];
  //   b[i] = sum/element_[idim_ + i];
  // }

  std::size_t i, ip, j;
  std::size_t ii { 0 };
  std::size_t idim_;
  T sum;
  std::size_t n { size() };
  for (i = 0; i < n; i++) {
    ip = indx[i];
    sum = b[ip];
    b[ip] = b[i];
    idim_ = i*dim_;
    if (ii != 0)
      for (j = ii-1; j < i; j++) sum -= element_[idim_ + j] * b[j];
    else if (sum != T(0))
      ii = i + 1;
    b[i] = sum;
  }
  for (i = n; i--;) {
    sum = b[i];
    idim_ = i*dim_;
    for (j = i + 1; j < n; j++) sum -= element_[idim_ + j] * b[j];
    b[i] = sum/element_[idim_ + i];
  }
}

template <class T>
std::complex<long double> Matrix<T>::lndet_LU_destroy ()
{
  std::vector<IndexU> indx(size());
  T d;
  std::complex<long double> lndet { 0 };

  (*this).ludcmp(indx, d);

  lndet = log(std::complex<long double>(d));

  for (IndexU j { 0 }; j < size(); j++) {
    lndet += log(std::complex<long double>(element_[j*dim_ + j]));
  }

  return lndet;
}


export template <class T>
std::ostream& operator<< (std::ostream& s, const Matrix<T>& M)
{
  for (IndexU i { 0 }; i < M.size(); ++i) {
    for (IndexU j { 0 }; j < M.size(); ++j) s << M(i,j) << "\t";
    s << std::endl;
  }
  return s;
}

////////////////////
// ↑ Class Matrix //
////////////////////



// Additional utilities


inline Real SIGN (const Real& a, const Real& b)
{
  return b >= Real(0) ? (a >= Real(0) ? a : -a) : (a >= Real(0) ? -a : a);
}

Real pythag(const Real a, const Real b)
{
  Real absa, absb;
  absa = std::abs(a);
  absb = std::abs(b);
  if (absa > absb) return absa * std::sqrt(Real(1) + (absb/absa)*(absb/absa));
  else return absb * std::sqrt(Real(1) + (absa/absb)*(absa/absb));
}


// Singular value decomposition
// Numerical Recipes section 2.6
// Restricted to square matrix due to use of Matrix class
export void svdcmp(Matrix<Real>& a, std::vector<Real>& w, Matrix<Real>& v)
{
  bool flag;
  int i, its, j, jj, k, l, nm;
  Real anorm, c, f, g, h, s, scale, x, y, z;

  int m { int(a.size()) };
  int n { int(a.size()) };

  std::vector<Real> rv1(n);
  g = Real(0); scale = Real(0); anorm = Real(0);

  for (i = 0; i < n; i++) {
    l = i+2;
    rv1[i] = scale*g;
    g = Real(0); s = Real(0); scale = Real(0);
    if (i < m) {
      for (k = i; k < m; k++) scale += std::abs(a(k,i));
      if (scale != Real(0)) {
	for (k = i; k < m; k++) {
	  a(k,i) /= scale;
	  s += a(k,i) * a(k,i);
	}
	f = a(i,i);
	g = -SIGN(std::sqrt(s), f);
	h = f*g-s;
	a(i,i) = f-g;
	for (j = l-1; j < n; j++) {
	  for (s=Real(0), k=i; k < m; k++) s += a(k,i)*a(k,j);
	  f=s/h;
	  for (k=i; k < m; k++) a(k,j) += f*a(k,i);
	}
	for (k=i; k < m; k++) a(k,i) *= scale;
      }
    }

    w[i] = scale*g;
    g = Real(0); s = Real(0); scale = Real(0);

    if (i+1 <= m && i != n) {
      for (k = l-1; k < n; k++) scale += std::abs(a(i,k));
      if (scale != Real(0)) {
	for (k=l-1; k < n; k++) {
	  a(i,k) /= scale;
	  s += a(i,k)*a(i,k);
	}
	f = a(i, l-1);
	g = -SIGN(std::sqrt(s), f);
	h = f*g - s;
	a(i, l-1) = f-g;
	for (k=l-1; k < n; k++) rv1[k] = a(i,k)/h;
	for (j = l-1; j < m; j++) {
	  for (s=Real(0), k=l-1; k < n; k++) s += a(j,k)*a(i,k);
	  for (k=l-1; k < n; k++) a(j,k) += s*rv1[k];
	}
	for (k=l-1; k < n; k++) a(i,k) *= scale;
      }
    }
    anorm = std::max(anorm, std::abs(w[i]) + std::abs(rv1[i]));
  }

  for (i = n-1; i >= 0; i--) { // accumulation of right-hand transforms
    if (i < n-1) {
      if (g != Real(0)) {
	for (j = l; j < n; j++) v(j,i) = (a(i,j)/a(i,l))/g;
	for (j=l; j < n; j++) {
	  for (s=Real(0), k=l; k < n; k++) s += a(i,k)*v(k,j);
	  for (k=l; k < n; k++) v(k,j) += s*v(k,i);
	}
      }
      for (j=l; j < n; j++) {
	v(i,j) = Real(0);
	v(j,i) = Real(0);
      }
    }
    v(i,i) = Real(1);
    g = rv1[i];
    l = i;
  }

  for (i = std::min(m,n) -1; i >= 0; i--) { // accumulation of left-hand transforms
    l = i+1;
    g = w[i];
    for (j=l; j < n; j++) a(i,j) = Real(0);
    if (g != Real(0)) {
      g = Real(1)/g;
      for (j=l; j < n; j++) {
	for (s=Real(0), k=l; k < m; k++) s += a(k,i) * a(k,j);
	f = (s/a(i,i))*g;
	for (k=i; k < m; k++) a(k,j) += f*a(k,i);
      }
      for (j=i; j < m; j++) a(j,i) *= g;
    }
    else for (j=i; j < m; j++) a(j,i) = Real(0);
    ++a(i,i);
  }

  for (k=n-1; k >= 0; k--) {
    for (its=0; its < 30; its++) {
      flag = true;
      for (l=k; l >= 0; l--) {
	nm = l-1;
	if (std::abs(rv1[l]) + anorm == anorm) {
	  flag = false;
	  break;
	}
	if (std::abs(w[nm]) + anorm == anorm) break;
      }
      if (flag) {
	c = Real(0);
	s = Real(1);
	for (i=l-1; i < k+1; i++) {
	  f = s*rv1[i];
	  rv1[i] = c*rv1[i];
	  if (std::abs(f) + anorm == anorm) break;
	  g = w[i];
	  h = pythag(f,g);
	  w[i] = h;
	  h = Real(1)/h;
	  c = g*h;
	  s = -f*h;
	  for (j=0; j < m; j++) {
	    y = a(j, nm);
	    z = a(j,i);
	    a(j, nm) = y*c + z*s;
	    a(j,i) = z*c - y*s;
	  }
	}
      }

      z = w[k];
      if (l == k) {
	if (z < Real(0)) {
	  w[k] = -z;
	  for (j=0; j < n; j++) v(j,k) = -v(j,k);
	}
	break;
      }

      if (its == 29) throw AbacusException("No convergence in 30 svdmp iterations");
      x = w[l];
      nm = k-1;
      y = w[nm];
      g = rv1[nm];
      h = rv1[k];
      f = ((y-z) * (y+z) + (g-h) * (g+h))/(2*h*y);
      g = pythag(f, Real(1));
      f = ((x-z)*(x+z) + h*((y/(f+SIGN(g,f)))-h))/x;
      c = Real(1); s = Real(1);  // Next QR transformation
      for (j=l; j <= nm; j++) {
	i = j+1;
	g = rv1[i];
	y = w[i];
	h = s*g;
	g = c*g;
	z = pythag(f,h);
	rv1[j] = z;
	c = f/z;
	s = h/z;
	f = x*c + g*s;
	g = g*c - x*s;
	h = y*s;
	y *= c;
	for (jj = 0; jj < n; jj++) {
	  x = v(jj, j);
	  z = v(jj, i);
	  v(jj, j) = x*c + z*s;
	  v(jj, i) = z*c - x*s;
	}
	z = pythag(f,h);
	w[j] = z;
	if (z) {
	  z = Real(1)/z;
	  c = f*z;
	  s = h*z;
	}
	f = c*g + s*y;
	x = c*y - s*g;
	for (jj=0; jj < m; jj++) {
	  y = a(jj, j);
	  z = a(jj, i);
	  a(jj, j) = y*c + z*s;
	  a(jj, i) = z*c - y*s;
	}
      }
      rv1[l] = Real(0);
      rv1[k] = f;
      w[k] = x;
    }
  }
}


// Householder reduction of a real symmetric matrix
// NR section 11.2
export void tred2 (Matrix<Real>& a, std::vector<Real>& d, std::vector<Real>& e)
{
  int l, k, j, i;
  Real scale, hh, h, g, f;

  int n { int(a.size()) };
  for (i = n-1; i > 0; i--) {
    l = i - 1;
    h = scale = Real(0);
    if (l > 0) {
      for (k = 0; k < l + 1; k++) scale += std::abs(a(i,k));
      if (scale == Real(0)) e[i] = a(i,l);
      else {
	for (k = 0; k < l + 1; k++) {
	  a(i,k) /= scale;
	  h += a(i,k) * a(i,k);
	}
	f = a(i,l);
	g = (f >= Real(0) ? -std::sqrt(h) : std::sqrt(h));
	e[i] = scale * g;
	h -= f * g;
	a(i,l) = f - g;
	f = Real(0);
	for (j = 0; j < l + 1; j++) {
	  a(j,i) = a(i,j)/h;
	  g = Real(0);
	  for (k = 0; k < j + 1; k++) g += a(j,k) * a(i,k);
	  for (k = j + 1; k < l + 1; k++) g += a(k,j) * a(i,k);
	  e[j] = g/h;
	  f += e[j] * a(i,j);
	}
	hh = f/(h +h);
	for (j = 0; j < l + 1; j++) {
	  f = a(i,j);
	  e[j] = g = e[j] - hh * f;
	  for (k = 0; k < j + 1; k++) a(j,k) -= (f * e[k] + g * a(i,k));
	}
      }
    }
    else e[i] = a(i,l);
    d[i] = h;
  }
  d[0] = Real(0);
  e[0] = Real(0);

  for (i = 0; i < n; i++) {
    l = i;
    if (d[i] != Real(0)) {
      for (j = 0; j < l; j++) {
	g = Real(0);
	for (k = 0; k < l; k++) g += a(i,k) * a(k,j);
	for (k = 0; k < l; k++) a(k,j) -= g * a(k,i);
      }
    }
    d[i] = a(i,i);
    a(i,i) = Real(1);
    for (j = 0; j < l; j++) a(j,i) = a(i,j) = Real(0);
  }
} // tred2



// QL algorithm with implicit shifts
// NR section 11.3
export void tqli (std::vector<Real>& d, std::vector<Real>& e, Matrix<Real>& z)
{
  int m, l, iter, i, k;
  Real s, r, p, g, f, dd, c, b;

  int n { int(d.size()) };
  for (i = 1; i < n; i++) e[i-1] = e[i];
  e[n-1] = Real(0);
  for (l = 0; l < n; l++) {
    iter = 0;
    do {
      for (m = l; m < n-1; m++) {
	dd = std::abs(d[m]) + std::abs(d[m+1]);
	if (std::abs(e[m]) + dd == dd) break;
      }
      if (m != l) {
	if (iter++ == 30) {
	  std::cout << "Too many iterations in tqli" << std::endl;
	  std::exit(1);
	}
	g = (d[l + 1] - d[l])/(2 * e[l]);
	r = pythag(g, Real(1));
	g = d[m] - d[l] + e[l]/(g + SIGN(r, g));
	s = c = Real(1);
	p = Real(0);
	for (i = m - 1; i >= l; i--) {
	  f = s * e[i];
	  b = c * e[i];
	  e[i + 1] = (r = pythag(f,g));
	  if (r == Real(0)) {
	    d[i + 1] -= p;
	    e[m] = Real(0);
	    break;
	  }
	  s = f/r;
	  c = g/r;
	  g = d[i + 1] - p;
	  r = (d[i] - g) * s + 2 * c * b;
	  d[i + 1] = g + (p = s * r);
	  g = c * r - b;
	  for (k = 0; k < n; k++) {
	    f = z(k, i + 1);
	    z(k, i + 1) = s * z(k, i) + c * f;
	    z(k, i) = c * z(k, i) - s * f;
	  }
	}
	if (r == Real(0) && i >= l) continue;
	d[l] -= p;
	e[l] = g;
	e[m] = Real(0);
      }
    } while (m != l);
  }
} // tqli




// Failed version using T**
//
// ////////////////////
// // ↓ Class Matrix //
// ////////////////////

// export template <class T>
// class Matrix {

// private:
//   std::size_t dim_;
//   T** element_;

// public:
//   Matrix (std::size_t dim);
//   // ~Matrix () = default;
//   ~Matrix ();

//   inline std::size_t size() { return dim_; };

//   // inline T* operator[] (const IndexU i);
//   // inline const T* operator[] (const IndexU i) const;

//   inline T operator() (IndexU i, IndexU j) const { return element_[i][j]; };
//   inline T& operator() (IndexU i, IndexU j) { return element_[i][j]; };

//   void setZero();

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

//   Matrix<T>& operator*= (const T& a);
//   Matrix<T>& operator/= (const T& a);

//   void ludcmp (std::vector<IndexU>& indx, T& d);
//   void lubksb (std::vector<IndexU>& indx, std::vector<T>& b);
//   std::complex<long double> lndet_LU_destroy ();
// };

// template <class T>
// Matrix<T>::Matrix (std::size_t dim)
//   : dim_(dim)
//   , element_(new T*[dim])
// {
//   // element_[0] = new T[dim_*dim_];
//   // for (IndexU i { 0 }; i + 1 < dim_; i++) element_[i+1] = element_[i] + dim_;
//   for (IndexU i { 0 }; i < dim_; ++i) element_[i] = new T[dim];
// }

// template <class T>
// Matrix<T>::~Matrix()
// {
//   // if (element_ != 0) {
//   //   delete[] (element_[0]);
//   //   delete[] (element_);
//   // }
//   for (IndexU i { 0 }; i < dim_; ++i) delete[] element_[i];
//   delete[] element_;
// }

// // template <class T>
// // inline T* Matrix<T>::operator[] (const IndexU i)
// // {
// //   return element_[i];
// // }

// // template <class T>
// // inline const T* Matrix<T>::operator[] (const IndexU i) const
// // {
// //   return element_[i];
// // }

// template <class T>
// void Matrix<T>::setZero ()
// {
//   for (IndexU i { 0 }; i < dim_; ++i)
//     for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] = T(0);
// }

// // template <class T>
// // Matrix<T>& Matrix<T>::operator= (const Matrix<T>& rhs)
// // {
// //   if (this != &rhs) {
// //     if (dim_ != rhs.dim_) {
// //       throw;
// //     }

// //     for (int i = 0; i < dim_; ++i)
// //       for (int j = 0; j < dim_; ++j) element_[i][j] = rhs.element_[i][j];
// //   }
// //   return *this;
// // }



// template <class T>
// Matrix<T>& Matrix<T>::operator*= (const T& a)
// {
//   for (IndexU i { 0 }; i < dim_; ++i)
//     for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] *= a;
//   return *this;
// }

// template <class T>
// Matrix<T>& Matrix<T>::operator/= (const T& a)
// {
//   T oneovera = T(1)/a;
//   for (IndexU i { 0 }; i < dim_; ++i)
//     for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] *= oneovera;
//   return *this;
// }


// template <class T>
// void Matrix<T>::ludcmp (std::vector<IndexU>& indx, T& d)
// {
//   IndexU i, j, k;
//   IndexU imax { 0 };
//   T big, dum, sum, temp;

//   IndexU n { size() };

//   std::vector<T> vv(n);
//   d = T(1);
//   for (i = 0; i < n; i++) {
//     big = T(0);
//     for (j = 0; j < n; j++) {
//       if ((std::fabs(temp = element_[i][j])) > std::fabs(big)) big = temp;
//     }
//     if (big == T(0)) throw;
//     vv[i] = T(1)/big;
//   }

//   for (j = 0; j < n; j++) {
//     for (i = 0; i < j; i++) {
//       sum = element_[i][j];
//       for (k = 0; k < i; k++) sum -= element_[i][k] * element_[k][j];
//       element_[i][j] = sum;
//     }
//     big = T(0);
//     for (i = j; i < n; i++) {
//       sum = element_[i][j];
//       for (k = 0; k < j; k++) sum -= element_[i][k] * element_[k][j];
//       element_[i][j] = sum;
//       if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) {
// 	big = dum;
// 	imax = i;
//       }
//     }
//     if (j != imax) {
//       for (k = 0; k < n; k++) {
// 	dum = element_[imax][k];
// 	element_[imax][k] = element_[j][k];
// 	element_[j][k] = dum;
//       }
//       d = -d;
//       vv[imax] = vv[j];
//     }
//     indx[j] = imax;
//     if (j !=n-1) {
//       dum = T(1)/(element_[j][j]);
//       for (i = j + 1; i < n; i++) element_[i][j] *= dum;
//     }
//   }
// }

// template <class T>
// void Matrix<T>::lubksb (std::vector<IndexU>& indx, std::vector<T>& b)
// {
//   int i, ip, j;
//   int ii { 0 };
//   T sum;

//   int n { int(size()) };
//   for (i = 0; i < n; i++) {
//     ip = indx[i];
//     sum = b[ip];
//     b[ip] = b[i];
//     if (ii != 0)
//       for (j = ii-1; j < i; j++) sum -= element_[i][j] * b[j];
//     else if (sum != T(0))
//       ii = i + 1;
//     b[i] = sum;
//   }
//   for (i = n - 1; i >= 0; i--) {
//     sum = b[i];
//     for (j = i + 1; j < n; j++) sum -= element_[i][j] * b[j];
//     b[i] = sum/element_[i][i];
//   }
// }

// template <class T>
// std::complex<long double> Matrix<T>::lndet_LU_destroy ()
// {
//   std::vector<IndexU> indx(size());
//   T d;
//   std::complex<long double> lndet { 0 };

//   (*this).ludcmp(indx, d);

//   lndet = log(std::complex<long double>(d));

//   for (IndexU j { 0 }; j < size(); j++) {
//     lndet += log(std::complex<long double>(element_[j][j]));
//   }

//   return lndet;
// }


// template <class T>
// std::ostream& operator<< (std::ostream& s, const Matrix<T>& M)
// {
//   for (IndexU i { 0 }; i < M.dim(); ++i) {
//     for (IndexU j { 0 }; j < M.dim(); ++j) s << M[i][j] << "\t";
//     s << std::endl;
//   }
//   return s;
// }

// ////////////////////
// // ↑ Class Matrix //
// ////////////////////



// // Version using vector of vector (not working)

// ////////////////////
// // ↓ Class Matrix //
// ////////////////////

// export template <class T>
// class Matrix {

// private:
//   std::size_t dim_;
//   std::vector<std::vector<T>> element_;

// public:
//   Matrix (std::size_t dim);
//   ~Matrix () = default;

//   inline std::size_t size() { return dim_; };
//   inline std::vector<std::vector<T>> rows () const { return element_; };
//   inline std::vector<T>& operator[] (const std::size_t i);
//   inline T operator() (IndexU i, IndexU j) const { return element_[i][j]; };
//   inline T& operator() (IndexU i, IndexU j) { return element_[i][j]; };

//   void setZero();

//   Matrix<T>& operator*= (const T& a);
//   Matrix<T>& operator/= (const T& a);

//   void ludcmp (std::vector<IndexU>& indx, T& d);
//   void lubksb (std::vector<IndexU>& indx, std::vector<T>& b);
//   std::complex<long double> lndet_LU_destroy ();
// };

// template <class T>
// inline std::vector<T>& Matrix<T>::operator[] (const std::size_t i)
// {
//   return element_[i];
// }

// template <class T>
// Matrix<T>::Matrix (std::size_t dim)
//   : dim_(dim)
//   , element_(std::vector<std::vector<T>>(dim_, std::vector<T>(dim_)))
// {
// }

// template <class T>
// Matrix<T>& Matrix<T>::operator*= (const T& a)
// {
//   for (auto &row : element_)
//     for (auto &col : row) col *= a;
//   return *this;
// }

// template <class T>
// Matrix<T>& Matrix<T>::operator/= (const T& a)
// {
//   for (auto &row : element_)
//     for (auto &col : row) col /= a;
//   return *this;
// }

// template <class T>
// void Matrix<T>::setZero ()
// {
//   for (auto &row : element_)
//     for (auto &col : row) col = T(0);
// }

// template <class T>
// void Matrix<T>::ludcmp (std::vector<IndexU>& indx, T& d)
// {
//   IndexU i, j, k;
//   IndexU imax { 0 };
//   T big, dum, sum, temp;

//   IndexU n { size() };

//   std::vector<T> vv(n);
//   d = T(1);
//   for (i = 0; i < n; i++) {
//     big = T(0);
//     for (j = 0; j < n; j++) {
//       if ((std::fabs(temp = element_[i][j])) > std::fabs(big)) big = temp;
//     }
//     if (big == T(0)) throw;
//     vv[i] = T(1)/big;
//   }

//   for (j = 0; j < n; j++) {
//     for (i = 0; i < j; i++) {
//       sum = element_[i][j];
//       for (k = 0; k < i; k++) sum -= element_[i][k] * element_[k][j];
//       element_[i][j] = sum;
//     }
//     big = T(0);
//     for (i = j; i < n; i++) {
//       sum = element_[i][j];
//       for (k = 0; k < j; k++) sum -= element_[i][k] * element_[k][j];
//       element_[i][j] = sum;
//       if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) {
// 	big = dum;
// 	imax = i;
//       }
//     }
//     if (j != imax) {
//       for (k = 0; k < n; k++) {
// 	dum = element_[imax][k];
// 	element_[imax][k] = element_[j][k];
// 	element_[j][k] = dum;
//       }
//       d = -d;
//       vv[imax] = vv[j];
//     }
//     indx[j] = imax;
//     if (j !=n-1) {
//       dum = T(1)/(element_[j][j]);
//       for (i = j + 1; i < n; i++) element_[i][j] *= dum;
//     }
//   }
// }

// template <class T>
// void Matrix<T>::lubksb (std::vector<IndexU>& indx, std::vector<T>& b)
// {
//   int i, ip, j;
//   int ii { 0 };
//   T sum;

//   int n { int(size()) };
//   for (i = 0; i < n; i++) {
//     ip = indx[i];
//     sum = b[ip];
//     b[ip] = b[i];
//     if (ii != 0)
//       for (j = ii-1; j < i; j++) sum -= element_[i][j] * b[j];
//     else if (sum != T(0))
//       ii = i + 1;
//     b[i] = sum;
//   }
//   for (i = n - 1; i >= 0; i--) {
//     sum = b[i];
//     for (j = i + 1; j < n; j++) sum -= element_[i][j] * b[j];
//     b[i] = sum/element_[i][i];
//   }
// }

// template <class T>
// std::complex<long double> Matrix<T>::lndet_LU_destroy ()
// {
//   std::vector<IndexU> indx(size());
//   T d;
//   std::complex<long double> lndet { 0 };

//   (*this).ludcmp(indx, d);

//   lndet = log(std::complex<long double>(d));

//   for (IndexU j { 0 }; j < size(); j++) {
//     lndet += log(std::complex<long double>(element_[j][j]));
//   }

//   return lndet;
// }


// template <class T>
// std::ostream& operator<< (std::ostream& s, const Matrix<T>& M)
// {
//   for (auto row : M.rows()) {
//     for (auto col : row) s << col << "\t";
//     s << std::endl;
//   }
//   return s;
// }

// ////////////////////
// // ↑ Class Matrix //
// ////////////////////