namespace CholeskyEskow
{

  template<typename T, typename A>
  void minmax_diag(A&m m, int j, int N, T& minval, T& maxval, int& i_min, int& i_max)
  {
    minval = maxval = m(j,j);

    for (int k = j+1; k < N; k++)
      {
	maxval = std::max(maxval, m(k,k));
	minval = std::min(minval, m(k,k));
      }
  
    for (int k = j; k < N; k++)
      {
	if (m(k,k) == minval)
	  i_min = k;
	if (m(k,k) == maxval)
	  i_max = k;
      }
  }

  template<typename T, typename A>
  void swap_rows(A& m, int N, int i0, int i1)
  {
    for (int r = 0; r < N; r++)
      std::swap(m(r,i0), m(r,i1));
  }

  template<typename T, typename A>
  void swap_cols(A& m, int N, int i0, int i1)
  {
    for (int c = 0; c < N; c++)
      std::swap(m(i0,c), m(i1,c));
  }

  template<typename T>
  T square(T x)
  {
    return x*x;
  }

  template<typename T, typename A>
  T max_row(A& m, int i, int j, int N)
  {
    T v = m(i,i) - square(m(i,j))/m(j,j);
  }


  template<typename T, typename A>
  void cholesky_eskow(A& m, int N)
  {
    T tau_bar = std::pow(mach_epsilon<T>(), 2./3);
    T tau = std::pow(mach_epsilon<T>(), 1./3);
    T mu = 0.1;
    bool phaseone = true;
    T gamma = max_diag(m, 0, N);
    int j;

    for (j = 0; j < N && phaseone; j++)
      {
	T minval, maxval;
      
	minmax_diag(m, j, N, minval, maxval, i_min, i_max);
	if (maxval < tau_bar*gamma || minval < -mu*maxval)
	  {
	    phaseone = false;
	    break;
	  }
      
	if (i_max != j)
	  {
	    swap_cols(m, N, i_max, j);
	    swap_rows(m, N, i_max, j);
	  }
      
	if (max_row(m, j, N) < -mu*gamma)
	  {
	    phaseone = false;
	    break;
	  }
      
	T L_jj = std::sqrt(m(j,j));

	m(j,j) = L_jj;
	for (int i = j+1; i < N; i++)
	  {
	    m(i,j) /= L_jj;
	    for (int k = j+1; k < i; k++)
	      m(i,k) -= m(i,j)*m(k,j);
	  }
      }

    if (!phaseone && j == N-1)
      {
	T A_nn = m(N-1,N-1);
	T delta = -A_nn + std::max(tau*(-A_nn)/(1-tau), tau_bar*gamma);
      
	m(N-1,N-1) = std::sqrt(m(N-1,N-1) + delta);	  	  
      }
    if (!phaseone && j < (N-1))
      {
      
      }
  }

};