From 8d4419d2fd6be46ce0f1de487867231f02306849 Mon Sep 17 00:00:00 2001 From: Your Name Date: Tue, 13 Dec 2011 17:48:37 -0500 Subject: [PATCH] Got a compiling version of the algorithm --- sample/CMakeLists.txt | 4 ++ sample/testEskow.cpp | 50 ++++++++++++++ src/eskow.hpp | 150 +++++++++++++++++++++++++++++++++++++++--- 3 files changed, 194 insertions(+), 10 deletions(-) create mode 100644 sample/testEskow.cpp diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index c2c835a..82eabcd 100644 --- a/sample/CMakeLists.txt +++ b/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}) \ No newline at end of file diff --git a/sample/testEskow.cpp b/sample/testEskow.cpp new file mode 100644 index 0000000..05a92d5 --- /dev/null +++ b/sample/testEskow.cpp @@ -0,0 +1,50 @@ +#include +#include +#include +#include +#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 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; +} diff --git a/src/eskow.hpp b/src/eskow.hpp index a4a8884..d18d56c 100644 --- a/src/eskow.hpp +++ b/src/eskow.hpp @@ -1,9 +1,29 @@ +#ifndef __ESKOW_CHOLESKY_HPP +#define __ESKOW_CHOLESKY_HPP + +#include +#include +#include "mach.hpp" + +/* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */ namespace CholeskyEskow { template - 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 + 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 - 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 + int g_max(const std::vector& 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 - void cholesky_eskow(A& m, int N) + void cholesky_eskow(A& m, int N, T& norm_E) { T tau_bar = std::pow(mach_epsilon(), 2./3); T tau = std::pow(mach_epsilon(), 1./3); T mu = 0.1; bool phaseone = true; - T gamma = max_diag(m, 0, N); + T gamma = max_diag(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(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(m, N, i_max, j); + swap_rows(m, N, i_max, j); } - if (max_row(m, j, N) < -mu*gamma) + if (min_row(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 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(m, N, i, j); + swap_rows(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