diff --git a/sample/Hartmann_Matrix.txt b/sample/Hartmann_Matrix.txt new file mode 100644 index 0000000..18df79d --- /dev/null +++ b/sample/Hartmann_Matrix.txt @@ -0,0 +1,7 @@ +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 diff --git a/sample/testEskow.cpp b/sample/testEskow.cpp index 05a92d5..2f173a6 100644 --- a/sample/testEskow.cpp +++ b/sample/testEskow.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -26,25 +27,39 @@ struct MatrixOp } }; -int main() +int main(int argc, char **argv) { MatrixOp M; double norm_E; + ifstream fi(argv[1]); + ofstream f("eskowed.txt"); + CholeskyEskow chol; - M.N = 6; + fi >> M.N; 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) << " "; + fi >> M(i,j); + if (j > i) + M(i,j) =0; } - cout << endl; + } + + chol.cholesky(M, M.N, norm_E); + + for (int i = 0; i < M.N; i++) + { + for (int j = 0; j < M.N; j++) + { + if (j > i) + f << "0 "; + else + f << setprecision(25) << M(i,j) << " "; + } + f << endl; } return 0; } diff --git a/src/eskow.hpp b/src/eskow.hpp index d18d56c..ac9b292 100644 --- a/src/eskow.hpp +++ b/src/eskow.hpp @@ -7,24 +7,49 @@ /* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */ -namespace CholeskyEskow +template +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 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 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 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 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 T square(T x) { return x*x; } - template 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 int g_max(const std::vector& g, int j, int N) { T a = g[j]; @@ -93,37 +113,45 @@ namespace CholeskyEskow return k; } - template - void cholesky_eskow(A& m, int N, T& norm_E) +public: + CholeskyEskow() + { + tau = std::pow(mach_epsilon(), 1./3); + tau_bar = std::pow(mach_epsilon(), 2./3); + mu=0.1; + } + + void cholesky(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); + + 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(m, N, i_max, j); - swap_rows(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(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 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(m, N, i, j); - swap_rows(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); } } diff --git a/src/h5_readFlash.cpp b/src/h5_readFlash.cpp index 0559a8f..0c78f58 100644 --- a/src/h5_readFlash.cpp +++ b/src/h5_readFlash.cpp @@ -179,13 +179,13 @@ void h5_read_flash3_particles (H5File* file, int* totalparticles, int* localnp, int* particle_offset, - float pos1[], - float pos2[], - float pos3[], - float vel1[], - float vel2[], - float vel3[], - int id[]) + float *pos1, + float *pos2, + float *pos3, + float *vel1, + float *vel2, + float *vel3, + int *id) { herr_t status; @@ -339,15 +339,27 @@ void h5_read_flash3_particles (H5File* file, /* convert buffer into particle struct */ + if (id) { for(p=0; p < (pcount); p++) { id[p+poffset] = (int) *(partBuffer+iptag-1+p*numProps); + } } + + if (pos1 && pos2 && pos3) { + for(p=0; p < (pcount); p++) { pos1[p+poffset] = (float) *(partBuffer+ipx-1+p*numProps); pos2[p+poffset] = (float) *(partBuffer+ipy-1+p*numProps); pos3[p+poffset] = (float) *(partBuffer+ipz-1+p*numProps); + } + } + + + if (vel1 && vel2 && vel3) { + for(p=0; p < (pcount); p++) { vel1[p+poffset] = (float) *(partBuffer+ipvx-1+p*numProps); vel2[p+poffset] = (float) *(partBuffer+ipvy-1+p*numProps); vel3[p+poffset] = (float) *(partBuffer+ipvz-1+p*numProps); } + } memspace.close(); //status = H5Sclose(memspace); diff --git a/src/h5_readFlash.hpp b/src/h5_readFlash.hpp index 273ca07..5187232 100644 --- a/src/h5_readFlash.hpp +++ b/src/h5_readFlash.hpp @@ -23,13 +23,13 @@ void h5_read_flash3_particles (H5File* file, int* totalparticles, int* localnp, int* particle_offset, - float pos1[], - float pos2[], - float pos3[], - float vel1[], - float vel2[], - float vel3[], - int id[]); + float *pos1, + float *pos2, + float *pos3, + float *vel1, + float *vel2, + float *vel3, + int *id); void h5_read_flash3_header_info(H5File* file, double* time, /* simulation time */ diff --git a/src/loadFlash.cpp b/src/loadFlash.cpp index cc58521..ada498d 100644 --- a/src/loadFlash.cpp +++ b/src/loadFlash.cpp @@ -51,29 +51,36 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags) data->Omega_Lambda = omegalambda; // particle data - for (int i = 0; i < 3; i++) { - data->Pos[i] = new float[data->NumPart]; - if (data->Pos[i] == 0) { - delete data; + if (loadflags& NEED_POSITION) { + for (int i = 0; i < 3; i++) { + data->Pos[i] = new float[data->NumPart]; + if (data->Pos[i] == 0) { + delete data; return 0; } - } + } } + + if (loadflags &NEED_VELOCITY) { for (int i = 0; i < 3; i++) { data->Vel[i] = new float[data->NumPart]; if (data->Vel[i] == 0) { delete data; return 0; } - } + } } + if (loadflags & NEED_GADGET_ID) { data->Id = new int[data->NumPart]; if (data->Id == 0) { delete data; return 0; } + } int offset = 0; + + if (loadflags &(NEED_GADGET_ID|NEED_POSITION|NEED_VELOCITY)) h5_read_flash3_particles(&file, &npart, &npart, &offset, data->Pos[0], data->Pos[1], data->Pos[2], data->Vel[0], data->Vel[1], data->Vel[2], @@ -81,8 +88,8 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags) for (int i = 0; i < 3; i++) { for (int n = 0; n < data->NumPart; n++) { - data->Pos[i][n] = data->Pos[i][n] * data->Hubble / kpc2cm; - data->Vel[i][n] = data->Vel[i][n] * data->Hubble / km2cm; + if (loadflags& NEED_POSITION) data->Pos[i][n] = data->Pos[i][n] * data->Hubble / kpc2cm; + if (loadflags&NEED_VELOCITY) data->Vel[i][n] = data->Vel[i][n] * data->Hubble / km2cm; } } diff --git a/src/mach.hpp b/src/mach.hpp index 60c98bf..abcb22a 100644 --- a/src/mach.hpp +++ b/src/mach.hpp @@ -1,5 +1,7 @@ -#ifndef __MACHINE_INFORMATION_HPP -#define __MACHINE_INFORMATION_HPP +#ifndef __COSMO_MACHINE_TEST_HPP +#define __COSMO_MACHINE_TEST_HPP + +#include template T mach_epsilon()