From e626ff6c930d469756323f36b19a40588ee6c921 Mon Sep 17 00:00:00 2001 From: Oli Date: Sun, 16 Jan 2022 19:05:42 -0300 Subject: [PATCH] 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; +}