diff options
| -rw-r--r-- | CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/analyze_correlations.cpp | 99 |
2 files changed, 73 insertions, 28 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ea6cce..59fe2eb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ endif() target_link_libraries(wolff_finite gsl m cblas fftw3 tcmalloc profiler) target_link_libraries(wolff_heisenberg gsl m cblas fftw3 tcmalloc profiler) target_link_libraries(wolff_planar gsl m cblas fftw3 tcmalloc profiler) -target_link_libraries(analyze_correlations gsl m cblas fftw3 tcmalloc profiler) +target_link_libraries(analyze_correlations gsl m cblas fftw3 profiler) install(TARGETS wolff_finite wolff_heisenberg wolff_planar analyze_correlations DESTINATION bin) diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp index 9f72498..0eca56d 100644 --- a/src/analyze_correlations.cpp +++ b/src/analyze_correlations.cpp @@ -4,32 +4,57 @@ #include <stdio.h> #include <stdlib.h> #include <getopt.h> +#include <fftw3.h> -template <class T> -double *compute_OO(count_t length, T *data, count_t N) { +double *compute_OO(count_t length, double *tmp_data, count_t N) { double *OO = (double *)calloc(2 + length, sizeof(double)); - OO[0] = (double)N; - for (count_t i = 0; i < N; i++) { + if (OO == NULL) { + printf("calloc failure\n"); + exit(EXIT_FAILURE); + } + + int better_N = pow(2, floor(log(N) / log(2))); + double *data = tmp_data + (N - better_N); + + OO[0] = (double)better_N; + + for (count_t i = 0; i < better_N; i++) { OO[1] += data[i]; - for (count_t j = 0; j < length; j++) { - if (i + j < N) { - OO[2 + j] += data[i] * data[i + j]; - } - } } - OO[1] /= N; + OO[1] /= better_N; + + fftw_plan plan = fftw_plan_r2r_1d(better_N, data, data, FFTW_R2HC, FFTW_ESTIMATE); + fftw_execute(plan); + + double *tmp = (double *)malloc(better_N * sizeof(double)); + + tmp[0] = data[0] * data[0]; + tmp[better_N / 2] = data[better_N/2] * data[better_N/2]; + + for (count_t i = 1; i < better_N / 2 - 1; i++) { + tmp[i] = data[i] * data[i] + data[better_N - i] * data[better_N - i]; + tmp[better_N - i] = 0; + } + + fftw_destroy_plan(plan); + + fftw_plan plan2 = fftw_plan_r2r_1d(better_N, tmp, tmp, FFTW_HC2R, FFTW_ESTIMATE); + fftw_execute(plan2); for (count_t j = 0; j < length; j++) { - OO[2 + j] /= N - j; + OO[2 + j] = tmp[j] / pow(better_N, 2); } + free(tmp); + fftw_destroy_plan(plan2); + return OO; } -float finite_energy(q_t nb, double *J, q_t q, double *H, v_t nv, v_t ne, uint32_t *bo, uint32_t *so) { - float energy = 0; +double finite_energy(q_t nb, double *J, q_t q, double *H, v_t nv, v_t ne, uint32_t *bo, uint32_t *so) { + double energy = 0; v_t tot = 0; for (q_t i = 0; i < nb - 1; i++) { @@ -120,11 +145,16 @@ int main (int argc, char *argv[]) { printf("Number of steps %lu is less than %" PRIcount ", nothing done.\n", N, drop); } else { - uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); + double *data_S = (double *)malloc(N * sizeof(double)); uint32_t *data_B = (uint32_t *)malloc((nb - 1) * N * sizeof(uint32_t)); uint32_t *data_M = (uint32_t *)malloc((q - 1) * N * sizeof(uint32_t)); - fread(data_S, N, sizeof(uint32_t), file_S); + uint32_t tmp; + for (count_t i = 0; i < N; i++) { + fread(&tmp, 1, sizeof(uint32_t), file_S); + data_S[i] = (double)tmp; + } + fread(data_B, N * (nb - 1), sizeof(uint32_t), file_B); fread(data_M, N * (q - 1), sizeof(uint32_t), file_M); @@ -133,11 +163,11 @@ int main (int argc, char *argv[]) { sprintf(filename_S, "wolff_%lu_S_OO.dat", id); FILE *file_S_new = fopen(filename_S, "wb"); - fwrite(OO_S, 2 + length, sizeof(double), file_S_new); + fwrite(OO_S, sizeof(double), 2 + length, file_S_new); fclose(file_S_new); free(data_S); - float *data_E = (float *)malloc(N * sizeof(float)); + double *data_E = (double *)malloc(N * sizeof(double)); for (count_t i = 0; i < N; i++) { data_E[i] = finite_energy(nb, J, q, H, nv, ne, data_B + (nb - 1) * i, data_M + (q - 1) * i); @@ -152,13 +182,18 @@ int main (int argc, char *argv[]) { sprintf(filename_B, "wolff_%lu_E_OO.dat", id); FILE *file_E = fopen(filename_B, "wb"); - fwrite(OO_E, 2 + length, sizeof(double), file_E); + if (file_E == NULL) { + printf("failed to open file\n"); + exit(EXIT_FAILURE); + } + fwrite(OO_E, sizeof(double), 2 + length, file_E); fclose(file_E); + printf("Correlation functions for %g steps written.\n", OO_S[0]); + free(OO_S); free(OO_E); - printf("Correlation functions for %lu steps written.\n", N - drop); } fclose(file_M); @@ -200,11 +235,17 @@ int main (int argc, char *argv[]) { printf("Number of steps %lu is less than %" PRIcount ", nothing done.\n", N, drop); } else { - uint32_t *data_S = (uint32_t *)malloc(N * sizeof(uint32_t)); - float *data_E = (float *)malloc(N * sizeof(float)); + double *data_S = (double *)malloc(N * sizeof(double)); + double *data_E = (double *)malloc(N * sizeof(double)); - fread(data_S, N, sizeof(uint32_t), file_S); - fread(data_E, N, sizeof(float), file_E); + uint32_t tmp; + float tmp2; + for (count_t i = 0; i < N; i++) { + fread(&tmp, 1, sizeof(uint32_t), file_S); + data_S[i] = (double)tmp; + fread(&tmp2, 1, sizeof(float), file_E); + data_E[i] = (double)tmp2; + } double *OO_S = compute_OO(length, data_S + drop, N - drop); double *OO_E = compute_OO(length, data_E + drop, N - drop); @@ -213,17 +254,19 @@ int main (int argc, char *argv[]) { sprintf(filename_E, "wolff_%lu_E_OO.dat", id); FILE *file_S_new = fopen(filename_S, "wb"); - FILE *file_E_new = fopen(filename_E, "wb"); - fwrite(OO_S, 1 + length, sizeof(double), file_S_new); - fwrite(OO_E, 1 + length, sizeof(double), file_E_new); + fwrite(OO_S, sizeof(double), 2 + length, file_S_new); fclose(file_S_new); + + FILE *file_E_new = fopen(filename_E, "wb"); + fwrite(OO_E, sizeof(double), 2 + length, file_E_new); fclose(file_E_new); + free(data_S); free(data_E); + printf("Correlation functions for %g steps written.\n", OO_S[0]); free(OO_S); free(OO_E); - printf("Correlation functions for %lu steps written.\n", N - drop); } free(filename_E); free(filename_S); @@ -233,8 +276,10 @@ int main (int argc, char *argv[]) { } } + free(model); fclose(metadata); + fftw_cleanup(); return 0; } |