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; }