From 8ab0002e44a38067543e6b62c89335973dd512f0 Mon Sep 17 00:00:00 2001 From: chortas Date: Thu, 13 Jan 2022 22:13:19 -0300 Subject: [PATCH 01/15] Remove logs --- fftma_module/gen/lib_src/Py_getvalues.c | 12 -- fftma_module/gen/lib_src/Py_kgeneration.c | 38 ----- fftma_module/gen/lib_src/build_real.c | 34 ----- fftma_module/gen/lib_src/cardsin.c | 35 ----- fftma_module/gen/lib_src/cgrid.c | 36 ----- fftma_module/gen/lib_src/clean_real.c | 36 ----- fftma_module/gen/lib_src/cov_value.c | 36 ----- fftma_module/gen/lib_src/covariance.c | 35 ----- fftma_module/gen/lib_src/cubic.c | 37 +---- fftma_module/gen/lib_src/exponential.c | 3 - fftma_module/gen/lib_src/fftma2.c | 35 ----- fftma_module/gen/lib_src/fourt.c | 36 +---- fftma_module/gen/lib_src/gammf.c | 35 ----- fftma_module/gen/lib_src/gasdev.c | 55 ------- fftma_module/gen/lib_src/gaussian.c | 2 - fftma_module/gen/lib_src/generate.c | 35 ----- fftma_module/gen/lib_src/length.c | 35 ----- fftma_module/gen/lib_src/log.c | 177 ---------------------- fftma_module/gen/lib_src/log.h | 49 ------ fftma_module/gen/lib_src/maxfactor.c | 35 ----- fftma_module/gen/lib_src/memory.c | 94 ------------ fftma_module/gen/lib_src/memory.h | 24 --- fftma_module/gen/lib_src/nor2log.c | 19 --- fftma_module/gen/lib_src/nugget.c | 3 - fftma_module/gen/lib_src/power.c | 2 - fftma_module/gen/lib_src/prebuild_gwn.c | 35 ----- fftma_module/gen/lib_src/ran2.c | 35 ----- fftma_module/gen/lib_src/spherical.c | 4 - fftma_module/gen/moduleFFTMA.c | 5 +- fftma_module/gen/setup.py | 2 - 30 files changed, 4 insertions(+), 1015 deletions(-) delete mode 100644 fftma_module/gen/lib_src/log.c delete mode 100644 fftma_module/gen/lib_src/log.h delete mode 100644 fftma_module/gen/lib_src/memory.c delete mode 100644 fftma_module/gen/lib_src/memory.h diff --git a/fftma_module/gen/lib_src/Py_getvalues.c b/fftma_module/gen/lib_src/Py_getvalues.c index f83dc5e..3c1c8db 100644 --- a/fftma_module/gen/lib_src/Py_getvalues.c +++ b/fftma_module/gen/lib_src/Py_getvalues.c @@ -1,6 +1,5 @@ #include "Py_py-api.h" #include "genlib.h" -#include "log.h" #include "geostat.h" #include "pressure.h" #include "simpio.h" @@ -44,8 +43,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario_mod* variogram, struct statistic_mod* stat, int* cores) { clock_t t = clock(); - - log_info("RESULT = in progress"); int i, varioNargs = 12, j = 0; PyObject* listvario; PyObject* vgr; @@ -59,7 +56,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario stat->variance = (double*)malloc(stat->nblock_var * sizeof(double)); if (stat->variance == NULL) { free(stat->mean); - log_error("RESULT = failed"); return 0; } @@ -80,7 +76,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario cores)) { free(stat->mean); free(stat->variance); - log_error("RESULT = failed"); return 0; } @@ -90,7 +85,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario if (variogram->var == NULL) { free(stat->mean); free(stat->variance); - log_error("RESULT = failed"); return 0; } variogram->vario = (int*)malloc(variogram->Nvario * sizeof(int)); @@ -98,7 +92,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario free(stat->mean); free(stat->variance); free(variogram->var); - log_error("RESULT = failed"); return 0; } variogram->alpha = (double*)malloc(variogram->Nvario * sizeof(double)); @@ -107,7 +100,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario free(stat->variance); free(variogram->var); free(variogram->vario); - log_error("RESULT = failed"); return 0; } variogram->scf = (double*)malloc(3 * variogram->Nvario * sizeof(double)); @@ -117,7 +109,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario free(variogram->var); free(variogram->vario); free(variogram->alpha); - log_error("RESULT = failed"); return 0; } variogram->ap = (double*)malloc(9 * variogram->Nvario * sizeof(double)); @@ -128,13 +119,11 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario free(variogram->vario); free(variogram->alpha); free(variogram->scf); - log_error("RESULT = failed"); return 0; } for (i = 0; i < variogram->Nvario; i++) { vgr = PyList_GetItem(listvario, i); if (PyTuple_Size(vgr) != 12) { - log_error("RESULT = failed"); return 0; } (variogram->var)[i] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); @@ -154,6 +143,5 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario t = clock() - t; double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - log_info("RESULT = success, ELAPSED = %f seconds", time_taken); return 1; } diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c index 3886db4..c79ad46 100644 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -4,8 +4,6 @@ #include "simpio.h" #include "toolsFFTMA.h" #include "toolsIO.h" -#include "log.h" -#include "memory.h" #include #include #include @@ -20,11 +18,6 @@ /* Y is the realization with mean and variance wanted */ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, struct vario_mod variogram, struct realization_mod* Z, struct realization_mod* Y, int n[3], int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, N; int typelog; @@ -34,15 +27,6 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, n[1] = 0; n[2] = 0; - log_info("RESULT = in progress, N = %d", N); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - generate(&seed, N, Z, cores); /*FFTMA*/ @@ -55,26 +39,4 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, typelog = stat.type + 2; nor2log(Y, typelog, Y); } - - printf("termino nor2log\n"); - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - printf("termino pykgeneration\n"); - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c index ece1c66..f563f85 100644 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -1,6 +1,4 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include #include @@ -22,21 +20,8 @@ /*ireal: vector defining the i-part */ void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, j, k, maille1; double temp; - log_info("RESULT = in progress, NTOT = %d, covar = %f, n[0] = %d, n[1] = %d, n[2] = %d", NTOT, *covar, n[0], n[1], n[2]); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } /*decomposition and multiplication in the spectral domain*/ for (k = 1; k <= n[2]; k++) { @@ -54,23 +39,4 @@ void build_real(int n[3], int NTOT, double* covar, double* realization, double* } } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, realization = %f, ireal = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", *realization, *ireal, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/cardsin.c b/fftma_module/gen/lib_src/cardsin.c index f335e1d..b6a1ba7 100644 --- a/fftma_module/gen/lib_src/cardsin.c +++ b/fftma_module/gen/lib_src/cardsin.c @@ -1,26 +1,10 @@ #include "genlib.h" -#include "log.h" -#include "memory.h" #include #include #include /*cardsin covariance function*/ double cardsin(double h, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, h = %f", h); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - float delta = 20.371; double z; @@ -31,24 +15,5 @@ double cardsin(double h, int cores) { z = 1.; } - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", z, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return z; } diff --git a/fftma_module/gen/lib_src/cgrid.c b/fftma_module/gen/lib_src/cgrid.c index ed8b981..252ae35 100644 --- a/fftma_module/gen/lib_src/cgrid.c +++ b/fftma_module/gen/lib_src/cgrid.c @@ -1,6 +1,4 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include @@ -13,23 +11,9 @@ /* X, Y and Z axes for the underlying grid */ /* i = [0 1 2] */ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3], int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, N; double D; - log_info("RESULT = in progress"); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - if (n == NULL || n[0] == 0 || n[1] == 0 || n[2] == 0) { for (i = 0; i < 3; i++) { switch (i) { @@ -50,27 +34,7 @@ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3], int cores } } else { if ((n[0] < grid.NX) || (n[1] < grid.NY) || (n[2] < grid.NZ)) { - log_error("RESULT = failed - Indicated dimensions are inappropriate in cgrid"); exit; } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, n[0] = %d, n[1] = %d, n[2] = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", n[0], n[1], n[2], time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/clean_real.c b/fftma_module/gen/lib_src/clean_real.c index e0da069..9917554 100644 --- a/fftma_module/gen/lib_src/clean_real.c +++ b/fftma_module/gen/lib_src/clean_real.c @@ -1,6 +1,4 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include #include @@ -10,11 +8,6 @@ #include void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, j, k, maille0, maille1; double NTOT; @@ -22,19 +15,9 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, /*is the output realization already allocated?*/ /*if not, memory allocation*/ - log_info("RESULT = in progress, NTOT = %f", NTOT); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - if (realout->vector == NULL || realout->n != realin->n) { realout->vector = (double*)malloc(realin->n * sizeof(double)); if (realout->vector == NULL) { - log_error("RESULT = failed - No memory available"); exit; } } @@ -52,23 +35,4 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, } } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/cov_value.c b/fftma_module/gen/lib_src/cov_value.c index 2acbda5..4e30f90 100644 --- a/fftma_module/gen/lib_src/cov_value.c +++ b/fftma_module/gen/lib_src/cov_value.c @@ -1,26 +1,10 @@ #include "genlib.h" #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include /*selection of model covariance*/ double cov_value(struct vario_mod variogram, double di, double dj, double dk, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, di = %f, dj = %f, dk = %f", di, dj, dk); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - double hx, hy, hz, h; double cov; int k; @@ -64,25 +48,5 @@ double cov_value(struct vario_mod variogram, double di, double dj, double dk, in break; } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, cov = %f, hx = %f, hy = %f, hz = %f, h = %f , ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", cov, hx, hy, hz, h, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return cov; } diff --git a/fftma_module/gen/lib_src/covariance.c b/fftma_module/gen/lib_src/covariance.c index 6afc796..1bc9064 100644 --- a/fftma_module/gen/lib_src/covariance.c +++ b/fftma_module/gen/lib_src/covariance.c @@ -1,28 +1,12 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include /*builds the sampled covariance function*/ /*dimensions are even*/ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, j, k, maille, n2[3], symmetric; double di, dj, dk; - log_info("RESULT = in progress"); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - for (i = 0; i < 3; i++) n2[i] = n[i] / 2; @@ -90,23 +74,4 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, } } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, di = %f, dj = %f, dk = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", di, dj, dk, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/cubic.c b/fftma_module/gen/lib_src/cubic.c index a844b53..ac36add 100644 --- a/fftma_module/gen/lib_src/cubic.c +++ b/fftma_module/gen/lib_src/cubic.c @@ -1,26 +1,10 @@ #include "genlib.h" -#include "log.h" -#include "memory.h" #include #include #include /*cubic covariance function*/ double cubic(double h, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, h = %f", h); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - double z; if (h >= 1.) { @@ -28,25 +12,6 @@ double cubic(double h, int cores) { } else { z = 1. - 7. * (double)(h * h) + (35. / 4.) * (double)(h * h * h) - 3.5 * (double)(h * h * h * h * h) + .75 * (double)(h * h * h * h * h * h * h); } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", z, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - + return z; } diff --git a/fftma_module/gen/lib_src/exponential.c b/fftma_module/gen/lib_src/exponential.c index ab06305..af1ef27 100644 --- a/fftma_module/gen/lib_src/exponential.c +++ b/fftma_module/gen/lib_src/exponential.c @@ -1,11 +1,8 @@ #include "genlib.h" -#include "log.h" #include #include /*exponential covariance function*/ double exponential(double h) { - log_info("RESULT = in progress, h = %f", h); - return (exp(-3. * (double)h)); } diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 6c24e72..d699498 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -1,6 +1,4 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include #include @@ -25,20 +23,6 @@ /*realout: structure defining a realization - */ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct realization_mod* realin, struct realization_mod* realout, int cores, long* seed) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress"); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz; int solver; double temp; @@ -99,9 +83,6 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r /* build realization in spectral domain */ build_real(n, NTOT, covar, realization, ireal, cores); - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - free(covar); /*backward fourier transform*/ @@ -113,20 +94,4 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r /*output realization*/ clean_real(realin, n, grid, realization, realout, cores); - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - log_info("RESULT = success, NTOT = %d, NMAX = %d, NDIM = %d, ntot = %d, nmax = %d, NXYZ = %d, nxyz = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", NTOT, NMAX, NDIM, ntot, nmax, NXYZ, nxyz, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index df3814b..fc4184c 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -1,5 +1,3 @@ -#include "log.h" -#include "memory.h" #include #include #include @@ -93,20 +91,6 @@ /* 10-06-2000, MLR */ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress"); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; @@ -598,23 +582,5 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp np1 = np2; nprev = n; } -L920: - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - log_info("RESULT = success, ELAPSED = %f, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); +L920: return; } diff --git a/fftma_module/gen/lib_src/gammf.c b/fftma_module/gen/lib_src/gammf.c index 8269609..90e5cea 100644 --- a/fftma_module/gen/lib_src/gammf.c +++ b/fftma_module/gen/lib_src/gammf.c @@ -1,50 +1,15 @@ #include "genlib.h" -#include "log.h" -#include "memory.h" #include #include #include /*gamma covariance function*/ double gammf(double h, double alpha, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, h = %f, alpha = %f", h, alpha); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - float delta; double z; delta = pow(20., 1. / alpha) - 1.; z = 1. / (double)(pow(1. + h * delta, alpha)); - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - log_info("RESULT = success, delta = %f, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", delta, z, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); return z; } diff --git a/fftma_module/gen/lib_src/gasdev.c b/fftma_module/gen/lib_src/gasdev.c index 4511bb6..7388087 100644 --- a/fftma_module/gen/lib_src/gasdev.c +++ b/fftma_module/gen/lib_src/gasdev.c @@ -1,6 +1,4 @@ #include "genlib.h" -#include "log.h" -#include "memory.h" #include #include @@ -10,20 +8,6 @@ double gasdev(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { /*returns a normally distributed deviate with 0 mean*/ /*and unit variance, using ran2(idum) as the source */ /*of uniform deviates */ - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, idum = %f, idum2 = %f, iy = %f", *idum, *idum2, *iy); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores); static int iset = 0; static double gset; @@ -40,48 +24,9 @@ double gasdev(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { gset = v1 * fac; iset = 1; - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, gset = %f, fac = %f, v1 = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", gset, fac, v1, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return (v2 * fac); } else { iset = 0; - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, gset = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", gset, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return gset; } } diff --git a/fftma_module/gen/lib_src/gaussian.c b/fftma_module/gen/lib_src/gaussian.c index d571f70..98bba6e 100644 --- a/fftma_module/gen/lib_src/gaussian.c +++ b/fftma_module/gen/lib_src/gaussian.c @@ -1,10 +1,8 @@ #include "genlib.h" -#include "log.h" #include #include /*gaussian covariance function*/ double gaussian(double h) { - log_info("RESULT = in progress, h = %f", h); return (exp(-3. * (double)(h * h))); } diff --git a/fftma_module/gen/lib_src/generate.c b/fftma_module/gen/lib_src/generate.c index 8d5ca42..84e4ed0 100644 --- a/fftma_module/gen/lib_src/generate.c +++ b/fftma_module/gen/lib_src/generate.c @@ -2,8 +2,6 @@ #include #include #include -#include "log.h" -#include "memory.h" /* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */ /*input: */ @@ -13,20 +11,6 @@ /* realization: structure defining the realization*/ void generate(long* seed, int n, struct realization_mod* realization, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - log_info("RESULT = in progress, n = %d", n); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - - clock_t t = clock(); - int i; long idum2 = 123456789, iy = 0; long* iv; @@ -55,24 +39,5 @@ void generate(long* seed, int n, struct realization_mod* realization, int cores) chunk_array_flush((*realization).vector_2);*/ - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - free(iv); - - log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/length.c b/fftma_module/gen/lib_src/length.c index 5b705fe..e47410a 100644 --- a/fftma_module/gen/lib_src/length.c +++ b/fftma_module/gen/lib_src/length.c @@ -1,24 +1,8 @@ -#include "log.h" -#include "memory.h" #include #include /* compute the length for one dimension*/ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, N = %d, i = %d, D = %f, Nvari = %d", N, i, D, Nvari); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - int maxfactor(int n, int cores); double temp1, temp2; int n, j, k, nmax; @@ -52,24 +36,5 @@ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores } } - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, n = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", n, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return n; } diff --git a/fftma_module/gen/lib_src/log.c b/fftma_module/gen/lib_src/log.c deleted file mode 100644 index 287ac1f..0000000 --- a/fftma_module/gen/lib_src/log.c +++ /dev/null @@ -1,177 +0,0 @@ -/* - * Copyright (c) 2020 rxi - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to - * deal in the Software without restriction, including without limitation the - * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or - * sell copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in - * all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS - * IN THE SOFTWARE. - */ - -#include -#include "log.h" -#include - -#define MAX_CALLBACKS 32 - -typedef struct { - log_LogFn fn; - void *udata; - int level; -} Callback; - -static struct { - void *udata; - log_LockFn lock; - int level; - bool quiet; - Callback callbacks[MAX_CALLBACKS]; -} L; - - -static const char *level_strings[] = { - "TRACE", "DEBUG", "INFO", "WARN", "ERROR", "FATAL" -}; - -#ifdef LOG_USE_COLOR -static const char *level_colors[] = { - "\x1b[94m", "\x1b[36m", "\x1b[32m", "\x1b[33m", "\x1b[31m", "\x1b[35m" -}; -#endif - - -static void stdout_callback(log_Event *ev) { - char buf[16]; - buf[strftime(buf, sizeof(buf), "%H:%M:%S", ev->time)] = '\0'; -#ifdef LOG_USE_COLOR - fprintf( - ev->udata, "%s %s%-5s\x1b[0m \x1b[0m%s:%d:\x1b[0m ", - buf, level_colors[ev->level], level_strings[ev->level], - ev->file, ev->line); -#else - fprintf( - ev->udata, "%s %-5s %s:%d: ", - buf, level_strings[ev->level], ev->file, ev->line); -#endif - vfprintf(ev->udata, ev->fmt, ev->ap); - fprintf(ev->udata, "\n"); - fflush(ev->udata); -} - - -static void file_callback(log_Event *ev) { - char buf[64]; - buf[strftime(buf, sizeof(buf), "%Y-%m-%d %H:%M:%S", ev->time)] = '\0'; - fprintf( - ev->udata, "%s %-5s %s:%d: ", - buf, level_strings[ev->level], ev->file, ev->line); - vfprintf(ev->udata, ev->fmt, ev->ap); - fprintf(ev->udata, "\n"); - fflush(ev->udata); -} - - -static void lock(void) { - if (L.lock) { L.lock(true, L.udata); } -} - - -static void unlock(void) { - if (L.lock) { L.lock(false, L.udata); } -} - - -const char* log_level_string(int level) { - return level_strings[level]; -} - - -void log_set_lock(log_LockFn fn, void *udata) { - L.lock = fn; - L.udata = udata; -} - - -void log_set_level(int level) { - L.level = level; -} - - -void log_set_quiet(bool enable) { - L.quiet = enable; -} - - -int log_add_callback(log_LogFn fn, void *udata, int level) { - for (int i = 0; i < MAX_CALLBACKS; i++) { - if (!L.callbacks[i].fn) { - L.callbacks[i] = (Callback) { fn, udata, level }; - return 0; - } - } - return -1; -} - - -int log_add_fp(FILE *fp, int level) { - return log_add_callback(file_callback, fp, level); -} - - -static void init_event(log_Event *ev, void *udata) { - if (!ev->time) { - time_t t = time(NULL); - ev->time = localtime(&t); - } - ev->udata = udata; -} - - -void log_log(int level, const char *file, int line, const char *fmt, ...) { - log_Event ev = { - .fmt = fmt, - .file = file, - .line = line, - .level = level, - }; - - char* env_var = getenv("ENV"); - if (env_var != NULL && strcmp("false", env_var) == 0) return; - - char* substr_mem = strstr(fmt, "MEM"); - char* substr_cpu = strstr(fmt, "CPU"); - if (env_var != NULL && strcmp("analysis", env_var) == 0 && substr_mem == NULL && substr_cpu == NULL) return; - - lock(); - - if (!L.quiet && level >= L.level) { - init_event(&ev, stderr); - va_start(ev.ap, fmt); - stdout_callback(&ev); - va_end(ev.ap); - } - - for (int i = 0; i < MAX_CALLBACKS && L.callbacks[i].fn; i++) { - Callback *cb = &L.callbacks[i]; - if (level >= cb->level) { - init_event(&ev, cb->udata); - va_start(ev.ap, fmt); - cb->fn(&ev); - va_end(ev.ap); - } - } - - unlock(); -} diff --git a/fftma_module/gen/lib_src/log.h b/fftma_module/gen/lib_src/log.h deleted file mode 100644 index b1fae24..0000000 --- a/fftma_module/gen/lib_src/log.h +++ /dev/null @@ -1,49 +0,0 @@ -/** - * Copyright (c) 2020 rxi - * - * This library is free software; you can redistribute it and/or modify it - * under the terms of the MIT license. See `log.c` for details. - */ - -#ifndef LOG_H -#define LOG_H - -#include -#include -#include -#include - -#define LOG_VERSION "0.1.0" - -typedef struct { - va_list ap; - const char *fmt; - const char *file; - struct tm *time; - void *udata; - int line; - int level; -} log_Event; - -typedef void (*log_LogFn)(log_Event *ev); -typedef void (*log_LockFn)(bool lock, void *udata); - -enum { LOG_TRACE, LOG_DEBUG, LOG_INFO, LOG_WARN, LOG_ERROR, LOG_FATAL }; - -#define log_trace(...) log_log(LOG_TRACE, __FILE__, __LINE__, __VA_ARGS__) -#define log_debug(...) log_log(LOG_DEBUG, __FILE__, __LINE__, __VA_ARGS__) -#define log_info(...) log_log(LOG_INFO, __FILE__, __LINE__, __VA_ARGS__) -#define log_warn(...) log_log(LOG_WARN, __FILE__, __LINE__, __VA_ARGS__) -#define log_error(...) log_log(LOG_ERROR, __FILE__, __LINE__, __VA_ARGS__) -#define log_fatal(...) log_log(LOG_FATAL, __FILE__, __LINE__, __VA_ARGS__) - -const char* log_level_string(int level); -void log_set_lock(log_LockFn fn, void *udata); -void log_set_level(int level); -void log_set_quiet(bool enable); -int log_add_callback(log_LogFn fn, void *udata, int level); -int log_add_fp(FILE *fp, int level); - -void log_log(int level, const char *file, int line, const char *fmt, ...); - -#endif diff --git a/fftma_module/gen/lib_src/maxfactor.c b/fftma_module/gen/lib_src/maxfactor.c index f3bc12b..95c285c 100644 --- a/fftma_module/gen/lib_src/maxfactor.c +++ b/fftma_module/gen/lib_src/maxfactor.c @@ -1,23 +1,7 @@ #include "genlib.h" -#include "log.h" -#include "memory.h" /*determines the greatest prime factor of an integer*/ int maxfactor(int n, int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress, n = %d", n); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - int test_fact(int* pnum, int fact, int* pmaxfac); int lnum, fact; int maxfac; @@ -47,24 +31,5 @@ int maxfactor(int n, int cores) { } } - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, maxfac = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", maxfac, time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); - return maxfac; } diff --git a/fftma_module/gen/lib_src/memory.c b/fftma_module/gen/lib_src/memory.c deleted file mode 100644 index 36c2746..0000000 --- a/fftma_module/gen/lib_src/memory.c +++ /dev/null @@ -1,94 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include "log.h" -#include "memory.h" - -static unsigned long long lastTotalUser, lastTotalUserLow, lastTotalSys, lastTotalIdle; -static clock_t lastCPU, lastSysCPU, lastUserCPU; -static int numProcessors; - -void getTotalVirtualMem(double* total_ram) { - const double megabyte = 1024 * 1024; - struct sysinfo si; - sysinfo(&si); - *total_ram = si.totalram / megabyte; -} - -void getVirtualMemUsed(double* used_ram) { - const double megabyte = 1024 * 1024; - struct sysinfo si; - sysinfo(&si); - *used_ram = (si.totalram - si.freeram) / megabyte; -} - -int parseLine(char* line) { - // This assumes that a digit will be found and the line ends in " Kb". - int i = strlen(line); - const char* p = line; - while (*p <'0' || *p > '9') p++; - line[i-3] = '\0'; - i = atoi(p); - return i; -} - -int getVirtualMemUsedByCurrentProcess() { - FILE* file = fopen("/proc/self/status", "r"); - int result = -1; - char line[128]; - - while (fgets(line, 128, file) != NULL){ - if (strncmp(line, "VmSize:", 7) == 0){ - result = parseLine(line); - break; - } - } - fclose(file); - return result / 1024; -} - -void skip_lines(FILE *fp, int numlines) { - int cnt = 0; - char ch; - while ((cnt < numlines) && ((ch = getc(fp)) != EOF)) { - if (ch == '\n') cnt++; - } -} - -void get_stats(struct cpustat *st, int cpunum) { - FILE *fp = fopen("/proc/stat", "r"); - int lskip = cpunum+1; - skip_lines(fp, lskip); - char cpun[255]; - fscanf(fp, "%s %d %d %d %d %d %d %d", cpun, &(st->t_user), &(st->t_nice), - &(st->t_system), &(st->t_idle), &(st->t_iowait), &(st->t_irq), - &(st->t_softirq)); - fclose(fp); - return; -} - -double calculate_load(struct cpustat *prev, struct cpustat *cur) { - int idle_prev = (prev->t_idle) + (prev->t_iowait); - int idle_cur = (cur->t_idle) + (cur->t_iowait); - - int nidle_prev = (prev->t_user) + (prev->t_nice) + (prev->t_system) + (prev->t_irq) + (prev->t_softirq); - int nidle_cur = (cur->t_user) + (cur->t_nice) + (cur->t_system) + (cur->t_irq) + (cur->t_softirq); - - int total_prev = idle_prev + nidle_prev; - int total_cur = idle_cur + nidle_cur; - - double totald = (double) total_cur - (double) total_prev; - - double idled = (double) idle_cur - (double) idle_prev; - - if (totald == 0 && idled == 0) return 0; - - double cpu_perc = (1000 * (totald - idled) / totald + 1) / 10; - - return cpu_perc; -} diff --git a/fftma_module/gen/lib_src/memory.h b/fftma_module/gen/lib_src/memory.h deleted file mode 100644 index a462325..0000000 --- a/fftma_module/gen/lib_src/memory.h +++ /dev/null @@ -1,24 +0,0 @@ -#include "sys/types.h" -#include "sys/sysinfo.h" -#include "stdlib.h" -#include "stdio.h" -#include "string.h" -#include "sys/times.h" -#include "sys/vtimes.h" - -void getTotalVirtualMem(); -void getVirtualMemUsed(); -int getVirtualMemUsedByCurrentProcess(); - -struct cpustat { - unsigned long t_user; - unsigned long t_nice; - unsigned long t_system; - unsigned long t_idle; - unsigned long t_iowait; - unsigned long t_irq; - unsigned long t_softirq; -}; - -void get_stats(struct cpustat *st, int cpunum); -double calculate_load(struct cpustat *prev, struct cpustat *cur); \ No newline at end of file diff --git a/fftma_module/gen/lib_src/nor2log.c b/fftma_module/gen/lib_src/nor2log.c index 2e96e68..14367f9 100644 --- a/fftma_module/gen/lib_src/nor2log.c +++ b/fftma_module/gen/lib_src/nor2log.c @@ -1,5 +1,4 @@ #include "geostat.h" -#include "log.h" #include #include #include @@ -15,12 +14,6 @@ /* lognormal numbers */ void nor2log(struct realization_mod* realin, int typelog, struct realization_mod* realout) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - - log_info("RESULT = in progress"); int i; double coeff; @@ -69,20 +62,8 @@ void nor2log(struct realization_mod* realin, int typelog, struct realization_mod (*realout).vector[i] = exp((*realin).vector[i] * coeff); break; default: - log_error("RESULT = failed - Unexpected case in nor2log"); return; break; } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/nugget.c b/fftma_module/gen/lib_src/nugget.c index 4aec020..292138e 100644 --- a/fftma_module/gen/lib_src/nugget.c +++ b/fftma_module/gen/lib_src/nugget.c @@ -1,12 +1,9 @@ #include "genlib.h" -#include "log.h" #include #include /*nugget covariance function*/ double nugget(double h) { - log_info("RESULT = in progress, h = %f", h); - if (h == 0) { return (1.); } else { diff --git a/fftma_module/gen/lib_src/power.c b/fftma_module/gen/lib_src/power.c index b0f3d40..5e6cdb8 100644 --- a/fftma_module/gen/lib_src/power.c +++ b/fftma_module/gen/lib_src/power.c @@ -1,10 +1,8 @@ #include "genlib.h" -#include "log.h" #include #include /*power covariance function*/ double power(double h, double alpha) { - log_info("RESULT = in progress, h = %f, alpha = %f", h, alpha); return pow(h, alpha); } diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c index 50ae486..02d9028 100644 --- a/fftma_module/gen/lib_src/prebuild_gwn.c +++ b/fftma_module/gen/lib_src/prebuild_gwn.c @@ -1,6 +1,4 @@ #include "geostat.h" -#include "log.h" -#include "memory.h" #include #include #include @@ -23,11 +21,6 @@ /*realization: structure defining a realization*/ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int i, j, k, maille0, maille1; int ntot; @@ -38,15 +31,6 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin if (*seed > 0.0) *seed = -(*seed); - log_info("RESULT = in progress, n[0] = %d, n[1] = %d, n[2] = %d, solver = %d", n[0], n[1], n[2], solver); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - ntot = n[0] * n[1] * n[2]; realization[0] = 0.; /*printf("Antes de llamar a chunkarray read\n"); @@ -76,23 +60,4 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin } } } - - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - - log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0); - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/ran2.c b/fftma_module/gen/lib_src/ran2.c index 93cd468..cd9dcdd 100644 --- a/fftma_module/gen/lib_src/ran2.c +++ b/fftma_module/gen/lib_src/ran2.c @@ -1,7 +1,5 @@ #include #include "genlib.h" -#include "log.h" -#include "memory.h" #define IM1 2147483563 #define IM2 2147483399 @@ -19,24 +17,10 @@ #define RNMX (1.0 - EPS) double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { - double* used_ram_t0 = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_t0); - - clock_t t = clock(); - int j; long k; double temp; - log_info("RESULT = in progress, idum = %f, idum2 = %f, iy = %f", *idum, *idum2, *iy); - - struct cpustat initial[cores]; - struct cpustat final[cores]; - - for (int i = 0; i < cores; i++) { - get_stats(&initial[i], i - 1); - } - if (*idum <= 0) { if (-(*idum) < 1) *idum = 1; @@ -68,29 +52,10 @@ double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { if (*iy < 1) (*iy) += IMM1; - t = clock() - t; - double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time - - for (int i = 0; i < cores; i++) { - get_stats(&final[i], i - 1); - } - - for (int i = 0; i < cores; i++) { - log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i])); - } - - double* used_ram_tf = malloc(sizeof(double)); - getVirtualMemUsed(used_ram_tf); - if ((temp = AM * (*iy)) > RNMX) { - log_info("RESULT = success, RNMX = %f, ELAPSED = %f, DIF USED VIRTUAL MEM = %5.1f MB", RNMX, time_taken, *used_ram_tf - *used_ram_t0); return (RNMX); } else { - log_info("RESULT = success, temp = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", temp, time_taken, *used_ram_tf - *used_ram_t0); return temp; } - - free(used_ram_t0); - free(used_ram_tf); } diff --git a/fftma_module/gen/lib_src/spherical.c b/fftma_module/gen/lib_src/spherical.c index a63c3ac..277fd67 100644 --- a/fftma_module/gen/lib_src/spherical.c +++ b/fftma_module/gen/lib_src/spherical.c @@ -1,12 +1,10 @@ #include "genlib.h" -#include "log.h" #include #include /*spherical covariance function*/ double spherical(double h) { double z; - log_info("RESULT = in progress, h = %f", h); if (h >= 1.) { z = 0.; @@ -14,7 +12,5 @@ double spherical(double h) { z = 1. - 1.5 * (double)h + 0.5 * (double)(h * h * h); } - log_info("RESULT = success, z = %f", z); - return z; } diff --git a/fftma_module/gen/moduleFFTMA.c b/fftma_module/gen/moduleFFTMA.c index bd64c3d..77337a6 100644 --- a/fftma_module/gen/moduleFFTMA.c +++ b/fftma_module/gen/moduleFFTMA.c @@ -7,7 +7,6 @@ #include "toolsFFTMA.h" #include "toolsFFTPSIM.h" #include "toolsIO.h" -#include "lib_src/log.h" #include "chunk_array.h" #include #include @@ -23,7 +22,7 @@ /* Y is the realization with mean and variance wanted */ static PyObject* genFunc(PyObject* self, PyObject* args) { - log_info("RESULT = in progress"); + //log_info("RESULT = in progress"); int n[3]; struct realization_mod Z, Y; struct grid_mod grid; @@ -77,7 +76,7 @@ static PyObject* genFunc(PyObject* self, PyObject* args) { /*chunk_array_free(Z.vector_2); remove("realization1.txt");*/ - log_info("RESULT = success"); + //log_info("RESULT = success"); return out_array; } diff --git a/fftma_module/gen/setup.py b/fftma_module/gen/setup.py index 9feb438..aa752ec 100644 --- a/fftma_module/gen/setup.py +++ b/fftma_module/gen/setup.py @@ -42,8 +42,6 @@ module_FFTMA = Extension( "./lib_src/clean_real.c", "./lib_src/testmemory.c", "./lib_src/genlib.c", - "./lib_src/log.c", - "./lib_src/memory.c", "./lib_src/chunk_array.c" ], ) From 0d40eae129630d03559a177c36f43ae84018ef5a Mon Sep 17 00:00:00 2001 From: chortas Date: Fri, 14 Jan 2022 22:40:26 -0300 Subject: [PATCH 02/15] Commit to debug --- fftma_module/gen/include/geostat.h | 4 +- fftma_module/gen/include/toolsFFTMA.h | 6 +- fftma_module/gen/lib_src/Py_kgeneration.c | 2 + fftma_module/gen/lib_src/build_real.c | 11 +- fftma_module/gen/lib_src/clean_real.c | 9 +- fftma_module/gen/lib_src/covariance.c | 29 ++- fftma_module/gen/lib_src/fftma2.c | 21 +-- fftma_module/gen/lib_src/fourt.c | 213 +++++++++++++++------- fftma_module/gen/lib_src/geostat.h | 4 +- fftma_module/gen/lib_src/prebuild_gwn.c | 11 +- 10 files changed, 206 insertions(+), 104 deletions(-) diff --git a/fftma_module/gen/include/geostat.h b/fftma_module/gen/include/geostat.h index 33cf755..f1d2922 100644 --- a/fftma_module/gen/include/geostat.h +++ b/fftma_module/gen/include/geostat.h @@ -327,7 +327,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid); /*variogram: structure defined above */ /*grid: structure defined above */ /*n: number of gridblocks along X,Y and Z*/ -void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); +void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ @@ -403,7 +403,7 @@ double exponential(double h); /*workr: utility real part vector for storage */ /*worki: utility imaginary part vector for storage */ /*The transformed data are returned in datar and datai*/ -void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); +void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); /*calculates F(x) = (1/a)*exp(-x*x/2)*/ double funtrun1(double x); diff --git a/fftma_module/gen/include/toolsFFTMA.h b/fftma_module/gen/include/toolsFFTMA.h index 1731ebb..42af6d4 100644 --- a/fftma_module/gen/include/toolsFFTMA.h +++ b/fftma_module/gen/include/toolsFFTMA.h @@ -51,7 +51,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r /* must be a Gaussian white noise */ /*realization: structure defining a realization*/ -void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed); +void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed); /* build_real */ /* build a realization in the spectral domain */ @@ -65,8 +65,8 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores); +void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores); -void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout); +void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores); #endif // define _TOOLSFFTMA_H diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c index c79ad46..a20469a 100644 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -27,6 +27,8 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, n[1] = 0; n[2] = 0; + printf("aaa"); + generate(&seed, N, Z, cores); /*FFTMA*/ diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c index f563f85..6274e3b 100644 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -1,4 +1,5 @@ #include "geostat.h" +#include "chunk_array.h" #include #include #include @@ -19,22 +20,26 @@ /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores) { +void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores) { int i, j, k, maille1; double temp; + chunk_array_read(realization); + /*decomposition and multiplication in the spectral domain*/ for (k = 1; k <= n[2]; k++) { for (j = 1; j <= n[1]; j++) { for (i = 1; i <= n[0]; i++) { maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; - temp = covar[maille1]; + chunk_aray_get(covar, maille1, &temp); if (temp > 0.) { temp = sqrt(temp) / (double)NTOT; } else if (temp < 0.) { temp = sqrt(-temp) / (double)NTOT; } - realization[maille1] *= temp; + double valuemaille1; + chunk_array_get(realization, maille1, &valuemaille1); + chunk_array_save(realization, maille1, temp * valuemaille1); ireal[maille1] *= temp; } } diff --git a/fftma_module/gen/lib_src/clean_real.c b/fftma_module/gen/lib_src/clean_real.c index 9917554..a6de2b7 100644 --- a/fftma_module/gen/lib_src/clean_real.c +++ b/fftma_module/gen/lib_src/clean_real.c @@ -1,4 +1,5 @@ #include "geostat.h" +#include "chunk_array.h" #include #include #include @@ -7,7 +8,7 @@ #include #include -void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout, int cores) { +void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores) { int i, j, k, maille0, maille1; double NTOT; @@ -15,6 +16,8 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, /*is the output realization already allocated?*/ /*if not, memory allocation*/ + chunk_array_read(vectorresult); + if (realout->vector == NULL || realout->n != realin->n) { realout->vector = (double*)malloc(realin->n * sizeof(double)); if (realout->vector == NULL) { @@ -31,7 +34,9 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; /* Modif du 18 juin 2003 */ /*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/ - realout->vector[maille0] = vectorresult[maille1]; + double valuemaille1; + chunk_array_get(vectorresult, maille1, &valuemaille1); + realout->vector[maille0] = valuemaille1; } } } diff --git a/fftma_module/gen/lib_src/covariance.c b/fftma_module/gen/lib_src/covariance.c index 1bc9064..97ef87a 100644 --- a/fftma_module/gen/lib_src/covariance.c +++ b/fftma_module/gen/lib_src/covariance.c @@ -1,15 +1,18 @@ #include "geostat.h" +#include "chunk_array.h" #include /*builds the sampled covariance function*/ /*dimensions are even*/ -void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { +void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { int i, j, k, maille, n2[3], symmetric; double di, dj, dk; for (i = 0; i < 3; i++) n2[i] = n[i] / 2; + chunk_array_read(covar); + for (i = 0; i <= n2[0]; i++) { for (j = 0; j <= n2[1]; j++) { for (k = 0; k <= n2[2]; k++) { @@ -19,12 +22,14 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, di = (double)i * mesh.DX; dj = (double)j * mesh.DY; dk = (double)k * mesh.DZ; - covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); + chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) { /*area 2*/ symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k)); - covar[symmetric] = covar[maille]; + double value; + chunk_array_get(covar, maille, &value); + chunk_array_save(covar, symmetric, value); } if (i > 0 && i < n2[0]) { @@ -33,13 +38,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, dj = (double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k); - covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); + chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); } if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { /*area 8*/ symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k)); - covar[symmetric] = covar[maille]; + double value; + chunk_array_get(covar, maille, &value); + chunk_array_save(covar, symmetric, value); } if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) { @@ -48,13 +55,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, dj = -(double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k); - covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); + chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); } if (k > 0 && k < n2[2]) { /*area 6*/ symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k)); - covar[symmetric] = covar[maille]; + double value; + chunk_array_get(covar, maille, &value); + chunk_array_save(covar, symmetric, value); } if (j > 0 && j < n2[1]) { @@ -63,13 +72,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, dj = -(double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + i + n[0] * (n[1] - j + n[1] * k); - covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); + chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); } if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { /*area 7*/ symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k)); - covar[symmetric] = covar[maille]; + double value; + chunk_array_get(covar, maille, &value); + chunk_array_save(covar, symmetric, value); } } } diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index d699498..5c8b224 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -1,4 +1,5 @@ #include "geostat.h" +#include "chunk_array.h" #include #include #include @@ -26,7 +27,10 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz; int solver; double temp; - double *ireal, *covar, *workr, *worki, *realization; + chunk_array_t *covar, *ireal, *realization; + double *workr, *worki; + + printf("Ayuda"); /*covariance axis normalization*/ axes(variogram.ap, variogram.scf, variogram.Nvario); @@ -50,14 +54,9 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r nxyz = NXYZ + 1; /*array initialization*/ - covar = (double*)malloc(ntot * sizeof(double)); - testmemory(covar); - - ireal = (double*)malloc(ntot * sizeof(double)); - testmemory(ireal); - - realization = (double*)malloc(ntot * sizeof(double)); - testmemory(realization); + covar = chunk_array_create("covar.txt", ntot, 512); + ireal = chunk_array_create("ireal.txt", ntot, 512); + realization = chunk_array_create("realization.txt", ntot, 512); workr = (double*)malloc(nmax * sizeof(double)); testmemory(workr); @@ -83,12 +82,12 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r /* build realization in spectral domain */ build_real(n, NTOT, covar, realization, ireal, cores); - free(covar); + chunk_array_free(covar); /*backward fourier transform*/ fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); - free(ireal); + chunk_array_free(ireal); free(workr); free(worki); diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index fc4184c..9d4ed0f 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -1,6 +1,7 @@ #include #include #include +#include "chunk_array.h" /*fast fourier transform */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ @@ -90,7 +91,7 @@ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* 10-06-2000, MLR */ -void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { +void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; @@ -99,6 +100,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp ntot *= nn[idim]; } + chunk_array_read(datar); + chunk_array_read(datai); + /*main loop for each dimension*/ np1 = 1; for (idim = 1; idim <= ndim; idim++) { @@ -184,8 +188,11 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp ntot /= 2; i = 1; for (j = 1; j <= ntot; j++) { - datar[j] = datar[i]; - datai[j] = datar[i + 1]; + double valuei, valuei1; + chunk_array_get(datar, i, &valuei); + chunk_array_get(datar, i+1, &valuei1); + chunk_array_save(datar, j, valuei); + chunk_array_save(datai, j, valuei1); i += 2; } @@ -202,12 +209,20 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp for (i1 = i2; i1 <= i1max; i1++) { for (i3 = i1; i3 <= ntot; i3 += np2) { j3 = j + i3 - i2; - tempr = datar[i3]; - tempi = datai[i3]; - datar[i3] = datar[j3]; - datai[i3] = datai[j3]; - datar[j3] = tempr; - datai[j3] = tempi; + double valueri3, valueii3, valuerj3, valueij3; + + chunk_array_get(datar, i3, &valueri3); + chunk_array_get(datai, i3, &valueii3); + chunk_array_get(datar, j3, &valuerj3); + chunk_array_get(datai, j3, &valueij3); + + tempr = valueri3; + tempi = valueii3; + + chunk_array_save(datar, i3, valuerj3); + chunk_array_save(datai, i3, valueij3); + chunk_array_save(datar, j3, tempr); + chunk_array_save(datai, j3, tempi); } } @@ -232,11 +247,14 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp for (i3 = i1; i3 <= ntot; i3 += np2) { j = i3; for (i = 1; i <= n; i++) { + double valuerj, valueij; + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); if (icase != 3) { - workr[i] = datar[j]; - worki[i] = datai[j]; + workr[i] = valuerj; + worki[i] = valueij; } else { - workr[i] = datar[j]; + workr[i] = valuerj; worki[i] = 0.; } ifp2 = np2; @@ -255,8 +273,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { - datar[i2] = workr[i]; - datai[i2] = worki[i]; + chunk_save(datar, i2, workr[i]); + chunk_save(datai, i2, worki[i]); i++; } } @@ -277,12 +295,16 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp L310: j = i1; for (i = imin; i <= ntot; i += istep) { - tempr = datar[i]; - tempi = datai[i]; - datar[i] = datar[j] - tempr; - datai[i] = datai[j] - tempi; - datar[j] = datar[j] + tempr; - datai[j] = datai[j] + tempi; + chunk_array_get(datar, i, &tempr); + chunk_array_get(datai, i, &tempi); + + double valuerj, valueij; + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); + chunk_array_save(datar, i, valuerj - tempr); + chunk_array_save(datai, i, valueij - tempi); + chunk_array_save(datar, j, valuerj + tempr); + chunk_array_save(datai, j, valueij + tempi); j += istep; } imin = 2 * imin - i1; @@ -301,16 +323,26 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { if (ifrwd != 0) { - tempr = datai[i]; - tempi = -datar[i]; + chunk_array_get(datai, i, &tempr); + chunk_array_get(datar, i, &tempi); + tempi = -tempi; } else { - tempr = -datai[i]; - tempi = datar[i]; + chunk_array_get(datai, i, &tempr); + chunk_array_get(datar, i, &tempi); } - datar[i] = datar[j] - tempr; - datai[i] = datai[j] - tempi; - datar[j] += tempr; - datai[j] += tempi; + double valuerj, valueij; + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); + + chunk_array_save(datar, i, valuerj - tempr); + chunk_array_save(datai, i, valueij - tempi); + + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); + + chunk_array_save(datar, j, valuerj + tempr); + chunk_array_save(datai, j, valueij + tempi); + j += istep; } @@ -347,12 +379,18 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp L510: j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { - tempr = datar[i] * wr - datai[i] * wi; - tempi = datar[i] * wi + datai[i] * wr; - datar[i] = datar[j] - tempr; - datai[i] = datai[j] - tempi; - datar[j] += tempr; - datai[j] += tempi; + double valueri, valueii, valuerj, valueij; + chunk_array_get(datar, i, &valueri); + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, i, &valueii); + chunk_array_get(datai, j, &valueij); + + tempr = valueri * wr - valueii * wi; + tempi = valueri * wi + valueii * wr; + chunk_data_save(datar, i, valuerj - tempr); + chunk_data_save(datai, i, valueij - tempi); + chunk_data_save(datar, j, valuerj + tempr); + chunk_data_save(datai, j, valueij + tempi); j += istep; } imin = 2 * imin - i1; @@ -405,23 +443,27 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { j = jmin + ifp2 - ifp1; - sr = datar[j]; - si = datai[j]; + chunk_array_get(datar, j, &sr); + chunk_array_get(datai, j, &si); oldsr = 0.; oldsi = 0.; j -= ifp1; L620: stmpr = sr; stmpi = si; - sr = twowr * sr - oldsr + datar[j]; - si = twowr * si - oldsi + datai[j]; + double valuerj, valueij; + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); + + sr = twowr * sr - oldsr + valuerj; + si = twowr * si - oldsi + valueij; oldsr = stmpr; oldsi = stmpi; j -= ifp1; if (j > jmin) goto L620; - workr[i] = wr * sr - wi * si - oldsr + datar[j]; - worki[i] = wi * sr + wr * si - oldsi + datai[j]; + workr[i] = wr * sr - wi * si - oldsr + valuerj; + worki[i] = wi * sr + wr * si - oldsi + valueij; jmin += ifp2; i++; } @@ -433,8 +475,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { - datar[j3] = workr[i]; - datai[j3] = worki[i]; + chunk_array_save(datar, j3, workr[i]); + chunk_array_save(datai, j3, worki[i]); i++; } } @@ -478,16 +520,22 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp L710: j = jmin; for (i = imin; i <= ntot; i += np2) { - sumr = (datar[i] + datar[j]) / 2.; - sumi = (datai[i] + datai[j]) / 2.; - difr = (datar[i] - datar[j]) / 2.; - difi = (datai[i] - datai[j]) / 2.; + double valueri, valueii, valuerj, valueij; + chunk_array_get(datar, i, &valueri); + chunk_array_get(datai, i, &valueii); + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); + + sumr = (valueri + valuerj) / 2.; + sumi = (valueii + valueij) / 2.; + difr = (valueri - valuerj) / 2.; + difi = (valueii - valueij) / 2.; tempr = wr * sumi + wi * difr; tempi = wi * sumi - wr * difr; - datar[i] = sumr + tempr; - datai[i] = difi + tempi; - datar[j] = sumr - tempr; - datai[j] = tempi - difi; + chunk_data_save(datar, i, sumr + tempr); + chunk_data_save(datai, i, difi + tempi); + chunk_data_save(datar, j, sumr - tempr); + chunk_data_save(datai, j, tempi - difi); j += np2; } imin++; @@ -504,7 +552,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp if (ifrwd == 0) goto L740; for (i = imin; i <= ntot; i += np2) { - datai[i] = -datai[i]; + double valueii; + chunk_data_get(datai, i, &valueii); + chunk_data_save(datai, i, -valueii); } L740: np2 *= 2; @@ -515,36 +565,56 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp imin = imax - nhalf; i = imin; goto L755; - L750: - datar[j] = datar[i]; - datai[j] = -datai[i]; + L750: ; + double valueri, valueii; + chunk_data_get(datar, i, &valueri); + chunk_data_get(datai, i, &valueii); + chunk_data_save(datar, j, valueri); + chunk_data_save(datai, j, -valueii); L755: i++; j--; if (i < imax) goto L750; - datar[j] = datar[imin] - datai[imin]; - datai[j] = 0.; + double valuerimin, valueiimin; + chunk_data_get(datar, imin, &valuerimin); + chunk_data_get(datai, imin, &valueiimin); + + chunk_data_save(datar, j, valuerimin - valueiimin); + chunk_data_save(datai, j, 0.); if (i >= j) { goto L780; } else { goto L770; } L765: - datar[j] = datar[i]; - datai[j] = datai[i]; + chunk_data_get(datar, i, &valueri); + chunk_data_get(datai, i, &valueii); + + chunk_data_save(datar, j, valueri); + chunk_data_save(datai, j, valueii); L770: i--; j--; if (i > imin) goto L765; - datar[j] = datar[imin] + datai[imin]; - datai[j] = 0.; + + chunk_data_get(datar, imin, &valuerimin); + chunk_data_get(datai, imin, &valueiimin); + + chunk_data_save(datar, j, valuerimin + valueiimin); + chunk_data_save(datai, j, 0.); + imax = imin; goto L745; - L780: - datar[1] += datai[1]; - datai[1] = 0.; + L780: ; + double valuei1, valuer1; + chunk_data_get(datai, 1, &valuei1); + chunk_data_get(datar, 1, &valuer1); + + + chunk_data_save(datar, 1, valuei1 + valuer1); + chunk_data_save(datai, 1, 0.); goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ @@ -562,15 +632,24 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp if (idim > 2) { j = jmax + np0; for (i = imin; i <= imax; i++) { - datar[i] = datar[j]; - datai[i] = -datai[j]; + double valuerj, valueij; + chunk_data_get(datar, j, &valuerj); + chunk_data_get(datai, j, &valueij); + + chunk_data_save(datar, i, valuerj); + chunk_data_save(datai, i, -valueij); j--; } } j = jmax; for (i = imin; i <= imax; i += np0) { - datar[i] = datar[j]; - datai[i] = -datai[j]; + double valuerj, valueij; + chunk_data_get(datar, j, &valuerj); + chunk_data_get(datai, j, &valueij); + + chunk_data_save(datar, i, valuerj); + chunk_data_save(datai, i, -valueij); + j -= np0; } } diff --git a/fftma_module/gen/lib_src/geostat.h b/fftma_module/gen/lib_src/geostat.h index 35d8bcc..c7244c9 100644 --- a/fftma_module/gen/lib_src/geostat.h +++ b/fftma_module/gen/lib_src/geostat.h @@ -328,7 +328,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid); /*variogram: structure defined above */ /*grid: structure defined above */ /*n: number of gridblocks along X,Y and Z*/ -void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); +void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ @@ -404,7 +404,7 @@ double exponential(double h); /*workr: utility real part vector for storage */ /*worki: utility imaginary part vector for storage */ /*The transformed data are returned in datar and datai*/ -void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); +void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); /*calculates F(x) = (1/a)*exp(-x*x/2)*/ double funtrun1(double x); diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c index 02d9028..a6382cf 100644 --- a/fftma_module/gen/lib_src/prebuild_gwn.c +++ b/fftma_module/gen/lib_src/prebuild_gwn.c @@ -1,4 +1,5 @@ #include "geostat.h" +#include "chunk_array.h" #include #include #include @@ -20,7 +21,7 @@ /* must be a Gaussian white noise */ /*realization: structure defining a realization*/ -void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed) { +void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed) { int i, j, k, maille0, maille1; int ntot; @@ -32,14 +33,14 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin *seed = -(*seed); ntot = n[0] * n[1] * n[2]; - realization[0] = 0.; + chunk_array_save(realization, 0, 0.); /*printf("Antes de llamar a chunkarray read\n"); chunk_array_read((*realin).vector_2); printf("Despues de llamar a chunkarray read\n");*/ if (solver == 1) { for (i = 0; i < ntot; i++) { double value = gasdev(seed, &idum2, &iy, iv, cores); - realization[i+1] = value; + chunk_array_save(realization, i+1, value); //chunk_array_get((*realin).vector_2, i, &realization[i + 1]); } } else { @@ -51,10 +52,10 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin //maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; //printf("Maille0 es %d", maille0); double value = gasdev(seed, &idum2, &iy, iv, cores); - realization[maille1] = value; + chunk_array_save(realization, maille1, value); //chunk_array_get((*realin).vector_2, maille0, &realization[maille1]); } else { - realization[maille1] = 0.; + chunk_array_save(realization, maille1, 0.); } } } From dcbaafdb1a301406f2bb6489a2192a261545ea10 Mon Sep 17 00:00:00 2001 From: chortas Date: Fri, 14 Jan 2022 23:02:22 -0300 Subject: [PATCH 03/15] At least compiling. Segfault --- fftma_module/gen/include/toolsFFTMA.h | 1 + fftma_module/gen/lib_src/build_real.c | 4 +- fftma_module/gen/lib_src/fourt.c | 80 +++++++++++++-------------- 3 files changed, 43 insertions(+), 42 deletions(-) diff --git a/fftma_module/gen/include/toolsFFTMA.h b/fftma_module/gen/include/toolsFFTMA.h index 42af6d4..ed7a6d1 100644 --- a/fftma_module/gen/include/toolsFFTMA.h +++ b/fftma_module/gen/include/toolsFFTMA.h @@ -1,5 +1,6 @@ #include "genlib.h" #include "geostat.h" +#include "chunk_array.h" #include #include diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c index 6274e3b..e46dd5f 100644 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -1,5 +1,4 @@ #include "geostat.h" -#include "chunk_array.h" #include #include #include @@ -7,6 +6,7 @@ #include #include #include +#include "chunk_array.h" /* build_real */ /* build a realization in the spectral domain */ @@ -31,7 +31,7 @@ void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realiza for (j = 1; j <= n[1]; j++) { for (i = 1; i <= n[0]; i++) { maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; - chunk_aray_get(covar, maille1, &temp); + chunk_array_get(covar, maille1, &temp); if (temp > 0.) { temp = sqrt(temp) / (double)NTOT; } else if (temp < 0.) { diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index 9d4ed0f..4d015b6 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -273,8 +273,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { - chunk_save(datar, i2, workr[i]); - chunk_save(datai, i2, worki[i]); + chunk_array_save(datar, i2, workr[i]); + chunk_array_save(datai, i2, worki[i]); i++; } } @@ -387,10 +387,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; - chunk_data_save(datar, i, valuerj - tempr); - chunk_data_save(datai, i, valueij - tempi); - chunk_data_save(datar, j, valuerj + tempr); - chunk_data_save(datai, j, valueij + tempi); + chunk_array_save(datar, i, valuerj - tempr); + chunk_array_save(datai, i, valueij - tempi); + chunk_array_save(datar, j, valuerj + tempr); + chunk_array_save(datai, j, valueij + tempi); j += istep; } imin = 2 * imin - i1; @@ -532,10 +532,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int difi = (valueii - valueij) / 2.; tempr = wr * sumi + wi * difr; tempi = wi * sumi - wr * difr; - chunk_data_save(datar, i, sumr + tempr); - chunk_data_save(datai, i, difi + tempi); - chunk_data_save(datar, j, sumr - tempr); - chunk_data_save(datai, j, tempi - difi); + chunk_array_save(datar, i, sumr + tempr); + chunk_array_save(datai, i, difi + tempi); + chunk_array_save(datar, j, sumr - tempr); + chunk_array_save(datai, j, tempi - difi); j += np2; } imin++; @@ -553,8 +553,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int goto L740; for (i = imin; i <= ntot; i += np2) { double valueii; - chunk_data_get(datai, i, &valueii); - chunk_data_save(datai, i, -valueii); + chunk_array_get(datai, i, &valueii); + chunk_array_save(datai, i, -valueii); } L740: np2 *= 2; @@ -567,54 +567,54 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int goto L755; L750: ; double valueri, valueii; - chunk_data_get(datar, i, &valueri); - chunk_data_get(datai, i, &valueii); - chunk_data_save(datar, j, valueri); - chunk_data_save(datai, j, -valueii); + chunk_array_get(datar, i, &valueri); + chunk_array_get(datai, i, &valueii); + chunk_array_save(datar, j, valueri); + chunk_array_save(datai, j, -valueii); L755: i++; j--; if (i < imax) goto L750; double valuerimin, valueiimin; - chunk_data_get(datar, imin, &valuerimin); - chunk_data_get(datai, imin, &valueiimin); + chunk_array_get(datar, imin, &valuerimin); + chunk_array_get(datai, imin, &valueiimin); - chunk_data_save(datar, j, valuerimin - valueiimin); - chunk_data_save(datai, j, 0.); + chunk_array_save(datar, j, valuerimin - valueiimin); + chunk_array_save(datai, j, 0.); if (i >= j) { goto L780; } else { goto L770; } L765: - chunk_data_get(datar, i, &valueri); - chunk_data_get(datai, i, &valueii); + chunk_array_get(datar, i, &valueri); + chunk_array_get(datai, i, &valueii); - chunk_data_save(datar, j, valueri); - chunk_data_save(datai, j, valueii); + chunk_array_save(datar, j, valueri); + chunk_array_save(datai, j, valueii); L770: i--; j--; if (i > imin) goto L765; - chunk_data_get(datar, imin, &valuerimin); - chunk_data_get(datai, imin, &valueiimin); + chunk_array_get(datar, imin, &valuerimin); + chunk_array_get(datai, imin, &valueiimin); - chunk_data_save(datar, j, valuerimin + valueiimin); - chunk_data_save(datai, j, 0.); + chunk_array_save(datar, j, valuerimin + valueiimin); + chunk_array_save(datai, j, 0.); imax = imin; goto L745; L780: ; double valuei1, valuer1; - chunk_data_get(datai, 1, &valuei1); - chunk_data_get(datar, 1, &valuer1); + chunk_array_get(datai, 1, &valuei1); + chunk_array_get(datar, 1, &valuer1); - chunk_data_save(datar, 1, valuei1 + valuer1); - chunk_data_save(datai, 1, 0.); + chunk_array_save(datar, 1, valuei1 + valuer1); + chunk_array_save(datai, 1, 0.); goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ @@ -633,22 +633,22 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j = jmax + np0; for (i = imin; i <= imax; i++) { double valuerj, valueij; - chunk_data_get(datar, j, &valuerj); - chunk_data_get(datai, j, &valueij); + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); - chunk_data_save(datar, i, valuerj); - chunk_data_save(datai, i, -valueij); + chunk_array_save(datar, i, valuerj); + chunk_array_save(datai, i, -valueij); j--; } } j = jmax; for (i = imin; i <= imax; i += np0) { double valuerj, valueij; - chunk_data_get(datar, j, &valuerj); - chunk_data_get(datai, j, &valueij); + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); - chunk_data_save(datar, i, valuerj); - chunk_data_save(datai, i, -valueij); + chunk_array_save(datar, i, valuerj); + chunk_array_save(datai, i, -valueij); j -= np0; } From c64a61e0718d9a5df142c1be52d7b96eed44920b Mon Sep 17 00:00:00 2001 From: chortas Date: Sat, 15 Jan 2022 11:12:00 -0300 Subject: [PATCH 04/15] Add logs to debug and fix seg fault --- fftma_module/gen/include/toolsFFTMA.h | 2 +- fftma_module/gen/lib_src/Py_kgeneration.c | 2 - fftma_module/gen/lib_src/build_real.c | 15 +++-- fftma_module/gen/lib_src/chunk_array.c | 4 +- fftma_module/gen/lib_src/fftma2.c | 17 ++++- fftma_module/gen/lib_src/fourt.c | 81 ++++++++++++++++------- fftma_module/gen/moduleFFTMA.c | 20 +----- 7 files changed, 88 insertions(+), 53 deletions(-) diff --git a/fftma_module/gen/include/toolsFFTMA.h b/fftma_module/gen/include/toolsFFTMA.h index ed7a6d1..4dfe82e 100644 --- a/fftma_module/gen/include/toolsFFTMA.h +++ b/fftma_module/gen/include/toolsFFTMA.h @@ -66,7 +66,7 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores); +void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, chunk_array_t* ireal, int cores); void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores); diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c index a20469a..c79ad46 100644 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -27,8 +27,6 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, n[1] = 0; n[2] = 0; - printf("aaa"); - generate(&seed, N, Z, cores); /*FFTMA*/ diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c index e46dd5f..9d31bed 100644 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -20,11 +20,13 @@ /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores) { +void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, chunk_array_t* ireal, int cores) { int i, j, k, maille1; double temp; chunk_array_read(realization); + chunk_array_read(ireal); + chunk_array_read(covar); /*decomposition and multiplication in the spectral domain*/ for (k = 1; k <= n[2]; k++) { @@ -37,10 +39,13 @@ void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realiza } else if (temp < 0.) { temp = sqrt(-temp) / (double)NTOT; } - double valuemaille1; - chunk_array_get(realization, maille1, &valuemaille1); - chunk_array_save(realization, maille1, temp * valuemaille1); - ireal[maille1] *= temp; + double valuerealizationmaille1, valueirealmaille1; + + chunk_array_get(realization, maille1, &valuerealizationmaille1); + chunk_array_get(ireal, maille1, &valueirealmaille1); + + chunk_array_save(realization, maille1, temp * valuerealizationmaille1); + chunk_array_save(ireal, maille1, temp * valueirealmaille1); } } } diff --git a/fftma_module/gen/lib_src/chunk_array.c b/fftma_module/gen/lib_src/chunk_array.c index 806ec88..81932a0 100644 --- a/fftma_module/gen/lib_src/chunk_array.c +++ b/fftma_module/gen/lib_src/chunk_array.c @@ -8,7 +8,6 @@ void chunk_array_free(chunk_array_t* chunk_array) { } bool chunk_array_update_read(chunk_array_t* chunk_array) { - printf("Llame a chunk array update\n"); size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); chunk_array->init_pos += newLen; } @@ -52,9 +51,12 @@ chunk_array_t* chunk_array_create(char* filename, size_t total_size, size_t chun } void chunk_array_read(chunk_array_t* chunk_array) { + printf("chunk_array 54\n"); rewind(chunk_array->fp); chunk_array->init_pos = 0; + printf("chunk_array 57\n"); size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); + printf("chunk_array 59\n"); } void chunk_array_write(chunk_array_t* chunk_array, char* filename) { diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 5c8b224..6659c69 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -30,8 +30,6 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r chunk_array_t *covar, *ireal, *realization; double *workr, *worki; - printf("Ayuda"); - /*covariance axis normalization*/ axes(variogram.ap, variogram.scf, variogram.Nvario); @@ -64,11 +62,15 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r worki = (double*)malloc(nmax * sizeof(double)); testmemory(worki); + printf("pre covariance\n"); /*covariance function creation*/ covariance(covar, variogram, grid, n, cores); + printf("post covariance\n"); /*power spectrum*/ + printf("pre fourt1\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); + printf("post fourt1\n"); /*organization of the input Gaussian white noise*/ printf("pre prebuild_gwn\n"); @@ -76,21 +78,32 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r prebuild_gwn(grid, n, realin, realization, solver, cores, seed); printf("post prebuild_gwn\n"); + printf("pre fourt2\n"); /*forward fourier transform of the GWN*/ fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); + printf("post fourt2\n"); + printf("pre build_real\n"); /* build realization in spectral domain */ build_real(n, NTOT, covar, realization, ireal, cores); + printf("post build_real\n"); chunk_array_free(covar); + remove("covar.txt"); + printf("pre fourt3\n"); /*backward fourier transform*/ fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); + printf("post fourt3\n"); chunk_array_free(ireal); + remove("ireal.txt"); + free(workr); free(worki); + printf("pre clean_real\n"); /*output realization*/ clean_real(realin, n, grid, realization, realout, cores); + printf("post clean_real\n"); } diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index 4d015b6..e438cec 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -94,15 +94,21 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; - + double valueri, valueri1, valueri3, valueii3, valuerj3, valueij3, valuerj, valueij, valueii, valuerimin, valueiimin, valuei1, valuer1; + ntot = 1; for (idim = 0; idim < ndim; idim++) { ntot *= nn[idim]; } + //printf("104\n"); + chunk_array_read(datar); + //printf("107\n"); chunk_array_read(datai); + //printf("110\n"); + /*main loop for each dimension*/ np1 = 1; for (idim = 1; idim <= ndim; idim++) { @@ -188,11 +194,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ntot /= 2; i = 1; for (j = 1; j <= ntot; j++) { - double valuei, valuei1; - chunk_array_get(datar, i, &valuei); - chunk_array_get(datar, i+1, &valuei1); - chunk_array_save(datar, j, valuei); - chunk_array_save(datai, j, valuei1); + //printf("196\n"); + chunk_array_get(datar, i, &valueri); + chunk_array_get(datar, i+1, &valueri1); + chunk_array_save(datar, j, valueri); + chunk_array_save(datai, j, valueri1); + //printf("201\n"); i += 2; } @@ -209,20 +216,23 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (i1 = i2; i1 <= i1max; i1++) { for (i3 = i1; i3 <= ntot; i3 += np2) { j3 = j + i3 - i2; - double valueri3, valueii3, valuerj3, valueij3; - + + //printf("219\n"); chunk_array_get(datar, i3, &valueri3); chunk_array_get(datai, i3, &valueii3); chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); + //printf("224\n"); tempr = valueri3; tempi = valueii3; + //printf("229\n"); chunk_array_save(datar, i3, valuerj3); chunk_array_save(datai, i3, valueij3); chunk_array_save(datar, j3, tempr); chunk_array_save(datai, j3, tempi); + //printf("234\n"); } } @@ -247,9 +257,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (i3 = i1; i3 <= ntot; i3 += np2) { j = i3; for (i = 1; i <= n; i++) { - double valuerj, valueij; + //printf("259\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + //printf("262\n"); if (icase != 3) { workr[i] = valuerj; worki[i] = valueij; @@ -273,8 +284,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { + //printf("286\n"); chunk_array_save(datar, i2, workr[i]); chunk_array_save(datai, i2, worki[i]); + //printf("289\n"); i++; } } @@ -295,16 +308,17 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L310: j = i1; for (i = imin; i <= ntot; i += istep) { + //printf("310\n"); chunk_array_get(datar, i, &tempr); chunk_array_get(datai, i, &tempi); - double valuerj, valueij; chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); + //printf("320\n"); j += istep; } imin = 2 * imin - i1; @@ -323,14 +337,19 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { if (ifrwd != 0) { + //printf("339\n"); chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); + //printf("342\n"); tempi = -tempi; } else { + //printf("345\n"); chunk_array_get(datai, i, &tempr); + tempr = -tempr; chunk_array_get(datar, i, &tempi); + //printf("349\n"); } - double valuerj, valueij; + //printf("351\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); @@ -342,7 +361,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - + //printf("363\n"); j += istep; } @@ -379,7 +398,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L510: j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { - double valueri, valueii, valuerj, valueij; + //printf("400\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, i, &valueii); @@ -391,6 +410,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); + //printf("412\n"); j += istep; } imin = 2 * imin - i1; @@ -443,17 +463,20 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { j = jmin + ifp2 - ifp1; + //printf("465\n"); chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); + //printf("468\n"); oldsr = 0.; oldsi = 0.; j -= ifp1; L620: stmpr = sr; stmpi = si; - double valuerj, valueij; + //printf("475\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + //printf("478\n"); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; @@ -475,8 +498,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { + //printf("500\n"); chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); + //printf("503\n"); i++; } } @@ -520,7 +545,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L710: j = jmin; for (i = imin; i <= ntot; i += np2) { - double valueri, valueii, valuerj, valueij; + //printf("547\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); chunk_array_get(datar, j, &valuerj); @@ -536,6 +561,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datai, i, difi + tempi); chunk_array_save(datar, j, sumr - tempr); chunk_array_save(datai, j, tempi - difi); + //printf("563\n"); j += np2; } imin++; @@ -552,9 +578,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int if (ifrwd == 0) goto L740; for (i = imin; i <= ntot; i += np2) { - double valueii; + //printf("580\n"); chunk_array_get(datai, i, &valueii); chunk_array_save(datai, i, -valueii); + //printf("583\n"); } L740: np2 *= 2; @@ -566,55 +593,61 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i = imin; goto L755; L750: ; - double valueri, valueii; + //printf("595\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); + //printf("600\n"); L755: i++; j--; if (i < imax) goto L750; - double valuerimin, valueiimin; + //printf("606\n"); chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); + //printf("612\n"); if (i >= j) { goto L780; } else { goto L770; } L765: + //printf("619\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); + //printf("625\n"); L770: i--; j--; if (i > imin) goto L765; - + + //printf("632\n"); chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); + //printf("638\n"); imax = imin; goto L745; L780: ; - double valuei1, valuer1; + //printf("643\n"); chunk_array_get(datai, 1, &valuei1); chunk_array_get(datar, 1, &valuer1); - chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); + //printf("649\n"); goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ @@ -632,23 +665,25 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int if (idim > 2) { j = jmax + np0; for (i = imin; i <= imax; i++) { - double valuerj, valueij; + //printf("667\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); + //printf("673\n"); j--; } } j = jmax; for (i = imin; i <= imax; i += np0) { - double valuerj, valueij; + //printf("679\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); + //printf("685\n"); j -= np0; } diff --git a/fftma_module/gen/moduleFFTMA.c b/fftma_module/gen/moduleFFTMA.c index 77337a6..e7ead11 100644 --- a/fftma_module/gen/moduleFFTMA.c +++ b/fftma_module/gen/moduleFFTMA.c @@ -36,7 +36,6 @@ static PyObject* genFunc(PyObject* self, PyObject* args) { if (!Py_getvalues(args, &seed, &grid, &variogram, &stat, &cores)) return NULL; - printf("deberia imprimir\n"); Py_kgeneration(seed, grid, stat, variogram, &Z, &Y, n, cores); out_dims[0] = grid.NZ; @@ -49,34 +48,17 @@ static PyObject* genFunc(PyObject* self, PyObject* args) { PyArray_ENABLEFLAGS(out_array, NPY_ARRAY_OWNDATA); - printf("me localizo\n"); - free(stat.mean); free(stat.variance); - - printf("termino stat\n"); - free(variogram.var); - - printf("1\n"); - printf("2\n"); free(variogram.alpha); - printf("3\n"); free(variogram.scf); - printf("4\n"); free(variogram.ap); - - //free(variogram.vario); - - - printf("Termino variogram\n"); - - printf("aca no deberia llegar\n"); + free(variogram.vario); /*chunk_array_free(Z.vector_2); remove("realization1.txt");*/ - //log_info("RESULT = success"); return out_array; } From 7ef3c2d2e5cc4485023d442bcd5bd475639d7acd Mon Sep 17 00:00:00 2001 From: chortas Date: Sat, 15 Jan 2022 17:40:22 -0300 Subject: [PATCH 05/15] Debugging fourt --- fftma_module/gen/include/chunk_array.h | 2 +- fftma_module/gen/lib_src/fftma2.c | 8 +- fftma_module/gen/lib_src/fourt.c | 119 +++++++++++++++---------- fftma_module/gen/test.py | 5 +- 4 files changed, 80 insertions(+), 54 deletions(-) diff --git a/fftma_module/gen/include/chunk_array.h b/fftma_module/gen/include/chunk_array.h index 2705b81..8190a89 100644 --- a/fftma_module/gen/include/chunk_array.h +++ b/fftma_module/gen/include/chunk_array.h @@ -7,7 +7,7 @@ #include #include -#define MAX_CHUNK_SIZE 512 +#define MAX_CHUNK_SIZE 1500 typedef struct chunk_array { size_t init_pos; diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 6659c69..7204d50 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -51,10 +51,12 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r NXYZ = grid.NX * grid.NY * grid.NZ; nxyz = NXYZ + 1; + printf("ntot es %d\n", ntot); + /*array initialization*/ - covar = chunk_array_create("covar.txt", ntot, 512); - ireal = chunk_array_create("ireal.txt", ntot, 512); - realization = chunk_array_create("realization.txt", ntot, 512); + covar = chunk_array_create("covar.txt", ntot, 1500); + ireal = chunk_array_create("ireal.txt", ntot, 1500); + realization = chunk_array_create("realization.txt", ntot, 1500); workr = (double*)malloc(nmax * sizeof(double)); testmemory(workr); diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index e438cec..2d9af7e 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -101,14 +101,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ntot *= nn[idim]; } - //printf("104\n"); - chunk_array_read(datar); - //printf("107\n"); chunk_array_read(datai); - //printf("110\n"); - /*main loop for each dimension*/ np1 = 1; for (idim = 1; idim <= ndim; idim++) { @@ -194,12 +189,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ntot /= 2; i = 1; for (j = 1; j <= ntot; j++) { - //printf("196\n"); chunk_array_get(datar, i, &valueri); + printf("datar[i] = %f\n", valueri); chunk_array_get(datar, i+1, &valueri1); + printf("datar[i1] = %f\n", valueri1); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); - //printf("201\n"); i += 2; } @@ -217,22 +212,23 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (i3 = i1; i3 <= ntot; i3 += np2) { j3 = j + i3 - i2; - //printf("219\n"); chunk_array_get(datar, i3, &valueri3); chunk_array_get(datai, i3, &valueii3); chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); - //printf("224\n"); + + printf("datar[i3] = %f\n", valueri3); + printf("datai[i3] = %f\n", valueii3); + printf("datar[j3] = %f\n", valuerj3); + printf("datai[j3] = %f\n", valueij3); tempr = valueri3; tempi = valueii3; - //printf("229\n"); chunk_array_save(datar, i3, valuerj3); chunk_array_save(datai, i3, valueij3); chunk_array_save(datar, j3, tempr); chunk_array_save(datai, j3, tempi); - //printf("234\n"); } } @@ -257,10 +253,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (i3 = i1; i3 <= ntot; i3 += np2) { j = i3; for (i = 1; i <= n; i++) { - //printf("259\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - //printf("262\n"); + + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + if (icase != 3) { workr[i] = valuerj; worki[i] = valueij; @@ -284,10 +282,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { - //printf("286\n"); chunk_array_save(datar, i2, workr[i]); chunk_array_save(datai, i2, worki[i]); - //printf("289\n"); i++; } } @@ -308,17 +304,20 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L310: j = i1; for (i = imin; i <= ntot; i += istep) { - //printf("310\n"); chunk_array_get(datar, i, &tempr); chunk_array_get(datai, i, &tempi); - chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + + printf("datar[i] = %f\n", tempr); + printf("datai[i] = %f\n", tempi); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - //printf("320\n"); j += istep; } imin = 2 * imin - i1; @@ -337,31 +336,37 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { if (ifrwd != 0) { - //printf("339\n"); chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); - //printf("342\n"); + + printf("datai[i] = %f\n", tempr); + printf("datar[i] = %f\n", tempi); + tempi = -tempi; } else { - //printf("345\n"); chunk_array_get(datai, i, &tempr); + printf("datai[i] = %f\n", tempr); tempr = -tempr; chunk_array_get(datar, i, &tempi); - //printf("349\n"); + printf("datar[i] = %f\n", tempi); } - //printf("351\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - //printf("363\n"); j += istep; } @@ -398,19 +403,22 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L510: j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { - //printf("400\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, i, &valueii); chunk_array_get(datai, j, &valueij); + printf("datar[i] = %f\n", valueri); + printf("datar[j] = %f\n", valuerj); + printf("datai[i] = %f\n", valueii); + printf("datai[j] = %f\n", valueij); + tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - //printf("412\n"); j += istep; } imin = 2 * imin - i1; @@ -463,20 +471,24 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { j = jmin + ifp2 - ifp1; - //printf("465\n"); chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); - //printf("468\n"); + + printf("datar[j] = %f\n", sr); + printf("datai[j] = %f\n", si); + oldsr = 0.; oldsi = 0.; j -= ifp1; L620: stmpr = sr; stmpi = si; - //printf("475\n"); + chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - //printf("478\n"); + + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; @@ -498,10 +510,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { - //printf("500\n"); chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); - //printf("503\n"); i++; } } @@ -545,12 +555,16 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int L710: j = jmin; for (i = imin; i <= ntot; i += np2) { - //printf("547\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + printf("datar[i] = %f\n", valueri); + printf("datai[i] = %f\n", valueii); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + sumr = (valueri + valuerj) / 2.; sumi = (valueii + valueij) / 2.; difr = (valueri - valuerj) / 2.; @@ -578,10 +592,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int if (ifrwd == 0) goto L740; for (i = imin; i <= ntot; i += np2) { - //printf("580\n"); chunk_array_get(datai, i, &valueii); + printf("datai[i] = %f\n", valueii); chunk_array_save(datai, i, -valueii); - //printf("583\n"); } L740: np2 *= 2; @@ -593,61 +606,67 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i = imin; goto L755; L750: ; - //printf("595\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); + + printf("datar[i] = %f\n", valueri); + printf("datai[i] = %f\n", valueii); + chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); - //printf("600\n"); L755: i++; j--; if (i < imax) goto L750; - //printf("606\n"); chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); + printf("datar[imin] = %f\n", valuerimin); + printf("datai[imin] = %f\n", valueiimin); + chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); - //printf("612\n"); if (i >= j) { goto L780; } else { goto L770; } L765: - //printf("619\n"); chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); + printf("datar[i] = %f\n", valueri); + printf("datai[i] = %f\n", valueii); + chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); - //printf("625\n"); L770: i--; j--; if (i > imin) goto L765; - //printf("632\n"); chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); + printf("datar[imin] = %f\n", valuerimin); + printf("datai[imin] = %f\n", valueiimin); + chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); - //printf("638\n"); imax = imin; goto L745; L780: ; - //printf("643\n"); chunk_array_get(datai, 1, &valuei1); chunk_array_get(datar, 1, &valuer1); + printf("datai[1] = %f\n", valuei1); + printf("datar[1] = %f\n", valuer1); + chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); - //printf("649\n"); goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ @@ -665,25 +684,27 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int if (idim > 2) { j = jmax + np0; for (i = imin; i <= imax; i++) { - //printf("667\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); - //printf("673\n"); j--; } } j = jmax; for (i = imin; i <= imax; i += np0) { - //printf("679\n"); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); + printf("datar[j] = %f\n", valuerj); + printf("datai[j] = %f\n", valueij); + chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); - //printf("685\n"); j -= np0; } diff --git a/fftma_module/gen/test.py b/fftma_module/gen/test.py index 7ca3b6d..660a103 100644 --- a/fftma_module/gen/test.py +++ b/fftma_module/gen/test.py @@ -29,6 +29,9 @@ variance=3.5682389 typ=3 k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8) -np.save(f"out_{N}.npy",k) + +#np.save(f"out_{N}.npy",k) +#k2 = np.load("out_2_16.npy") + del k gc.collect() \ No newline at end of file From 4a90597609438aeb894dfd03e9e304c306fd875c Mon Sep 17 00:00:00 2001 From: chortas Date: Sat, 15 Jan 2022 18:24:03 -0300 Subject: [PATCH 06/15] Add numeration to prints and remove extra prints --- fftma_module/gen/lib_src/Py_kgeneration.c | 2 - fftma_module/gen/lib_src/chunk_array.c | 4 - fftma_module/gen/lib_src/fftma2.c | 16 ---- fftma_module/gen/lib_src/fourt.c | 95 +++++++++++------------ 4 files changed, 47 insertions(+), 70 deletions(-) diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c index c79ad46..c633ce9 100644 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -30,9 +30,7 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, generate(&seed, N, Z, cores); /*FFTMA*/ - printf("pre fftma2\n"); FFTMA2(variogram, grid, n, Z, Y, cores, &seed); - printf("post fftma2\n"); /* make a log normal realization */ if (stat.type == 1 || stat.type == 2) { diff --git a/fftma_module/gen/lib_src/chunk_array.c b/fftma_module/gen/lib_src/chunk_array.c index 81932a0..58d0c74 100644 --- a/fftma_module/gen/lib_src/chunk_array.c +++ b/fftma_module/gen/lib_src/chunk_array.c @@ -51,18 +51,14 @@ chunk_array_t* chunk_array_create(char* filename, size_t total_size, size_t chun } void chunk_array_read(chunk_array_t* chunk_array) { - printf("chunk_array 54\n"); rewind(chunk_array->fp); chunk_array->init_pos = 0; - printf("chunk_array 57\n"); size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); - printf("chunk_array 59\n"); } void chunk_array_write(chunk_array_t* chunk_array, char* filename) { chunk_array->fp = fopen(filename, "w"); if (chunk_array->fp == NULL) { - printf("ke"); fclose(chunk_array->fp); chunk_array->fp = fopen(filename, "w"); } diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 7204d50..dcca2ff 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -51,8 +51,6 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r NXYZ = grid.NX * grid.NY * grid.NZ; nxyz = NXYZ + 1; - printf("ntot es %d\n", ntot); - /*array initialization*/ covar = chunk_array_create("covar.txt", ntot, 1500); ireal = chunk_array_create("ireal.txt", ntot, 1500); @@ -64,39 +62,27 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r worki = (double*)malloc(nmax * sizeof(double)); testmemory(worki); - printf("pre covariance\n"); /*covariance function creation*/ covariance(covar, variogram, grid, n, cores); - printf("post covariance\n"); /*power spectrum*/ - printf("pre fourt1\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); - printf("post fourt1\n"); /*organization of the input Gaussian white noise*/ - printf("pre prebuild_gwn\n"); solver = 0; prebuild_gwn(grid, n, realin, realization, solver, cores, seed); - printf("post prebuild_gwn\n"); - printf("pre fourt2\n"); /*forward fourier transform of the GWN*/ fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); - printf("post fourt2\n"); - printf("pre build_real\n"); /* build realization in spectral domain */ build_real(n, NTOT, covar, realization, ireal, cores); - printf("post build_real\n"); chunk_array_free(covar); remove("covar.txt"); - printf("pre fourt3\n"); /*backward fourier transform*/ fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); - printf("post fourt3\n"); chunk_array_free(ireal); remove("ireal.txt"); @@ -104,8 +90,6 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r free(workr); free(worki); - printf("pre clean_real\n"); /*output realization*/ clean_real(realin, n, grid, realization, realout, cores); - printf("post clean_real\n"); } diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index 2d9af7e..c340173 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -190,9 +190,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i = 1; for (j = 1; j <= ntot; j++) { chunk_array_get(datar, i, &valueri); - printf("datar[i] = %f\n", valueri); + printf("[1] datar[i] = %f\n", valueri); chunk_array_get(datar, i+1, &valueri1); - printf("datar[i1] = %f\n", valueri1); + printf("[2] datar[i1] = %f\n", valueri1); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); i += 2; @@ -217,10 +217,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); - printf("datar[i3] = %f\n", valueri3); - printf("datai[i3] = %f\n", valueii3); - printf("datar[j3] = %f\n", valuerj3); - printf("datai[j3] = %f\n", valueij3); + printf("[3] datar[i3] = %f\n", valueri3); + printf("[4] datai[i3] = %f\n", valueii3); + printf("[5] datar[j3] = %f\n", valuerj3); + printf("[6] datai[j3] = %f\n", valueij3); tempr = valueri3; tempi = valueii3; @@ -256,8 +256,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[7] datar[j] = %f\n", valuerj); + printf("[8] datai[j] = %f\n", valueij); if (icase != 3) { workr[i] = valuerj; @@ -309,10 +309,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[i] = %f\n", tempr); - printf("datai[i] = %f\n", tempi); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[9] datar[i] = %f\n", tempr); + printf("[10] datai[i] = %f\n", tempi); + printf("[11] datar[j] = %f\n", valuerj); + printf("[12] datai[j] = %f\n", valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); @@ -339,22 +339,22 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); - printf("datai[i] = %f\n", tempr); - printf("datar[i] = %f\n", tempi); + printf("[13] datai[i] = %f\n", tempr); + printf("[14] datar[i] = %f\n", tempi); tempi = -tempi; } else { chunk_array_get(datai, i, &tempr); - printf("datai[i] = %f\n", tempr); + printf("[15] datai[i] = %f\n", tempr); tempr = -tempr; chunk_array_get(datar, i, &tempi); - printf("datar[i] = %f\n", tempi); + printf("[16] datar[i] = %f\n", tempi); } chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[17] datar[j] = %f\n", valuerj); + printf("[18] datai[j] = %f\n", valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); @@ -362,8 +362,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[19] datar[j] = %f\n", valuerj); + printf("[20] datai[j] = %f\n", valueij); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); @@ -408,10 +408,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &valueii); chunk_array_get(datai, j, &valueij); - printf("datar[i] = %f\n", valueri); - printf("datar[j] = %f\n", valuerj); - printf("datai[i] = %f\n", valueii); - printf("datai[j] = %f\n", valueij); + printf("[21] datar[i] = %f\n", valueri); + printf("[22] datar[j] = %f\n", valuerj); + printf("[23] datai[i] = %f\n", valueii); + printf("[24] datai[j] = %f\n", valueij); tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; @@ -474,8 +474,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); - printf("datar[j] = %f\n", sr); - printf("datai[j] = %f\n", si); + printf("[25] datar[j] = %f\n", sr); + printf("[26] datai[j] = %f\n", si); oldsr = 0.; oldsi = 0.; @@ -487,8 +487,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[27] datar[j] = %f\n", valuerj); + printf("[28] datai[j] = %f\n", valueij); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; @@ -560,10 +560,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[i] = %f\n", valueri); - printf("datai[i] = %f\n", valueii); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[29] datar[i] = %f\n", valueri); + printf("[30] datai[i] = %f\n", valueii); + printf("[31] datar[j] = %f\n", valuerj); + printf("[32] datai[j] = %f\n", valueij); sumr = (valueri + valuerj) / 2.; sumi = (valueii + valueij) / 2.; @@ -575,7 +575,6 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datai, i, difi + tempi); chunk_array_save(datar, j, sumr - tempr); chunk_array_save(datai, j, tempi - difi); - //printf("563\n"); j += np2; } imin++; @@ -593,7 +592,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int goto L740; for (i = imin; i <= ntot; i += np2) { chunk_array_get(datai, i, &valueii); - printf("datai[i] = %f\n", valueii); + printf("[33] datai[i] = %f\n", valueii); chunk_array_save(datai, i, -valueii); } L740: @@ -609,8 +608,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("datar[i] = %f\n", valueri); - printf("datai[i] = %f\n", valueii); + printf("[34] datar[i] = %f\n", valueri); + printf("[35] datai[i] = %f\n", valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); @@ -622,8 +621,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("datar[imin] = %f\n", valuerimin); - printf("datai[imin] = %f\n", valueiimin); + printf("[36] datar[imin] = %f\n", valuerimin); + printf("[37] datai[imin] = %f\n", valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); @@ -636,8 +635,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("datar[i] = %f\n", valueri); - printf("datai[i] = %f\n", valueii); + printf("[38] datar[i] = %f\n", valueri); + printf("[39] datai[i] = %f\n", valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); @@ -650,8 +649,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("datar[imin] = %f\n", valuerimin); - printf("datai[imin] = %f\n", valueiimin); + printf("[40] datar[imin] = %f\n", valuerimin); + printf("[41] datai[imin] = %f\n", valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); @@ -662,8 +661,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, 1, &valuei1); chunk_array_get(datar, 1, &valuer1); - printf("datai[1] = %f\n", valuei1); - printf("datar[1] = %f\n", valuer1); + printf("[42] datai[1] = %f\n", valuei1); + printf("[43] datar[1] = %f\n", valuer1); chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); @@ -687,8 +686,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[44] datar[j] = %f\n", valuerj); + printf("[45] datai[j] = %f\n", valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); @@ -700,8 +699,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("datar[j] = %f\n", valuerj); - printf("datai[j] = %f\n", valueij); + printf("[46] datar[j] = %f\n", valuerj); + printf("[47] datai[j] = %f\n", valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); From 30a14002aea76294728ff26c5aca3238e70bd97a Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 15:37:29 -0300 Subject: [PATCH 07/15] Add indexes to print --- fftma_module/gen/lib_src/fourt.c | 90 ++++++++++++++++---------------- 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index c340173..d75bf1c 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -190,9 +190,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i = 1; for (j = 1; j <= ntot; j++) { chunk_array_get(datar, i, &valueri); - printf("[1] datar[i] = %f\n", valueri); + printf("[1] datar[%d] = %f\n", i, valueri); chunk_array_get(datar, i+1, &valueri1); - printf("[2] datar[i1] = %f\n", valueri1); + printf("[2] datar[%d] = %f\n", i+1, valueri1); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); i += 2; @@ -217,10 +217,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); - printf("[3] datar[i3] = %f\n", valueri3); - printf("[4] datai[i3] = %f\n", valueii3); - printf("[5] datar[j3] = %f\n", valuerj3); - printf("[6] datai[j3] = %f\n", valueij3); + printf("[3] datar[%d] = %f\n", i3, valueri3); + printf("[4] datai[%d] = %f\n", i3, valueii3); + printf("[5] datar[%d] = %f\n", j3, valuerj3); + printf("[6] datai[%d] = %f\n", j3, valueij3); tempr = valueri3; tempi = valueii3; @@ -256,8 +256,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[7] datar[j] = %f\n", valuerj); - printf("[8] datai[j] = %f\n", valueij); + printf("[7] datar[%d] = %f\n", j, valuerj); + printf("[8] datai[%d] = %f\n", j, valueij); if (icase != 3) { workr[i] = valuerj; @@ -309,10 +309,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[9] datar[i] = %f\n", tempr); - printf("[10] datai[i] = %f\n", tempi); - printf("[11] datar[j] = %f\n", valuerj); - printf("[12] datai[j] = %f\n", valueij); + printf("[9] datar[%d] = %f\n", i, tempr); + printf("[10] datai[%d] = %f\n", i, tempi); + printf("[11] datar[%d] = %f\n", j, valuerj); + printf("[12] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); @@ -339,22 +339,22 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); - printf("[13] datai[i] = %f\n", tempr); - printf("[14] datar[i] = %f\n", tempi); + printf("[13] datai[%d] = %f\n", i, tempr); + printf("[14] datar[%d] = %f\n", i, tempi); tempi = -tempi; } else { chunk_array_get(datai, i, &tempr); - printf("[15] datai[i] = %f\n", tempr); + printf("[15] datai[%d] = %f\n", i, tempr); tempr = -tempr; chunk_array_get(datar, i, &tempi); - printf("[16] datar[i] = %f\n", tempi); + printf("[16] datar[%d] = %f\n", i, tempi); } chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[17] datar[j] = %f\n", valuerj); - printf("[18] datai[j] = %f\n", valueij); + printf("[17] datar[%d] = %f\n", j, valuerj); + printf("[18] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); @@ -362,8 +362,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[19] datar[j] = %f\n", valuerj); - printf("[20] datai[j] = %f\n", valueij); + printf("[19] datar[%d] = %f\n", j, valuerj); + printf("[20] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); @@ -408,10 +408,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &valueii); chunk_array_get(datai, j, &valueij); - printf("[21] datar[i] = %f\n", valueri); - printf("[22] datar[j] = %f\n", valuerj); - printf("[23] datai[i] = %f\n", valueii); - printf("[24] datai[j] = %f\n", valueij); + printf("[21] datar[%d] = %f\n", i, valueri); + printf("[22] datar[%d] = %f\n", j, valuerj); + printf("[23] datai[%d] = %f\n", i, valueii); + printf("[24] datai[%d] = %f\n", j, valueij); tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; @@ -474,8 +474,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); - printf("[25] datar[j] = %f\n", sr); - printf("[26] datai[j] = %f\n", si); + printf("[25] datar[%d] = %f\n", j, sr); + printf("[26] datai[%d] = %f\n", j, si); oldsr = 0.; oldsi = 0.; @@ -487,8 +487,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[27] datar[j] = %f\n", valuerj); - printf("[28] datai[j] = %f\n", valueij); + printf("[27] datar[%d] = %f\n", j, valuerj); + printf("[28] datai[%d] = %f\n", j, valueij); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; @@ -560,10 +560,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[29] datar[i] = %f\n", valueri); - printf("[30] datai[i] = %f\n", valueii); - printf("[31] datar[j] = %f\n", valuerj); - printf("[32] datai[j] = %f\n", valueij); + printf("[29] datar[%d] = %f\n", i, valueri); + printf("[30] datai[%d] = %f\n", i, valueii); + printf("[31] datar[%d] = %f\n", j, valuerj); + printf("[32] datai[%d] = %f\n", j, valueij); sumr = (valueri + valuerj) / 2.; sumi = (valueii + valueij) / 2.; @@ -592,7 +592,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int goto L740; for (i = imin; i <= ntot; i += np2) { chunk_array_get(datai, i, &valueii); - printf("[33] datai[i] = %f\n", valueii); + printf("[33] datai[%d] = %f\n", i, valueii); chunk_array_save(datai, i, -valueii); } L740: @@ -608,8 +608,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("[34] datar[i] = %f\n", valueri); - printf("[35] datai[i] = %f\n", valueii); + printf("[34] datar[%d] = %f\n", i, valueri); + printf("[35] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); @@ -621,8 +621,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("[36] datar[imin] = %f\n", valuerimin); - printf("[37] datai[imin] = %f\n", valueiimin); + printf("[36] datar[%d] = %f\n", imin, valuerimin); + printf("[37] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); @@ -635,8 +635,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("[38] datar[i] = %f\n", valueri); - printf("[39] datai[i] = %f\n", valueii); + printf("[38] datar[%d] = %f\n", i, valueri); + printf("[39] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); @@ -649,8 +649,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("[40] datar[imin] = %f\n", valuerimin); - printf("[41] datai[imin] = %f\n", valueiimin); + printf("[40] datar[%d] = %f\n", imin, valuerimin); + printf("[41] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); @@ -686,8 +686,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[44] datar[j] = %f\n", valuerj); - printf("[45] datai[j] = %f\n", valueij); + printf("[44] datar[%d] = %f\n", j, valuerj); + printf("[45] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); @@ -699,8 +699,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[46] datar[j] = %f\n", valuerj); - printf("[47] datai[j] = %f\n", valueij); + printf("[46] datar[%d] = %f\n", j, valuerj); + printf("[47] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); From 9d92fa5be14bb3af5a0afd009c30a7051c8908d8 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 15:49:05 -0300 Subject: [PATCH 08/15] Add corrida --- fftma_module/gen/lib_src/fftma2.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index dcca2ff..207deb3 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -66,6 +66,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r covariance(covar, variogram, grid, n, cores); /*power spectrum*/ + printf("Running with covar and ireal\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); /*organization of the input Gaussian white noise*/ @@ -73,6 +74,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r prebuild_gwn(grid, n, realin, realization, solver, cores, seed); /*forward fourier transform of the GWN*/ + printf("Running with realization and ireal\n"); fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); /* build realization in spectral domain */ @@ -82,6 +84,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r remove("covar.txt"); /*backward fourier transform*/ + printf("Running with realization and ireal\n"); fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); chunk_array_free(ireal); From 0a7f2fe1ba150883445ecb3f1b63ab6d7f71d914 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 17:14:13 -0300 Subject: [PATCH 09/15] Add prints for saving as well --- fftma_module/gen/lib_src/fourt.c | 67 ++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index d75bf1c..25faaf9 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -193,6 +193,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int printf("[1] datar[%d] = %f\n", i, valueri); chunk_array_get(datar, i+1, &valueri1); printf("[2] datar[%d] = %f\n", i+1, valueri1); + printf("[48] Saving in datar the value %f in pos %d\n", valueri, j); + printf("[49] Saving in datai the value %f in pos %d\n", valueri1, j); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); i += 2; @@ -225,6 +227,11 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int tempr = valueri3; tempi = valueii3; + printf("[50] Saving in datar the value %f in pos %d\n", valuerj3, i3); + printf("[51] Saving in datai the value %f in pos %d\n", valueij3, i3); + printf("[52] Saving in datar the value %f in pos %d\n", tempr, j3); + printf("[53] Saving in datai the value %f in pos %d\n", tempi, j3); + chunk_array_save(datar, i3, valuerj3); chunk_array_save(datai, i3, valueij3); chunk_array_save(datar, j3, tempr); @@ -282,6 +289,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { + printf("[54] Saving in datar the value %f in pos %d\n", workr[i], i2); + printf("[55] Saving in datai the value %f in pos %d\n", worki[i], i2); + chunk_array_save(datar, i2, workr[i]); chunk_array_save(datai, i2, worki[i]); i++; @@ -318,6 +328,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); + + printf("[56] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + printf("[57] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + printf("[58] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + printf("[59] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + j += istep; } imin = 2 * imin - i1; @@ -359,6 +375,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); + printf("[60] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + printf("[61] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); @@ -367,6 +386,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); + + printf("[62] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + printf("[63] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + j += istep; } @@ -415,10 +438,17 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; + chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); + + printf("[64] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + printf("[65] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + printf("[66] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + printf("[67] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + j += istep; } imin = 2 * imin - i1; @@ -510,6 +540,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { + printf("[68] Saving in datar the value %f in pos %d\n", workr[i], j3); + printf("[69] Saving in datai the value %f in pos %d\n", worki[i], j3); + chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); i++; @@ -571,10 +604,17 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int difi = (valueii - valueij) / 2.; tempr = wr * sumi + wi * difr; tempi = wi * sumi - wr * difr; + chunk_array_save(datar, i, sumr + tempr); chunk_array_save(datai, i, difi + tempi); chunk_array_save(datar, j, sumr - tempr); chunk_array_save(datai, j, tempi - difi); + + printf("[70] Saving in datar the value %f in pos %d\n", sumr + tempr, i); + printf("[71] Saving in datai the value %f in pos %d\n", difi + tempi, i); + printf("[72] Saving in datar the value %f in pos %d\n", sumr - tempr, j); + printf("[73] Saving in datai the value %f in pos %d\n", tempi - difi, j); + j += np2; } imin++; @@ -594,6 +634,7 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &valueii); printf("[33] datai[%d] = %f\n", i, valueii); chunk_array_save(datai, i, -valueii); + printf("[73] Saving in datai the value %f in pos %d\n", -valueii, i); } L740: np2 *= 2; @@ -613,6 +654,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); + + printf("[74] Saving in datar the value %f in pos %d\n", valueri, j); + printf("[75] Saving in datai the value %f in pos %d\n", -valueii, j); + L755: i++; j--; @@ -626,6 +671,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); + + printf("[75] Saving in datar the value %f in pos %d\n", valuerimin - valueiimin, j); + printf("[76] Saving in datai the value %f in pos %d\n", 0., j); + if (i >= j) { goto L780; } else { @@ -640,6 +689,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); + + printf("[77] Saving in datar the value %f in pos %d\n", valueri, j); + printf("[78] Saving in datai the value %f in pos %d\n", valueii, j); + L770: i--; j--; @@ -655,6 +708,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); + printf("[75] Saving in datar the value %f in pos %d\n", valuerimin + valueiimin, j); + printf("[76] Saving in datai the value %f in pos %d\n", 0., j); + imax = imin; goto L745; L780: ; @@ -666,6 +722,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); + + printf("[77] Saving in datar the value %f in pos %d\n", valuei1 + valuer1, 1); + printf("[78] Saving in datai the value %f in pos %d\n", 0., 1); + goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ @@ -691,6 +751,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); + + printf("[79] Saving in datar the value %f in pos %d\n", valuerj, i); + printf("[80] Saving in datai the value %f in pos %d\n", -valueij, i); + j--; } } @@ -705,6 +769,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); + printf("[81] Saving in datar the value %f in pos %d\n", valuerj, i); + printf("[82] Saving in datai the value %f in pos %d\n", -valueij, i); + j -= np0; } } From 90f76bb61d7a337a643113dec2c56b703c30c974 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 17:43:25 -0300 Subject: [PATCH 10/15] Add more prints.. --- fftma_module/gen/lib_src/fourt.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index 25faaf9..b0757cb 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -527,8 +527,13 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int j -= ifp1; if (j > jmin) goto L620; + workr[i] = wr * sr - wi * si - oldsr + valuerj; worki[i] = wi * sr + wr * si - oldsi + valueij; + + printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i); + printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i); + jmin += ifp2; i++; } @@ -540,8 +545,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { - printf("[68] Saving in datar the value %f in pos %d\n", workr[i], j3); - printf("[69] Saving in datai the value %f in pos %d\n", worki[i], j3); + printf("[68] Saving in datar the value %f in pos %d, i = %d\n", workr[i], j3, i); + printf("[69] Saving in datai the value %f in pos %d, i = %d\n", worki[i], j3, i); chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); From 44a16ce9b97803c10beb2f1d9f28cb9645937549 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 18:22:53 -0300 Subject: [PATCH 11/15] Add more prints.. --- fftma_module/gen/lib_src/fourt.c | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index b0757cb..188a358 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -319,8 +319,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[9] datar[%d] = %f\n", i, tempr); - printf("[10] datai[%d] = %f\n", i, tempi); + printf("[9] tempr = %f\n", i, tempr); + printf("[10] tempi = %f\n", i, tempi); printf("[11] datar[%d] = %f\n", j, valuerj); printf("[12] datai[%d] = %f\n", j, valueij); @@ -525,12 +525,18 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int oldsr = stmpr; oldsi = stmpi; j -= ifp1; - if (j > jmin) + if (j > jmin) goto L620; + + chunk_array_get(datar, j, &valuerj); + chunk_array_get(datai, j, &valueij); workr[i] = wr * sr - wi * si - oldsr + valuerj; worki[i] = wi * sr + wr * si - oldsi + valueij; + printf("[85] wr = %f, sr = %f, wi = %f, si = %f, oldsr = %f, datar[j] = %f\n", wr, sr, wi, si, oldsr, valuerj); + printf("[86] wi = %f, sr = %f, wr = %f, si = %f, oldsi = %f, datai[j] = %f\n", wi, sr, wr, si, oldsi, valueij); + printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i); printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i); From e626ff6c930d469756323f36b19a40588ee6c921 Mon Sep 17 00:00:00 2001 From: Oli Date: Sun, 16 Jan 2022 19:05:42 -0300 Subject: [PATCH 12/15] covar test --- fftma_module/gen/lib_src/fftma2.c | 2 +- fftma_module/gen/lib_src/fourt.c | 184 +++---- fftma_module/gen/lib_src/fourt_covar.c | 645 +++++++++++++++++++++++++ 3 files changed, 738 insertions(+), 93 deletions(-) create mode 100644 fftma_module/gen/lib_src/fourt_covar.c diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 207deb3..9ead660 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -81,7 +81,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r build_real(n, NTOT, covar, realization, ireal, cores); chunk_array_free(covar); - remove("covar.txt"); + //remove("covar.txt"); /*backward fourier transform*/ printf("Running with realization and ireal\n"); diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index 188a358..caf36b6 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -190,11 +190,11 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i = 1; for (j = 1; j <= ntot; j++) { chunk_array_get(datar, i, &valueri); - printf("[1] datar[%d] = %f\n", i, valueri); + ////printf("[1] datar[%d] = %f\n", i, valueri); chunk_array_get(datar, i+1, &valueri1); - printf("[2] datar[%d] = %f\n", i+1, valueri1); - printf("[48] Saving in datar the value %f in pos %d\n", valueri, j); - printf("[49] Saving in datai the value %f in pos %d\n", valueri1, j); + ////printf("[2] datar[%d] = %f\n", i+1, valueri1); + ////printf("[48] Saving in datar the value %f in pos %d\n", valueri, j); + ////printf("[49] Saving in datai the value %f in pos %d\n", valueri1, j); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); i += 2; @@ -219,18 +219,18 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); - printf("[3] datar[%d] = %f\n", i3, valueri3); - printf("[4] datai[%d] = %f\n", i3, valueii3); - printf("[5] datar[%d] = %f\n", j3, valuerj3); - printf("[6] datai[%d] = %f\n", j3, valueij3); + //printf("[3] datar[%d] = %f\n", i3, valueri3); + //printf("[4] datai[%d] = %f\n", i3, valueii3); + //printf("[5] datar[%d] = %f\n", j3, valuerj3); + //printf("[6] datai[%d] = %f\n", j3, valueij3); tempr = valueri3; tempi = valueii3; - printf("[50] Saving in datar the value %f in pos %d\n", valuerj3, i3); - printf("[51] Saving in datai the value %f in pos %d\n", valueij3, i3); - printf("[52] Saving in datar the value %f in pos %d\n", tempr, j3); - printf("[53] Saving in datai the value %f in pos %d\n", tempi, j3); + //printf("[50] Saving in datar the value %f in pos %d\n", valuerj3, i3); + //printf("[51] Saving in datai the value %f in pos %d\n", valueij3, i3); + //printf("[52] Saving in datar the value %f in pos %d\n", tempr, j3); + //printf("[53] Saving in datai the value %f in pos %d\n", tempi, j3); chunk_array_save(datar, i3, valuerj3); chunk_array_save(datai, i3, valueij3); @@ -263,8 +263,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[7] datar[%d] = %f\n", j, valuerj); - printf("[8] datai[%d] = %f\n", j, valueij); + //printf("[7] datar[%d] = %f\n", j, valuerj); + //printf("[8] datai[%d] = %f\n", j, valueij); if (icase != 3) { workr[i] = valuerj; @@ -289,8 +289,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { - printf("[54] Saving in datar the value %f in pos %d\n", workr[i], i2); - printf("[55] Saving in datai the value %f in pos %d\n", worki[i], i2); + //printf("[54] Saving in datar the value %f in pos %d\n", workr[i], i2); + //printf("[55] Saving in datai the value %f in pos %d\n", worki[i], i2); chunk_array_save(datar, i2, workr[i]); chunk_array_save(datai, i2, worki[i]); @@ -319,20 +319,20 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[9] tempr = %f\n", i, tempr); - printf("[10] tempi = %f\n", i, tempi); - printf("[11] datar[%d] = %f\n", j, valuerj); - printf("[12] datai[%d] = %f\n", j, valueij); + //printf("[9] tempr = %f\n", i, tempr); + //printf("[10] tempi = %f\n", i, tempi); + //printf("[11] datar[%d] = %f\n", j, valuerj); + //printf("[12] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - printf("[56] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); - printf("[57] Saving in datai the value %f in pos %d\n", valueij - tempi, i); - printf("[58] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); - printf("[59] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + //printf("[56] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + //printf("[57] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + //printf("[58] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + //printf("[59] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } @@ -355,40 +355,40 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); - printf("[13] datai[%d] = %f\n", i, tempr); - printf("[14] datar[%d] = %f\n", i, tempi); + //printf("[13] datai[%d] = %f\n", i, tempr); + //printf("[14] datar[%d] = %f\n", i, tempi); tempi = -tempi; } else { chunk_array_get(datai, i, &tempr); - printf("[15] datai[%d] = %f\n", i, tempr); + //printf("[15] datai[%d] = %f\n", i, tempr); tempr = -tempr; chunk_array_get(datar, i, &tempi); - printf("[16] datar[%d] = %f\n", i, tempi); + //printf("[16] datar[%d] = %f\n", i, tempi); } chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[17] datar[%d] = %f\n", j, valuerj); - printf("[18] datai[%d] = %f\n", j, valueij); + //printf("[17] datar[%d] = %f\n", j, valuerj); + //printf("[18] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); - printf("[60] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); - printf("[61] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + //printf("[60] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + //printf("[61] Saving in datai the value %f in pos %d\n", valueij - tempi, i); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[19] datar[%d] = %f\n", j, valuerj); - printf("[20] datai[%d] = %f\n", j, valueij); + //printf("[19] datar[%d] = %f\n", j, valuerj); + //printf("[20] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - printf("[62] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); - printf("[63] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + //printf("[62] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + //printf("[63] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } @@ -431,10 +431,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, i, &valueii); chunk_array_get(datai, j, &valueij); - printf("[21] datar[%d] = %f\n", i, valueri); - printf("[22] datar[%d] = %f\n", j, valuerj); - printf("[23] datai[%d] = %f\n", i, valueii); - printf("[24] datai[%d] = %f\n", j, valueij); + //printf("[21] datar[%d] = %f\n", i, valueri); + //printf("[22] datar[%d] = %f\n", j, valuerj); + //printf("[23] datai[%d] = %f\n", i, valueii); + //printf("[24] datai[%d] = %f\n", j, valueij); tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; @@ -444,10 +444,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); - printf("[64] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); - printf("[65] Saving in datai the value %f in pos %d\n", valueij - tempi, i); - printf("[66] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); - printf("[67] Saving in datai the value %f in pos %d\n", valueij + tempi, j); + //printf("[64] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); + //printf("[65] Saving in datai the value %f in pos %d\n", valueij - tempi, i); + //printf("[66] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); + //printf("[67] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } @@ -504,8 +504,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); - printf("[25] datar[%d] = %f\n", j, sr); - printf("[26] datai[%d] = %f\n", j, si); + //printf("[25] datar[%d] = %f\n", j, sr); + //printf("[26] datai[%d] = %f\n", j, si); oldsr = 0.; oldsi = 0.; @@ -517,8 +517,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[27] datar[%d] = %f\n", j, valuerj); - printf("[28] datai[%d] = %f\n", j, valueij); + //printf("[27] datar[%d] = %f\n", j, valuerj); + //printf("[28] datai[%d] = %f\n", j, valueij); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; @@ -534,11 +534,11 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int workr[i] = wr * sr - wi * si - oldsr + valuerj; worki[i] = wi * sr + wr * si - oldsi + valueij; - printf("[85] wr = %f, sr = %f, wi = %f, si = %f, oldsr = %f, datar[j] = %f\n", wr, sr, wi, si, oldsr, valuerj); - printf("[86] wi = %f, sr = %f, wr = %f, si = %f, oldsi = %f, datai[j] = %f\n", wi, sr, wr, si, oldsi, valueij); + //printf("[85] wr = %f, sr = %f, wi = %f, si = %f, oldsr = %f, datar[j] = %f\n", wr, sr, wi, si, oldsr, valuerj); + //printf("[86] wi = %f, sr = %f, wr = %f, si = %f, oldsi = %f, datai[j] = %f\n", wi, sr, wr, si, oldsi, valueij); - printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i); - printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i); + //printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i); + //printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i); jmin += ifp2; i++; @@ -551,8 +551,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { - printf("[68] Saving in datar the value %f in pos %d, i = %d\n", workr[i], j3, i); - printf("[69] Saving in datai the value %f in pos %d, i = %d\n", worki[i], j3, i); + //printf("[68] Saving in datar the value %f in pos %d, i = %d\n", workr[i], j3, i); + //printf("[69] Saving in datai the value %f in pos %d, i = %d\n", worki[i], j3, i); chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); @@ -604,10 +604,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[29] datar[%d] = %f\n", i, valueri); - printf("[30] datai[%d] = %f\n", i, valueii); - printf("[31] datar[%d] = %f\n", j, valuerj); - printf("[32] datai[%d] = %f\n", j, valueij); + //printf("[29] datar[%d] = %f\n", i, valueri); + //printf("[30] datai[%d] = %f\n", i, valueii); + //printf("[31] datar[%d] = %f\n", j, valuerj); + //printf("[32] datai[%d] = %f\n", j, valueij); sumr = (valueri + valuerj) / 2.; sumi = (valueii + valueij) / 2.; @@ -621,10 +621,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_save(datar, j, sumr - tempr); chunk_array_save(datai, j, tempi - difi); - printf("[70] Saving in datar the value %f in pos %d\n", sumr + tempr, i); - printf("[71] Saving in datai the value %f in pos %d\n", difi + tempi, i); - printf("[72] Saving in datar the value %f in pos %d\n", sumr - tempr, j); - printf("[73] Saving in datai the value %f in pos %d\n", tempi - difi, j); + //printf("[70] Saving in datar the value %f in pos %d\n", sumr + tempr, i); + //printf("[71] Saving in datai the value %f in pos %d\n", difi + tempi, i); + //printf("[72] Saving in datar the value %f in pos %d\n", sumr - tempr, j); + //printf("[73] Saving in datai the value %f in pos %d\n", tempi - difi, j); j += np2; } @@ -643,9 +643,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int goto L740; for (i = imin; i <= ntot; i += np2) { chunk_array_get(datai, i, &valueii); - printf("[33] datai[%d] = %f\n", i, valueii); + //printf("[33] datai[%d] = %f\n", i, valueii); chunk_array_save(datai, i, -valueii); - printf("[73] Saving in datai the value %f in pos %d\n", -valueii, i); + //printf("[73] Saving in datai the value %f in pos %d\n", -valueii, i); } L740: np2 *= 2; @@ -660,14 +660,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("[34] datar[%d] = %f\n", i, valueri); - printf("[35] datai[%d] = %f\n", i, valueii); + //printf("[34] datar[%d] = %f\n", i, valueri); + //printf("[35] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); - printf("[74] Saving in datar the value %f in pos %d\n", valueri, j); - printf("[75] Saving in datai the value %f in pos %d\n", -valueii, j); + //printf("[74] Saving in datar the value %f in pos %d\n", valueri, j); + //printf("[75] Saving in datai the value %f in pos %d\n", -valueii, j); L755: i++; @@ -677,14 +677,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("[36] datar[%d] = %f\n", imin, valuerimin); - printf("[37] datai[%d] = %f\n", imin, valueiimin); + //printf("[36] datar[%d] = %f\n", imin, valuerimin); + //printf("[37] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); - printf("[75] Saving in datar the value %f in pos %d\n", valuerimin - valueiimin, j); - printf("[76] Saving in datai the value %f in pos %d\n", 0., j); + //printf("[75] Saving in datar the value %f in pos %d\n", valuerimin - valueiimin, j); + //printf("[76] Saving in datai the value %f in pos %d\n", 0., j); if (i >= j) { goto L780; @@ -695,14 +695,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); - printf("[38] datar[%d] = %f\n", i, valueri); - printf("[39] datai[%d] = %f\n", i, valueii); + //printf("[38] datar[%d] = %f\n", i, valueri); + //printf("[39] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); - printf("[77] Saving in datar the value %f in pos %d\n", valueri, j); - printf("[78] Saving in datai the value %f in pos %d\n", valueii, j); + //printf("[77] Saving in datar the value %f in pos %d\n", valueri, j); + //printf("[78] Saving in datai the value %f in pos %d\n", valueii, j); L770: i--; @@ -713,14 +713,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); - printf("[40] datar[%d] = %f\n", imin, valuerimin); - printf("[41] datai[%d] = %f\n", imin, valueiimin); + //printf("[40] datar[%d] = %f\n", imin, valuerimin); + //printf("[41] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); - printf("[75] Saving in datar the value %f in pos %d\n", valuerimin + valueiimin, j); - printf("[76] Saving in datai the value %f in pos %d\n", 0., j); + //printf("[75] Saving in datar the value %f in pos %d\n", valuerimin + valueiimin, j); + //printf("[76] Saving in datai the value %f in pos %d\n", 0., j); imax = imin; goto L745; @@ -728,14 +728,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datai, 1, &valuei1); chunk_array_get(datar, 1, &valuer1); - printf("[42] datai[1] = %f\n", valuei1); - printf("[43] datar[1] = %f\n", valuer1); + //printf("[42] datai[1] = %f\n", valuei1); + //printf("[43] datar[1] = %f\n", valuer1); chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); - printf("[77] Saving in datar the value %f in pos %d\n", valuei1 + valuer1, 1); - printf("[78] Saving in datai the value %f in pos %d\n", 0., 1); + //printf("[77] Saving in datar the value %f in pos %d\n", valuei1 + valuer1, 1); + //printf("[78] Saving in datai the value %f in pos %d\n", 0., 1); goto L900; @@ -757,14 +757,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[44] datar[%d] = %f\n", j, valuerj); - printf("[45] datai[%d] = %f\n", j, valueij); + //printf("[44] datar[%d] = %f\n", j, valuerj); + //printf("[45] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); - printf("[79] Saving in datar the value %f in pos %d\n", valuerj, i); - printf("[80] Saving in datai the value %f in pos %d\n", -valueij, i); + //printf("[79] Saving in datar the value %f in pos %d\n", valuerj, i); + //printf("[80] Saving in datai the value %f in pos %d\n", -valueij, i); j--; } @@ -774,14 +774,14 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); - printf("[46] datar[%d] = %f\n", j, valuerj); - printf("[47] datai[%d] = %f\n", j, valueij); + //printf("[46] datar[%d] = %f\n", j, valuerj); + //printf("[47] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); - printf("[81] Saving in datar the value %f in pos %d\n", valuerj, i); - printf("[82] Saving in datai the value %f in pos %d\n", -valueij, i); + //printf("[81] Saving in datar the value %f in pos %d\n", valuerj, i); + //printf("[82] Saving in datai the value %f in pos %d\n", -valueij, i); j -= np0; } diff --git a/fftma_module/gen/lib_src/fourt_covar.c b/fftma_module/gen/lib_src/fourt_covar.c new file mode 100644 index 0000000..867c511 --- /dev/null +++ b/fftma_module/gen/lib_src/fourt_covar.c @@ -0,0 +1,645 @@ +#include +#include +#include +#include "chunk_array.h" + +/*fast fourier transform */ +/* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ +/* EVALUATES COMPLEX FOURIER SERIES FOR COMPLEX OR REAL FUNCTIONS. */ +/* THAT IS, IT COMPUTES */ +/* FTRAN(J1,J2,...)=SUM(DATA(I1,I2,...)*W1**(I1-1)*(J1-1) */ +/* *W2**(I2-1)*(J2-1)*...), */ +/* WHERE W1=EXP(-2*PI*SQRT(-1)/NN(1)), W2=EXP(-2*PI*SQRT(-1)/NN(2)), */ +/* ETC. AND I1 AND J1 RUN FROM 1 TO NN(1), I2 AND J2 RUN FROM 1 TO */ +/* NN(2), ETC. THERE IS NO LIMIT ON THE DIMENSIONALITY (NUMBER OF */ +/* SUBSCRIPTS) OF THE ARRAY OF DATA. THE PROGRAM WILL PERFORM */ +/* A THREE-DIMENSIONAL FOURIER TRANSFORM AS EASILY AS A ONE-DIMEN- */ +/* SIONAL ONE, THO IN A PROPORTIONATELY GREATER TIME. AN INVERSE */ +/* TRANSFORM CAN BE PERFORMED, IN WHICH THE SIGN IN THE EXPONENTIALS */ +/* IS +, INSTEAD OF -. IF AN INVERSE TRANSFORM IS PERFORMED UPON */ +/* AN ARRAY OF TRANSFORMED DATA, THE ORIGINAL DATA WILL REAPPEAR, */ +/* MULTIPLIED BY NN(1)*NN(2)*... THE ARRAY OF INPUT DATA MAY BE */ +/* REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF */ +/* ABOUT THIRTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX. */ +/* (FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.) */ +/* THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE */ +/* ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA. THE LENGTH */ +/* OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER. THE */ +/* PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS */ +/* PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO. */ +/* TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA. LET NTOT BE THE */ +/* TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT */ +/* IS, NTOT=NN(1)*NN(2)*... DECOMPOSE NTOT INTO ITS PRIME FACTORS, */ +/* SUCH AS 2**K2 * 3**K3 * 5**K5 * ... LET SUM2 BE THE SUM OF ALL */ +/* THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2. LET SUMF BE */ +/* THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+.. */ +/* THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA */ +/* IS T = T0 + T1*NTOT + T2*NTOT*SUM2 + T3*NTOT*SUMF. FOR THE PAR- */ +/* TICULAR IMPLEMENTATION FORTRAN 32 ON THE CDC 3300 (FLOATING POINT */ +/* ADD TIME = SIX MICROSECONDS), */ +/* T = 3000 + 600*NTOT + 50*NTOT*SUM2 + 175*NTOT*SUMF MICROSECONDS */ +/* ON COMPLEX DATA. */ +/* IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME */ +/* PROPORTIONAL TO NTOT**2. FOR HIGHLY COMPOSITE NTOT, THE SAVINGS */ +/* OFFERED BY COOLEY-TUKEY CAN BE DRAMATIC. A MATRIX 100 BY 100 WILL */ +/* BE TRANSFORMED IN TIME PROPORTIONAL TO 10000*(2+2+2+2+5+5+5+5) = */ +/* 280,000 (ASSUMING T2 AND T3 TO BE ROUGHLY COMPARABLE) VERSUS */ +/* 10000**2 = 100,000,000 FOR THE STRAIGHTFORWARD TECHNIQUE. */ +/* THE COOLEY-TUKEY ALGORITHM PLACES TWO RESTRICTIONS UPON THE */ +/* NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT */ +/* THE DATA FROM ONE CYCLE OF A PERIODIC FUNCTION. THEY ARE-- */ +/* 1. THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES */ +/* MUST BE THE SAME. */ +/* 2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN, */ +/* THEY MUST BE EQUI-SPACED AT INTERVALS OF DT. FURTHER, THE TRANS- */ +/* FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI- */ +/* SPACED FROM 0 TO 2*PI*(NN(I)-1)/(NN(I)*DT) AT INTERVALS OF */ +/* 2*PI/(NN(I)*DT) FOR EACH DIMENSION OF LENGTH NN(I). OF COURSE, */ +/* DT NEED NOT BE THE SAME FOR EVERY DIMENSION. */ + +/* THE CALLING SEQUENCE IS-- */ +/* CALL FOURT(DATAR,DATAI,NN,NDIM,IFRWD,ICPLX,WORKR,WORKI) */ + +/* DATAR AND DATAI ARE THE ARRAYS USED TO HOLD THE REAL AND IMAGINARY */ +/* PARTS OF THE INPUT DATA ON INPUT AND THE TRANSFORM VALUES ON */ +/* OUTPUT. THEY ARE FLOATING POINT ARRAYS, MULTIDIMENSIONAL WITH */ +/* IDENTICAL DIMENSIONALITY AND EXTENT. THE EXTENT OF EACH DIMENSION */ +/* IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM. THAT IS, */ +/* NDIM IS THE DIMENSIONALITY OF THE ARRAYS DATAR AND DATAI. */ +/* IFRWD IS AN INTEGER USED TO INDICATE THE DIRECTION OF THE FOURIER */ +/* TRANSFORM. IT IS NON-ZERO TO INDICATE A FORWARD TRANSFORM */ +/* (EXPONENTIAL SIGN IS -) AND ZERO TO INDICATE AN INVERSE TRANSFORM */ +/* (SIGN IS +). ICPLX IS AN INTEGER TO INDICATE WHETHER THE DATA */ +/* ARE REAL OR COMPLEX. IT IS NON-ZERO FOR COMPLEX, ZERO FOR REAL. */ +/* IF IT IS ZERO (REAL) THE CONTENTS OF ARRAY DATAI WILL BE ASSUMED */ +/* TO BE ZERO, AND NEED NOT BE EXPLICITLY SET TO ZERO. AS EXPLAINED */ +/* ABOVE, THE TRANSFORM RESULTS ARE ALWAYS COMPLEX AND ARE STORED */ +/* IN DATAR AND DATAI ON RETURN. WORKR AND WORKI ARE ARRAYS USED */ +/* FOR WORKING STORAGE. THEY ARE NOT NECESSARY IF ALL THE DIMENSIONS */ +/* OF THE DATA ARE POWERS OF TWO. IN THIS CASE, THE ARRAYS MAY BE */ +/* REPLACED BY THE NUMBER 0 IN THE CALLING SEQUENCE. THUS, USE OF */ +/* POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE. IF ANY DIMENSION */ +/* IS NOT A POWER OF TWO, THESE ARRAYS MUST BE SUPPLIED. THEY ARE */ +/* FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO THE LARGEST */ +/* ARRAY DIMENSION, THAT IS, TO THE LARGEST VALUE OF NN(I). */ +/* WORKR AND WORKI, IF SUPPLIED, MUST NOT BE THE SAME ARRAYS AS DATAR */ +/* OR DATAI. ALL SUBSCRIPTS OF ALL ARRAYS BEGIN AT 1. */ + +/* THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM. THE */ +/* PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE. */ + +/* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ +/* 10-06-2000, MLR */ + +void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { + int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; + double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; + double value1, valuei, valuej, valuei1, valueimin, valuei3, valuej3; + + ntot = 1; + for (idim = 0; idim < ndim; idim++) { + ntot *= nn[idim]; + } + + chunk_array_read(datar); + + /*main loop for each dimension*/ + np1 = 1; + for (idim = 1; idim <= ndim; idim++) { + n = nn[idim - 1]; + np2 = np1 * n; + + if (n < 1) { + goto L920; + } else if (n == 1) { + goto L900; + } + + /*is n a power of 2 and if not, what are its factors*/ + m = n; + ntwo = np1; + iff = 1; + idiv = 2; + + L10: + iquot = m / idiv; + irem = m - idiv * iquot; + if (iquot < idiv) + goto L50; + if (irem == 0) { + ntwo *= 2; + ifact[iff] = idiv; + iff++; + m = iquot; + goto L10; + } + idiv = 3; + inon2 = iff; + + L30: + iquot = m / idiv; + irem = m - idiv * iquot; + if (iquot < idiv) + goto L60; + if (irem == 0) { + ifact[iff] = idiv; + iff++; + m = iquot; + goto L30; + } + + idiv += 2; + goto L30; + + L50: + inon2 = iff; + if (irem != 0) + goto L60; + ntwo *= 2; + goto L70; + + L60: + ifact[iff] = m; + + L70: + non2p = np2 / ntwo; + + /*SEPARATE FOUR CASES-- + 1. COMPLEX TRANSFORM + 2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD: TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CONJUGATE SYMMETRY. + 3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD: SET THE IMAGINARY PARTS TO ZERO. + 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD: TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS ARE THE ODD-NUMBERED REAL VALUES. UNSCRAMBLE AND SUPPLY THE SECOND HALF BY CONJUGATE SYMMETRY. */ + + icase = 1; + ifmin = 1; + if (icplx != 0) + goto L100; + icase = 2; + if (idim > 1) + goto L100; + icase = 3; + if (ntwo <= np1) + goto L100; + icase = 4; + ifmin = 2; + ntwo /= 2; + n /= 2; + np2 /= 2; + ntot /= 2; + i = 1; + for (j = 1; j <= ntot; j++) { + chunk_array_get(datar, i, &valuei); + chunk_array_get(datar, i, &valuei1); + chunk_array_save(datar, j, valuei); + + //datar[j] = datar[i]; + datai[j] = valuei1; + i += 2; + } + + /*shuffle data by bit reversal, since n = 2^k. As the shuffling can be done by simple interchange, no working array is needed*/ + L100: + if (non2p > 1) + goto L200; + np2hf = np2 / 2; + j = 1; + for (i2 = 1; i2 <= np2; i2 += np1) { + if (j >= i2) + goto L130; + i1max = i2 + np1 - 1; + for (i1 = i2; i1 <= i1max; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + j3 = j + i3 - i2; + //tempr = datar[i3]; + tempi = datai[i3]; + //datar[i3] = datar[j3]; + datai[i3] = datai[j3]; + //datar[j3] = tempr; + datai[j3] = tempi; + + chunk_array_get(datar, i3, &valuei3); + chunk_array_get(datar, j3, &valuej3); + chunk_array_save(datar, i3, valuej3); + chunk_array_save(datar, j3, valuei3); + } + } + + L130: + m = np2hf; + + L140: + if (j <= m) { + j += m; + } else { + j -= m; + m /= 2; + if (m >= np1) + goto L140; + } + } + goto L300; + + /*shuffle data by digit reversal for general n*/ + L200: + for (i1 = 1; i1 <= np1; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + j = i3; + for (i = 1; i <= n; i++) { + if (icase != 3) { + //workr[i] = datar[j]; + chunk_array_get(datar, j, &workr[i]); + worki[i] = datai[j]; + } else { + chunk_array_get(datar, j, &workr[i]); + //workr[i] = datar[j]; + worki[i] = 0.; + } + ifp2 = np2; + iff = ifmin; + L250: + ifp1 = ifp2 / ifact[iff]; + j += ifp1; + if (j >= i3 + ifp2) { + j -= ifp2; + ifp2 = ifp1; + iff += 1; + if (ifp2 > np1) + goto L250; + } + } + i2max = i3 + np2 - np1; + i = 1; + for (i2 = i3; i2 <= i2max; i2 += np1) { + chunk_array_save(datar, i2, workr[i]); + //datar[i2] = workr[i]; + datai[i2] = worki[i]; + i++; + } + } + } + + /*special case-- W=1*/ + L300: + i1rng = np1; + if (icase == 2) + i1rng = np0 * (1 + nprev / 2); + if (ntwo <= np1) + goto L600; + for (i1 = 1; i1 <= i1rng; i1++) { + imin = np1 + i1; + istep = 2 * np1; + goto L330; + + L310: + j = i1; + for (i = imin; i <= ntot; i += istep) { + //tempr = datar[i]; + tempi = datai[i]; + //datar[i] = datar[j] - tempr; + datai[i] = datai[j] - tempi; + //datar[j] = datar[j] + tempr; + datai[j] = datai[j] + tempi; + + chunk_array_get(datar, i, &valuei); + chunk_array_get(datar, j, &valuej); + + chunk_array_save(datar, i, valuej - valuei); + chunk_array_save(datar, j, valuej + valuei); + + j += istep; + } + imin = 2 * imin - i1; + istep *= 2; + + L330: + if (istep <= ntwo) + goto L310; + + /*special case-- W = -sqrt(-1)*/ + imin = 3 * np1 + i1; + istep = 4 * np1; + goto L420; + + L400: + j = imin - istep / 2; + for (i = imin; i <= ntot; i += istep) { + if (ifrwd != 0) { + tempr = datai[i]; + //tempi = -datar[i]; + chunk_array_get(datar, i, &tempi); + tempi = -tempi; + } else { + tempr = -datai[i]; + //tempi = datar[i]; + chunk_array_get(datar, i, &tempi); + } + + chunk_array_get(datar, j, &valuej); + chunk_array_save(datar, i, valuej - tempr); + chunk_array_save(datar, j, valuej - tempr); + + //datar[i] = datar[j] - tempr; + datai[i] = datai[j] - tempi; + //datar[j] += tempr; + datai[j] += tempi; + j += istep; + } + + imin = 2 * imin - i1; + istep *= 2; + + L420: + if (istep <= ntwo) + goto L400; + } + + /*main loop for factors of 2. W=EXP(-2*PI*SQRT(-1)*m/mmax) */ + theta = -TWOPI / 8.; + wstpr = 0.; + wstpi = -1.; + if (ifrwd == 0) { + theta = -theta; + wstpi = 1.; + } + mmax = 8 * np1; + goto L540; + + L500: + wminr = cos(theta); + wmini = sin(theta); + wr = wminr; + wi = wmini; + mmin = mmax / 2 + np1; + mstep = np1 * 2; + for (m = mmin; m <= mmax; m += mstep) { + for (i1 = 1; i1 <= i1rng; i1++) { + istep = mmax; + imin = m + i1; + L510: + j = imin - istep / 2; + for (i = imin; i <= ntot; i += istep) { + double valuei, valuej; + chunk_array_get(datar, i, &valuei); + chunk_array_get(datar, j, &valuej); + tempr = valuei * wr - datai[i] * wi; + tempi = valuei * wi + datai[i] * wr; + chunk_array_save(datar, i, valuej - tempr); + //datar[i] = valuej - tempr; + datai[i] = datai[j] - tempi; + chunk_array_save(datar, i, valuej + tempr); + //datar[j] += tempr; + datai[j] += tempi; + j += istep; + } + imin = 2 * imin - i1; + istep *= 2; + if (istep <= ntwo) + goto L510; + } + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; + } + wstpr = wminr; + wstpi = wmini; + theta /= 2.; + mmax += mmax; + L540: + if (mmax <= ntwo) + goto L500; + + /*main loop for factors not equal to 2-- W=EXP(-2*PI*SQRT(-1)*(j2-i3)/ifp2)*/ + L600: + if (non2p <= 1) + goto L700; + ifp1 = ntwo; + iff = inon2; + L610: + ifp2 = ifact[iff] * ifp1; + theta = -TWOPI / (double)ifact[iff]; + if (ifrwd == 0) + theta = -theta; + thetm = theta / (double)(ifp1 / np1); + wstpr = cos(theta); + wstpi = sin(theta); + wmstr = cos(thetm); + wmsti = sin(thetm); + wminr = 1.; + wmini = 0.; + + for (j1 = 1; j1 <= ifp1; j1 += np1) { + i1max = j1 + i1rng - 1; + for (i1 = j1; i1 <= i1max; i1++) { + for (i3 = i1; i3 <= ntot; i3 += np2) { + i = 1; + wr = wminr; + wi = wmini; + j2max = i3 + ifp2 - ifp1; + for (j2 = i3; j2 <= j2max; j2 += ifp1) { + twowr = 2. * wr; + jmin = i3; + j3max = j2 + np2 - ifp2; + for (j3 = j2; j3 <= j3max; j3 += ifp2) { + j = jmin + ifp2 - ifp1; + //sr = datar[j]; + chunk_array_get(datar, j, &sr); + si = datai[j]; + oldsr = 0.; + oldsi = 0.; + j -= ifp1; + L620: + stmpr = sr; + stmpi = si; + chunk_array_get(datar, j, &valuej); + sr = twowr * sr - oldsr + valuej; + si = twowr * si - oldsi + datai[j]; + oldsr = stmpr; + oldsi = stmpi; + j -= ifp1; + if (j > jmin) + goto L620; + workr[i] = wr * sr - wi * si - oldsr + valuej; + worki[i] = wi * sr + wr * si - oldsi + datai[j]; + jmin += ifp2; + i++; + } + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; + } + i = 1; + for (j2 = i3; j2 <= j2max; j2 += ifp1) { + j3max = j2 + np2 - ifp2; + for (j3 = j2; j3 <= j3max; j3 += ifp2) { + //datar[j3] = workr[i]; + chunk_array_save(datar, j3, workr[i]); + datai[j3] = worki[i]; + i++; + } + } + } + } + wtemp = wminr * wmsti; + wminr = wminr * wmstr - wmini * wmsti; + wmini = wmini * wmstr + wtemp; + } + iff++; + ifp1 = ifp2; + if (ifp1 < np2) + goto L610; + + /*complete a real transform in the 1st dimension, n even, by conjugate symmetries*/ + L700: + switch (icase) { + case 1: + goto L900; + break; + case 2: + goto L800; + break; + case 3: + goto L900; + break; + } + + nhalf = n; + n += n; + theta = -TWOPI / (double)n; + if (ifrwd == 0) + theta = -theta; + wstpr = cos(theta); + wstpi = sin(theta); + wr = wstpr; + wi = wstpi; + imin = 2; + jmin = nhalf; + goto L725; + L710: + j = jmin; + for (i = imin; i <= ntot; i += np2) { + double valuei, valuej; + chunk_array_get(datar, i, &valuei); + chunk_array_get(datar, j, &valuej); + sumr = (valuei + valuej) / 2.; + sumi = (datai[i] + datai[j]) / 2.; + difr = (valuei - valuej) / 2.; + difi = (datai[i] - datai[j]) / 2.; + tempr = wr * sumi + wi * difr; + tempi = wi * sumi - wr * difr; + chunk_array_save(datar, i, sumr + tempr); + //datar[i] = sumr + tempr; + datai[i] = difi + tempi; + chunk_array_save(datar, j, sumr - tempr); + //datar[j] = sumr - tempr; + datai[j] = tempi - difi; + j += np2; + } + imin++; + jmin--; + wtemp = wr * wstpi; + wr = wr * wstpr - wi * wstpi; + wi = wi * wstpr + wtemp; + L725: + if (imin < jmin) { + goto L710; + } else if (imin > jmin) { + goto L740; + } + if (ifrwd == 0) + goto L740; + for (i = imin; i <= ntot; i += np2) { + datai[i] = -datai[i]; + } + L740: + np2 *= 2; + ntot *= 2; + j = ntot + 1; + imax = ntot / 2 + 1; + L745: + imin = imax - nhalf; + i = imin; + goto L755; + L750: + //datar[j] = datar[i]; + chunk_array_get(datar, i, &valuei); + chunk_array_save(datar, j, valuei); + datai[j] = -datai[i]; + L755: + i++; + j--; + if (i < imax) + goto L750; + + chunk_array_get(datar, imin, &valueimin); + chunk_array_save(datar, j, valueimin - datai[imin]); + //datar[j] = datar[imin] - datai[imin]; + datai[j] = 0.; + if (i >= j) { + goto L780; + } else { + goto L770; + } + L765: + //datar[j] = datar[i]; + chunk_array_get(datar, i, &valuei); + chunk_array_save(datar, j, valuei); + datai[j] = datai[i]; + L770: + i--; + j--; + if (i > imin) + goto L765; + //datar[j] = datar[imin] + datai[imin]; + chunk_array_get(datar, imin, &valueimin); + chunk_array_save(datar, j, valueimin - datai[imin]); + datai[j] = 0.; + imax = imin; + goto L745; + L780: + chunk_array_get(datar, 1, &value1); + chunk_array_save(datar, 1, value1 + datai[1]); + //datar[1] += datai[1]; + datai[1] = 0.; + goto L900; + + /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ + L800: + if (nprev <= 2) + goto L900; + for (i3 = 1; i3 <= ntot; i3 += np2) { + i2max = i3 + np2 - np1; + for (i2 = i3; i2 <= i2max; i2 += np1) { + imax = i2 + np1 - 1; + imin = i2 + i1rng; + jmax = 2 * i3 + np1 - imin; + if (i2 > i3) + jmax += np2; + if (idim > 2) { + j = jmax + np0; + for (i = imin; i <= imax; i++) { + //datar[i] = datar[j]; + chunk_array_get(datar, j, &valuej); + chunk_array_save(datar, i, valuej); + datai[i] = -datai[j]; + j--; + } + } + j = jmax; + for (i = imin; i <= imax; i += np0) { + //datar[i] = datar[j]; + chunk_array_get(datar, j, &valuej); + chunk_array_save(datar, i, valuej); + datai[i] = -datai[j]; + j -= np0; + } + } + } + + /*end of loop on each dimension*/ + L900: + np0 = np1; + np1 = np2; + nprev = n; + } +L920: return; +} From 897690a1e46b5abb52b4dbf829bf4dd732f19d68 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 19:44:15 -0300 Subject: [PATCH 13/15] wip --- fftma_module/gen/lib_src/fftma2.c | 27 +-------------------------- 1 file changed, 1 insertion(+), 26 deletions(-) diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 9ead660..341774b 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -69,30 +69,5 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r printf("Running with covar and ireal\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); - /*organization of the input Gaussian white noise*/ - solver = 0; - prebuild_gwn(grid, n, realin, realization, solver, cores, seed); - - /*forward fourier transform of the GWN*/ - printf("Running with realization and ireal\n"); - fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); - - /* build realization in spectral domain */ - build_real(n, NTOT, covar, realization, ireal, cores); - - chunk_array_free(covar); - //remove("covar.txt"); - - /*backward fourier transform*/ - printf("Running with realization and ireal\n"); - fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); - - chunk_array_free(ireal); - remove("ireal.txt"); - - free(workr); - free(worki); - - /*output realization*/ - clean_real(realin, n, grid, realization, realout, cores); + chunk_array_flush(covar); } From 6b4bb901e2449297a44a4d90e9fee3f607c14cfc Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 19:48:38 -0300 Subject: [PATCH 14/15] Revert "wip" This reverts commit 897690a1e46b5abb52b4dbf829bf4dd732f19d68. --- fftma_module/gen/lib_src/fftma2.c | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 341774b..9ead660 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -69,5 +69,30 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r printf("Running with covar and ireal\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); - chunk_array_flush(covar); + /*organization of the input Gaussian white noise*/ + solver = 0; + prebuild_gwn(grid, n, realin, realization, solver, cores, seed); + + /*forward fourier transform of the GWN*/ + printf("Running with realization and ireal\n"); + fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); + + /* build realization in spectral domain */ + build_real(n, NTOT, covar, realization, ireal, cores); + + chunk_array_free(covar); + //remove("covar.txt"); + + /*backward fourier transform*/ + printf("Running with realization and ireal\n"); + fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); + + chunk_array_free(ireal); + remove("ireal.txt"); + + free(workr); + free(worki); + + /*output realization*/ + clean_real(realin, n, grid, realization, realout, cores); } From 83affab77b74a22f906b60395175c6b6439c9678 Mon Sep 17 00:00:00 2001 From: chortas Date: Sun, 16 Jan 2022 19:58:45 -0300 Subject: [PATCH 15/15] Working for 16 --- fftma_module/gen/lib_src/chunk_array.c | 27 ++++++++++++++++++++++--- fftma_module/gen/lib_src/fftma2.c | 5 +---- fftma_module/gen/lib_src/prebuild_gwn.c | 7 ------- fftma_module/gen/test.py | 4 ++-- 4 files changed, 27 insertions(+), 16 deletions(-) diff --git a/fftma_module/gen/lib_src/chunk_array.c b/fftma_module/gen/lib_src/chunk_array.c index 58d0c74..607bdbf 100644 --- a/fftma_module/gen/lib_src/chunk_array.c +++ b/fftma_module/gen/lib_src/chunk_array.c @@ -7,14 +7,20 @@ void chunk_array_free(chunk_array_t* chunk_array) { free(chunk_array); } -bool chunk_array_update_read(chunk_array_t* chunk_array) { +bool chunk_array_update_read(chunk_array_t* chunk_array, size_t pos) { + + int init_pos = pos/chunk_array->chunk_size; + + fseek(chunk_array->fp, init_pos * chunk_array->chunk_size * sizeof(double), SEEK_SET); + size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); chunk_array->init_pos += newLen; } +/* bool chunk_array_get(chunk_array_t* chunk_array, size_t pos, double *valor) { if (pos>((chunk_array->init_pos + chunk_array->chunk_size)-1)) { - chunk_array_update_read(chunk_array); + chunk_array_update_read(chunk_array, pos); } *valor=chunk_array->data[pos%chunk_array->chunk_size]; return true; @@ -27,6 +33,7 @@ bool chunk_array_save(chunk_array_t* chunk_array, size_t pos, double valor) { chunk_array->data[pos%chunk_array->chunk_size]=valor; return true; } +*/ chunk_array_t* chunk_array_create(char* filename, size_t total_size, size_t chunk_size) { chunk_array_t* chunk_array = (chunk_array_t*)malloc(sizeof(chunk_array_t)); @@ -56,6 +63,7 @@ void chunk_array_read(chunk_array_t* chunk_array) { size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); } +/* void chunk_array_write(chunk_array_t* chunk_array, char* filename) { chunk_array->fp = fopen(filename, "w"); if (chunk_array->fp == NULL) { @@ -63,9 +71,22 @@ void chunk_array_write(chunk_array_t* chunk_array, char* filename) { chunk_array->fp = fopen(filename, "w"); } chunk_array->init_pos = 0; -} +}*/ void chunk_array_flush(chunk_array_t* chunk_array) { size_t newLen = fwrite(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp); chunk_array->init_pos += newLen; } + + +bool chunk_array_get(chunk_array_t* chunk_array, size_t pos, double *valor) { + fseek(chunk_array->fp, pos * sizeof(double), SEEK_SET); + fread(valor, sizeof(double), 1, chunk_array->fp); + return true; +} + +bool chunk_array_save(chunk_array_t* chunk_array, size_t pos, double valor) { + fseek(chunk_array->fp, pos * sizeof(double), SEEK_SET); + fwrite(&valor, sizeof(double), 1, chunk_array->fp); + return true; +} \ No newline at end of file diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 9ead660..dcca2ff 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -66,7 +66,6 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r covariance(covar, variogram, grid, n, cores); /*power spectrum*/ - printf("Running with covar and ireal\n"); fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); /*organization of the input Gaussian white noise*/ @@ -74,17 +73,15 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r prebuild_gwn(grid, n, realin, realization, solver, cores, seed); /*forward fourier transform of the GWN*/ - printf("Running with realization and ireal\n"); fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); /* build realization in spectral domain */ build_real(n, NTOT, covar, realization, ireal, cores); chunk_array_free(covar); - //remove("covar.txt"); + remove("covar.txt"); /*backward fourier transform*/ - printf("Running with realization and ireal\n"); fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); chunk_array_free(ireal); diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c index a6382cf..e228890 100644 --- a/fftma_module/gen/lib_src/prebuild_gwn.c +++ b/fftma_module/gen/lib_src/prebuild_gwn.c @@ -34,14 +34,10 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin ntot = n[0] * n[1] * n[2]; chunk_array_save(realization, 0, 0.); - /*printf("Antes de llamar a chunkarray read\n"); - chunk_array_read((*realin).vector_2); - printf("Despues de llamar a chunkarray read\n");*/ if (solver == 1) { for (i = 0; i < ntot; i++) { double value = gasdev(seed, &idum2, &iy, iv, cores); chunk_array_save(realization, i+1, value); - //chunk_array_get((*realin).vector_2, i, &realization[i + 1]); } } else { for (k = 1; k <= n[2]; k++) { @@ -49,11 +45,8 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin for (i = 1; i <= n[0]; i++) { maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) { - //maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; - //printf("Maille0 es %d", maille0); double value = gasdev(seed, &idum2, &iy, iv, cores); chunk_array_save(realization, maille1, value); - //chunk_array_get((*realin).vector_2, maille0, &realization[maille1]); } else { chunk_array_save(realization, maille1, 0.); } diff --git a/fftma_module/gen/test.py b/fftma_module/gen/test.py index 660a103..7b23a40 100644 --- a/fftma_module/gen/test.py +++ b/fftma_module/gen/test.py @@ -29,9 +29,9 @@ variance=3.5682389 typ=3 k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8) +k2 = np.load("out_2_16.npy") -#np.save(f"out_{N}.npy",k) -#k2 = np.load("out_2_16.npy") +print(k - k2) del k gc.collect() \ No newline at end of file