Got a compiling version of the algorithm

sample/CMakeLists.txt

@@ -26,3 +26,7 @@ if (HDF5_FOUND)
 	  add_executable(testReadFlash testReadFlash.cpp)
 	  target_link_libraries(testReadFlash ${tolink})
 endif (HDF5_FOUND)
+
+
+add_executable(testEskow testEskow.cpp)
+target_link_libraries(testEskow ${tolink})

sample/testEskow.cpp

@@ -0,0 +1,50 @@
+#include <cstring>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include "eskow.hpp"
+
+using namespace std;
+
+double Hartmann_Matrix[6][6] = {
+  { 14.8253,   -6.4243,   7.8746,  -1.2498,  10.2733,  10.2733 },
+  { -6.4243,   15.1024,  -1.1155,  -0.2761,  -8.2117,  -8.2117 },
+  {  7.8746,   -1.1155,  51.8519, -23.3482,  12.5902,  12.5902 },
+  { -1.2498,   -0.2761, -23.3482,  22.7962,  -9.8958,  -9.8958 },
+  { 10.2733,   -8.2117,  12.5902,  -9.8958,  21.0656,  21.0656 },
+  { 10.2733,   -8.2117,  12.5902,  -9.8958,  21.0656,  21.0656 }
+};
+
+struct MatrixOp
+{
+  vector<double> M;
+  int N;
+
+  double& operator()(int i, int j)
+  {
+    return M[i*N + j];
+  }
+};
+
+int main()
+{
+  MatrixOp M;
+  double norm_E;
+
+  M.N = 6;
+  M.M.resize(M.N*M.N);
+
+  memcpy(&M.M[0], &Hartmann_Matrix[0][0], sizeof(double)*36);
+
+  CholeskyEskow::cholesky_eskow(M, M.N, norm_E);
+
+  for (int i = 0; i < M.N; i++)
+    {
+      for (int j = 0; j < M.N; j++)
+	{
+	  cout << setprecision(25) << M(i,j) << " ";
+	}
+      cout << endl;
+    }
+  return 0;
+}

src/eskow.hpp

@@ -1,9 +1,29 @@
+#ifndef __ESKOW_CHOLESKY_HPP
+#define __ESKOW_CHOLESKY_HPP
+
+#include <cmath>
+#include <vector>
+#include "mach.hpp"
+
+/* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */
 
 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)
+  T max_diag(A& m, int j, int N)
+  {
+    T maxval = m(j,j);
+
+    for (int k = j+1; k < N; k++)
+      {
+	maxval = std::max(maxval, m(k,k));
+      }
+    return maxval;
+  }
+
+  template<typename T, typename A>
+  void minmax_diag(A& m, int j, int N, T& minval, T& maxval, int& i_min, int& i_max)
   {
     minval = maxval = m(j,j);
 
@@ -43,27 +63,54 @@ namespace CholeskyEskow
   }
 
   template<typename T, typename A>
-  T max_row(A& m, int i, int j, int N)
+  T min_row(A& m, int j, int N)
   {
-    T v = m(i,i) - square(m(i,j))/m(j,j);
+    T a = 1/m(j,j);
+    T v = m(j+1,j+1) - square(m(j+1,j))*a;
+    
+    for (int i = j+2; i < N; i++)
+      {
+	v = std::max(v, m(i, i) - square(m(i,j))*a);
       }
 
+    return v;
+  }
+
+  template<typename T>
+  int g_max(const std::vector<T>& g, int j, int N)
+  {
+    T a = g[j];
+    int k = j;
+
+    for (int i = j+1; i < N; i++)
+      {
+	if (a < g[i])
+	  {
+	    a = g[i];
+	    k = i;
+	  }
+      }
+    return k;
+  }
 
   template<typename T, typename A>
-  void cholesky_eskow(A& m, int N)
+  void cholesky_eskow(A& m, int N, T& norm_E)
   {
     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);
+    T gamma = max_diag<T,A>(m, 0, N);
     int j;
 
