guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Fixed decomposition. Add sample matrix.

Browse source
This commit is contained in:
Guilhem Lavaux 2011-12-14 17:40:10 -05:00
parent cf279bbb19
commit 196e17c242
3 changed files with 89 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

94
src/eskow.hpp
View file

@ -7,24 +7,49 @@
/* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */
namespace CholeskyEskow
template<typename T, typename A>
class CholeskyEskow
{
private:
static const bool verbose_eskow = true;
T tau, tau_bar, mu;
void print_matrix(A& m, int N)
{
using std::cout;
using std::endl;
using std::setprecision;
if (verbose_eskow)
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
cout.width(6);
cout << setprecision(5) << m(i,j) << " ";
}
cout << endl;
}
cout << endl;
}
}
template<typename T, typename A>
T max_diag(A& m, int j, int N)
{
T maxval = m(j,j);
T maxval = std::abs(m(j,j));
for (int k = j+1; k < N; k++)
{
maxval = std::max(maxval, m(k,k));
maxval = std::max(maxval, std::abs(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)
{
i_min = i_max = j;
minval = maxval = m(j,j);
for (int k = j+1; k < N; k++)
@ -35,34 +60,30 @@ namespace CholeskyEskow
for (int k = j; k < N; k++)
{
if (m(k,k) == minval)
if (m(k,k) == minval && i_min < 0)
i_min = k;
if (m(k,k) == maxval)
if (m(k,k) == maxval && i_max < 0)
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 min_row(A& m, int j, int N)
{
T a = 1/m(j,j);
@ -70,13 +91,12 @@ namespace CholeskyEskow
for (int i = j+2; i < N; i++)
{
v = std::max(v, m(i, i) - square(m(i,j))*a);
v = std::min(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];
@ -93,37 +113,45 @@ namespace CholeskyEskow
return k;
}
template<typename T, typename A>
void cholesky_eskow(A& m, int N, T& norm_E)
public:
CholeskyEskow()
{
tau = std::pow(mach_epsilon<T>(), 1./3);
tau_bar = std::pow(mach_epsilon<T>(), 2./3);
mu=0.1;
}
void cholesky(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<T,A>(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<T,A>(m, j, N, minval, maxval, i_min, i_max);
print_matrix(m, N);
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)
if (i_max != j)
{
swap_cols<T,A>(m, N, i_max, j);
swap_rows<T,A>(m, N, i_max, j);
std::cout << "Have to swap i=" << i_max << " and j=" << j << std::endl;
swap_cols(m, N, i_max, j);
swap_rows(m, N, i_max, j);
}
if (min_row<T,A>(m, j, N) < -mu*gamma)
if (min_row(m, j, N) < -mu*gamma)
{
phaseone = false;
break;
@ -135,7 +163,7 @@ namespace CholeskyEskow
for (int i = j+1; i < N; i++)
{
m(i,j) /= L_jj;
for (int k = j+1; k < i; k++)
for (int k = j+1; k <= i; k++)
m(i,k) -= m(i,j)*m(k,j);
}
}
@ -151,9 +179,10 @@ namespace CholeskyEskow
if (!phaseone && j < (N-1))
{
std::cout << "Phase two ! (j=" << j << ")" << std::endl;
int k = j-1;
std::vector<T> g(N);
@ -173,12 +202,15 @@ namespace CholeskyEskow
int i = g_max(g, j, N);
T norm_j;
print_matrix(m, N);
if (i != j)
{
swap_cols<T,A>(m, N, i, j);
swap_rows<T,A>(m, N, 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));
@ -208,7 +240,7 @@ namespace CholeskyEskow
for (int i = j+1; i < N; i++)
{
m(i,j) /= L_jj;
for (int k = j+1; k < i; k++)
for (int k = j+1; k <= i; k++)
m(i,k) -= m(i,j)*m(k,j);
}
}