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