+    norm_E = 0;
+
     for (j = 0; j < N && phaseone; j++)
       {
 	T minval, maxval;
+	int i_min, i_max;
       
-	minmax_diag(m, j, N, minval, maxval, i_min, i_max);
+	minmax_diag<T,A>(m, j, N, minval, maxval, i_min, i_max);
 	if (maxval < tau_bar*gamma || minval < -mu*maxval)
 	  {
 	    phaseone = false;
@@ -72,11 +119,11 @@ namespace CholeskyEskow
       
 	if (i_max != j)
 	  {
-	    swap_cols(m, N, i_max, j);
-	    swap_rows(m, N, i_max, j);
+	    swap_cols<T,A>(m, N, i_max, j);
+	    swap_rows<T,A>(m, N, i_max, j);
 	  }
       
-	if (max_row(m, j, N) < -mu*gamma)
+	if (min_row<T,A>(m, j, N) < -mu*gamma)
 	  {
 	    phaseone = false;
 	    break;
@@ -93,6 +140,7 @@ namespace CholeskyEskow
 	  }
       }
 
+
     if (!phaseone && j == N-1)
       {
 	T A_nn = m(N-1,N-1);
@@ -100,10 +148,92 @@ namespace CholeskyEskow
       
 	m(N-1,N-1) = std::sqrt(m(N-1,N-1) + delta);	  	  
       }
+
+    
+
+
     if (!phaseone && j < (N-1))
       {
+	int k = j-1;
+	std::vector<T> g(N);
 
+	for (int i = k+1; i < N; i++)
+	  {
+	    g[i] = m(i,i);
+	    for (int j = k+1; j < i; j++)
+	      g[i] -= std::abs(m(i,j));
+	    for (int j = i+1; j < N; j++)
+	      g[i] -= std::abs(m(j,i));
+	  }
+
+	T delta, delta_prev = 0;
+	
+	for (int j = k+1; j < N-2; j++)
+	  {
+	    int i = g_max(g, j, N);
+	    T norm_j;
+
+	    if (i != j)
+	      {
+		swap_cols<T,A>(m, N, i, j);
+		swap_rows<T,A>(m, N, i, j);		
+	      }
+
+	    for (int i = j+1; j < N; j++)
+	      {
+		norm_j += std::abs(m(i,j));
+	      }
+	    
+	    delta = std::max(delta_prev, std::max((T)0, -m(j,j) + std::max(norm_j,tau_bar*gamma)));
+	    if (delta > 0)
+	      {
+		m(j,j) += delta;
+		delta_prev = delta;
+	      }
+	    
+	    if (m(j,j) != norm_j)
+	      {
+		T temp = 1 - norm_j/m(j,j);
+		
+		for (int i = j+1; j < N; j++)
+		  {
+		    g[i] += std::abs(m(i,j))*temp;
+		  }
+	      }
+
+	    // Now we do the classic cholesky iteration
+	    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);
+	      }
+	  }
+
+	// The final 2x2 submatrix is special
+	T A00 = m(N-2, N-2), A01 = m(N-2, N-1), A11 = m(N-1,N-1);
+	T sq_DELTA = std::sqrt(square(A00-A11) + square(A01));
+	T lambda_hi = 0.5*((A00+A11) + sq_DELTA);
+	T lambda_lo = 0.5*((A00+A11) - sq_DELTA);
+       
+	delta = std::max(std::max((T)0, -lambda_lo + std::max(tau*sq_DELTA/(1-tau), tau_bar*gamma)),delta_prev);
+	if (delta > 0)
+	  {
+	    m(N-1,N-1) += delta;
+	    m(N,N) += delta;
+	    delta_prev = delta;
+	  }
+	m(N-2,N-2) = A00 = std::sqrt(A00);
+	m(N-1,N-2) = (A01 /= A00);
+	m(N-1,N-1) = std::sqrt(A11-A01*A01);
+	norm_E = delta_prev;
       }
   }
 
 };
+
+
+#endif