summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt2
-rw-r--r--src/analyze_correlations.cpp99
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;
}