Incomplete implementatiof Eskow algorithm

src/eskow.cpp

@@ -0,0 +1,109 @@
+
+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))
+      {
+      
+      }
+  }
+
+};