From cceba9a83c90643b2aac2cbd5c684eaa32e2c3fb Mon Sep 17 00:00:00 2001 From: Oli Date: Thu, 20 Jan 2022 20:38:52 -0300 Subject: [PATCH] rm cores param --- fftma_module/gen/example.py | 2 +- fftma_module/gen/include/Py_py-api.h | 4 +- fftma_module/gen/include/condor.h | 2 +- fftma_module/gen/include/geostat.h | 22 +- fftma_module/gen/include/toolsFFTMA.h | 8 +- fftma_module/gen/lib_src/Py_getvalues.c | 7 +- fftma_module/gen/lib_src/Py_kgeneration.c | 6 +- fftma_module/gen/lib_src/build_real.c | 2 +- fftma_module/gen/lib_src/cardsin.c | 2 +- fftma_module/gen/lib_src/cgrid.c | 4 +- fftma_module/gen/lib_src/clean_real.c | 2 +- fftma_module/gen/lib_src/cov_value.c | 8 +- fftma_module/gen/lib_src/covariance.c | 10 +- fftma_module/gen/lib_src/cubic.c | 2 +- fftma_module/gen/lib_src/fftma2.c | 18 +- fftma_module/gen/lib_src/fourt.c | 2 +- fftma_module/gen/lib_src/fourt_covar.c | 645 ---------------------- fftma_module/gen/lib_src/gammf.c | 2 +- fftma_module/gen/lib_src/gasdev.c | 8 +- fftma_module/gen/lib_src/generate.c | 2 +- fftma_module/gen/lib_src/geostat.h | 22 +- fftma_module/gen/lib_src/length.c | 8 +- fftma_module/gen/lib_src/maxfactor.c | 2 +- fftma_module/gen/lib_src/prebuild_gwn.c | 6 +- fftma_module/gen/lib_src/ran2.c | 2 +- fftma_module/gen/moduleFFTMA.c | 5 +- tests/stages/generation/test.py | 2 +- 27 files changed, 79 insertions(+), 726 deletions(-) delete mode 100644 fftma_module/gen/lib_src/fourt_covar.c diff --git a/fftma_module/gen/example.py b/fftma_module/gen/example.py index c41a73c..057c0f5 100644 --- a/fftma_module/gen/example.py +++ b/fftma_module/gen/example.py @@ -27,5 +27,5 @@ mean=15.3245987 variance=3.5682389 typ=3 -k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8) +k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ) np.save(f"out_{N}.npy",k) \ No newline at end of file diff --git a/fftma_module/gen/include/Py_py-api.h b/fftma_module/gen/include/Py_py-api.h index 4506469..8bf84a9 100644 --- a/fftma_module/gen/include/Py_py-api.h +++ b/fftma_module/gen/include/Py_py-api.h @@ -11,5 +11,5 @@ #include #include -int Py_getvalues(PyObject*, long*, struct grid_mod*, struct vario_mod*, struct statistic_mod*, int*); -void Py_kgeneration(long, struct grid_mod, struct statistic_mod, struct vario_mod, struct realization_mod*, struct realization_mod*, int[3], int); +int Py_getvalues(PyObject*, long*, struct grid_mod*, struct vario_mod*, struct statistic_mod*); +void Py_kgeneration(long, struct grid_mod, struct statistic_mod, struct vario_mod, struct realization_mod*, struct realization_mod*, int[3]); diff --git a/fftma_module/gen/include/condor.h b/fftma_module/gen/include/condor.h index 869144b..310a594 100644 --- a/fftma_module/gen/include/condor.h +++ b/fftma_module/gen/include/condor.h @@ -102,7 +102,7 @@ void FFTMA_gradient(struct vario_mod variogram, struct grid_mod grid, int n[3], /* n: number of components in the vector */ /*output: */ /* realization: structure defining the realization*/ -void generate(long* seed, int n, struct realization_mod* realization, int cores); +void generate(long* seed, int n, struct realization_mod* realization); void ikrig(int option, struct statfacies_mod facies, int kk, struct welldata_mod data, struct variotable_mod variogram, struct grid_mod grid, float* krigout); diff --git a/fftma_module/gen/include/geostat.h b/fftma_module/gen/include/geostat.h index 9147e82..3a6c69c 100644 --- a/fftma_module/gen/include/geostat.h +++ b/fftma_module/gen/include/geostat.h @@ -299,7 +299,7 @@ struct realization_mod { void axes(double* ap, double* scf, int N); /*cardsin covariance value for lag h*/ -double cardsin(double h, int cores); +double cardsin(double h); /*Cholesky decomposition of matrix C */ /* C : symetric positive-definite matrix recorded */ @@ -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(file_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); +void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3]); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ @@ -357,10 +357,10 @@ void cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well, /*dj: distance along the Y axis */ /*dk: distance along the Z axis */ /* The returned value is the computed covariance value */ -double cov_value(struct vario_mod variogram, double di, double dj, double dk, int cores); +double cov_value(struct vario_mod variogram, double di, double dj, double dk); /*cubic covariance value for lag h*/ -double cubic(double h, int cores); +double cubic(double h); /*truncation of the power spectrum to remove */ /*high frequencies - isotropic case */ @@ -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(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); +void fourt(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki); /*calculates F(x) = (1/a)*exp(-x*x/2)*/ double funtrun1(double x); @@ -412,7 +412,7 @@ double funtrun1(double x); float G(float x); /*gamma covariance value for lag h and exponent alpha*/ -double gammf(double h, double alpha, int cores); +double gammf(double h, double alpha); /*returns the value ln(G(x))*/ float gammln(float xx); @@ -424,7 +424,7 @@ float gammp(float a, float x); /*and unit variance, using ran1(idum) as the source */ /*of uniform deviates */ /*idum: seed */ -double gasdev(long* idum, long* idum2, long* iy, long* iv, int cores); +double gasdev(long* idum, long* idum2, long* iy, long* iv); /*gaussian covariance value for lag h*/ double gaussian(double h); @@ -479,7 +479,7 @@ void gradual(struct grad_mod grad, float* Zo, float* Z, float* Zfinal, int n, st /*n: vector with the number of cells along the */ /* 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); +void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]); /*incomplete gamma function evaluated by its series*/ /*representation as gamser, also returns ln(G(a)) */ @@ -504,7 +504,7 @@ void krig_stat(float* b, int n, struct vario_mod variogram, struct welldata_mod /*i: considered direction */ /*scf: correlation length */ /*ap: normalized anisotropy axes */ -int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores); +int length(int N, int i, double* scf, double* ap, double D, int Nvari); /*calculates L.Z/ /* L : lower triangular matrix recorded */ @@ -518,7 +518,7 @@ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores void LtimeZ(double* L, float* Z, float* b, int n); /*determines the greatest prime factor of an integer*/ -int maxfactor(int n, int cores); +int maxfactor(int n); /*metrop returns a boolean varible that issues a */ /*verdict on whether to accept a reconfiguration */ @@ -544,7 +544,7 @@ double power(double h, double alpha); /*generates uniform deviates between 0 and 1*/ /*idum: seed */ -double ran2(long* idum, long* idum2, long* iy, long* iv, int cores); +double ran2(long* idum, long* idum2, long* iy, long* iv); /*calculates bt.b */ /* b : vector, bi, i = [0...n-1] */ diff --git a/fftma_module/gen/include/toolsFFTMA.h b/fftma_module/gen/include/toolsFFTMA.h index d5c4791..5d952a4 100644 --- a/fftma_module/gen/include/toolsFFTMA.h +++ b/fftma_module/gen/include/toolsFFTMA.h @@ -37,7 +37,7 @@ /*realout4: structure defining a yvelocity field */ /*realout5: structure defining a zvelocity field */ -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); +void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct realization_mod* realin, struct realization_mod* realout, long* seed); /* prebuild_gwn */ /* Produce a first construction in real space of the Gaussian white noise */ @@ -52,7 +52,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, file_array_t* realization, int solver, int cores, long* seed); +void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, file_array_t* realization, int solver, long* seed); /* build_real */ /* build a realization in the spectral domain */ @@ -66,8 +66,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, file_array_t* covar, file_array_t* realization, file_array_t* ireal, int cores); +void build_real(int n[3], int NTOT, file_array_t* covar, file_array_t* realization, file_array_t* ireal); -void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, file_array_t* vectorresult, struct realization_mod* realout, int cores); +void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, file_array_t* vectorresult, struct realization_mod* realout); #endif // define _TOOLSFFTMA_H diff --git a/fftma_module/gen/lib_src/Py_getvalues.c b/fftma_module/gen/lib_src/Py_getvalues.c index 3c1c8db..7bb857c 100644 --- a/fftma_module/gen/lib_src/Py_getvalues.c +++ b/fftma_module/gen/lib_src/Py_getvalues.c @@ -41,7 +41,7 @@ #endif #endif -int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario_mod* variogram, struct statistic_mod* stat, int* cores) { +int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario_mod* variogram, struct statistic_mod* stat) { clock_t t = clock(); int i, varioNargs = 12, j = 0; PyObject* listvario; @@ -59,7 +59,7 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario return 0; } - if (!PyArg_ParseTuple(args, "iiidddlO!ddii", /*"iiidddslO!ddi",*/ + if (!PyArg_ParseTuple(args, "iiidddlO!ddi", /*"iiidddslO!ddi",*/ &(grid->NX), &(grid->NY), &(grid->NZ), @@ -72,8 +72,7 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario &listvario, stat->mean + 0, stat->variance + 0, - &(stat->type), - cores)) { + &(stat->type))) { free(stat->mean); free(stat->variance); return 0; diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c index c633ce9..af65bc7 100644 --- a/fftma_module/gen/lib_src/Py_kgeneration.c +++ b/fftma_module/gen/lib_src/Py_kgeneration.c @@ -17,7 +17,7 @@ /* Z is the GWN with 0-mean and 1-variance */ /* 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) { +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 i, N; int typelog; @@ -27,10 +27,10 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat, n[1] = 0; n[2] = 0; - generate(&seed, N, Z, cores); + generate(&seed, N, Z); /*FFTMA*/ - FFTMA2(variogram, grid, n, Z, Y, cores, &seed); + FFTMA2(variogram, grid, n, Z, Y, &seed); /* make a log normal realization */ if (stat.type == 1 || stat.type == 2) { diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c index f4b6b15..4b58c82 100644 --- a/fftma_module/gen/lib_src/build_real.c +++ b/fftma_module/gen/lib_src/build_real.c @@ -20,7 +20,7 @@ /*realization: vector defining the real part */ /*ireal: vector defining the i-part */ -void build_real(int n[3], int NTOT, file_array_t* covar, file_array_t* realization, file_array_t* ireal, int cores) { +void build_real(int n[3], int NTOT, file_array_t* covar, file_array_t* realization, file_array_t* ireal) { int i, j, k, maille1; double temp; diff --git a/fftma_module/gen/lib_src/cardsin.c b/fftma_module/gen/lib_src/cardsin.c index b6a1ba7..e0b1ce6 100644 --- a/fftma_module/gen/lib_src/cardsin.c +++ b/fftma_module/gen/lib_src/cardsin.c @@ -4,7 +4,7 @@ #include /*cardsin covariance function*/ -double cardsin(double h, int cores) { +double cardsin(double h) { float delta = 20.371; double z; diff --git a/fftma_module/gen/lib_src/cgrid.c b/fftma_module/gen/lib_src/cgrid.c index 252ae35..e1b0086 100644 --- a/fftma_module/gen/lib_src/cgrid.c +++ b/fftma_module/gen/lib_src/cgrid.c @@ -10,7 +10,7 @@ /*n: vector with the number of cells along the */ /* 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) { +void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]) { int i, N; double D; @@ -30,7 +30,7 @@ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3], int cores D = grid.DZ; break; } - n[i] = length(N, i, variogram.scf, variogram.ap, D, variogram.Nvario, cores); + n[i] = length(N, i, variogram.scf, variogram.ap, D, variogram.Nvario); } } else { if ((n[0] < grid.NX) || (n[1] < grid.NY) || (n[2] < grid.NZ)) { diff --git a/fftma_module/gen/lib_src/clean_real.c b/fftma_module/gen/lib_src/clean_real.c index 4514d69..79ee8c0 100644 --- a/fftma_module/gen/lib_src/clean_real.c +++ b/fftma_module/gen/lib_src/clean_real.c @@ -8,7 +8,7 @@ #include #include -void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, file_array_t* vectorresult, struct realization_mod* realout, int cores) { +void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, file_array_t* vectorresult, struct realization_mod* realout) { int i, j, k, maille0, maille1; double NTOT; diff --git a/fftma_module/gen/lib_src/cov_value.c b/fftma_module/gen/lib_src/cov_value.c index 4e30f90..abe0c3f 100644 --- a/fftma_module/gen/lib_src/cov_value.c +++ b/fftma_module/gen/lib_src/cov_value.c @@ -4,7 +4,7 @@ #include /*selection of model covariance*/ -double cov_value(struct vario_mod variogram, double di, double dj, double dk, int cores) { +double cov_value(struct vario_mod variogram, double di, double dj, double dk) { double hx, hy, hz, h; double cov; int k; @@ -29,16 +29,16 @@ double cov_value(struct vario_mod variogram, double di, double dj, double dk, in cov += variogram.var[k] * spherical(h); break; case 4: - cov += variogram.var[k] * cardsin(h, cores); + cov += variogram.var[k] * cardsin(h); break; case 5: cov += variogram.var[k] * stable(h, variogram.alpha[k]); break; case 6: - cov += variogram.var[k] * gammf(h, variogram.alpha[k], cores); + cov += variogram.var[k] * gammf(h, variogram.alpha[k]); break; case 7: - cov += variogram.var[k] * cubic(h, cores); + cov += variogram.var[k] * cubic(h); break; case 8: cov += variogram.var[k] * nugget(h); diff --git a/fftma_module/gen/lib_src/covariance.c b/fftma_module/gen/lib_src/covariance.c index 0923d58..e02ecd9 100644 --- a/fftma_module/gen/lib_src/covariance.c +++ b/fftma_module/gen/lib_src/covariance.c @@ -4,7 +4,7 @@ /*builds the sampled covariance function*/ /*dimensions are even*/ -void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { +void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3]) { int i, j, k, maille, n2[3], symmetric; double di, dj, dk; @@ -22,7 +22,7 @@ void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod di = (double)i * mesh.DX; dj = (double)j * mesh.DY; dk = (double)k * mesh.DZ; - file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); + file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk)); if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) { /*area 2*/ @@ -38,7 +38,7 @@ void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod dj = (double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k); - file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); + file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk)); } if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { @@ -55,7 +55,7 @@ void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod dj = -(double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k); - file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); + file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk)); } if (k > 0 && k < n2[2]) { @@ -72,7 +72,7 @@ void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod dj = -(double)j * mesh.DY; dk = (double)k * mesh.DZ; maille = 1 + i + n[0] * (n[1] - j + n[1] * k); - file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); + file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk)); } if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { diff --git a/fftma_module/gen/lib_src/cubic.c b/fftma_module/gen/lib_src/cubic.c index ac36add..a371e76 100644 --- a/fftma_module/gen/lib_src/cubic.c +++ b/fftma_module/gen/lib_src/cubic.c @@ -4,7 +4,7 @@ #include /*cubic covariance function*/ -double cubic(double h, int cores) { +double cubic(double h) { double z; if (h >= 1.) { diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c index 9148f92..0bb86ea 100644 --- a/fftma_module/gen/lib_src/fftma2.c +++ b/fftma_module/gen/lib_src/fftma2.c @@ -23,7 +23,7 @@ /*output: */ /*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) { +void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct realization_mod* realin, struct realization_mod* realout, long* seed) { int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz; int solver; double temp; @@ -34,7 +34,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r axes(variogram.ap, variogram.scf, variogram.Nvario); /*pseudo-grid definition*/ - cgrid(variogram, grid, n, cores); + cgrid(variogram, grid, n); /*constant definition*/ NTOT = n[0] * n[1] * n[2]; @@ -63,26 +63,26 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r testmemory(worki); /*covariance function creation*/ - covariance(covar, variogram, grid, n, cores); + covariance(covar, variogram, grid, n); /*power spectrum*/ - fourt(covar, ireal, n, NDIM, 1, 0, workr, worki, cores); + fourt(covar, ireal, n, NDIM, 1, 0, workr, worki); /*organization of the input Gaussian white noise*/ solver = 0; - prebuild_gwn(grid, n, realin, realization, solver, cores, seed); + prebuild_gwn(grid, n, realin, realization, solver, seed); /*forward fourier transform of the GWN*/ - fourt(realization, ireal, n, NDIM, 1, 0, workr, worki, cores); + fourt(realization, ireal, n, NDIM, 1, 0, workr, worki); /* build realization in spectral domain */ - build_real(n, NTOT, covar, realization, ireal, cores); + build_real(n, NTOT, covar, realization, ireal); file_array_free(covar); remove("covar.txt"); /*backward fourier transform*/ - fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); + fourt(realization, ireal, n, NDIM, 0, 1, workr, worki); file_array_free(ireal); remove("ireal.txt"); @@ -91,5 +91,5 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r free(worki); /*output realization*/ - clean_real(realin, n, grid, realization, realout, cores); + clean_real(realin, n, grid, realization, realout); } diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c index bf50a13..5f3f79f 100644 --- a/fftma_module/gen/lib_src/fourt.c +++ b/fftma_module/gen/lib_src/fourt.c @@ -91,7 +91,7 @@ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* 10-06-2000, MLR */ -void fourt(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { +void fourt(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki) { 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; diff --git a/fftma_module/gen/lib_src/fourt_covar.c b/fftma_module/gen/lib_src/fourt_covar.c deleted file mode 100644 index db5724f..0000000 --- a/fftma_module/gen/lib_src/fourt_covar.c +++ /dev/null @@ -1,645 +0,0 @@ -#include -#include -#include -#include "file_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(file_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]; - } - - file_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++) { - file_array_get(datar, i, &valuei); - file_array_get(datar, i, &valuei1); - file_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; - - file_array_get(datar, i3, &valuei3); - file_array_get(datar, j3, &valuej3); - file_array_save(datar, i3, valuej3); - file_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]; - file_array_get(datar, j, &workr[i]); - worki[i] = datai[j]; - } else { - file_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) { - file_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; - - file_array_get(datar, i, &valuei); - file_array_get(datar, j, &valuej); - - file_array_save(datar, i, valuej - valuei); - file_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]; - file_array_get(datar, i, &tempi); - tempi = -tempi; - } else { - tempr = -datai[i]; - //tempi = datar[i]; - file_array_get(datar, i, &tempi); - } - - file_array_get(datar, j, &valuej); - file_array_save(datar, i, valuej - tempr); - file_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; - file_array_get(datar, i, &valuei); - file_array_get(datar, j, &valuej); - tempr = valuei * wr - datai[i] * wi; - tempi = valuei * wi + datai[i] * wr; - file_array_save(datar, i, valuej - tempr); - //datar[i] = valuej - tempr; - datai[i] = datai[j] - tempi; - file_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]; - file_array_get(datar, j, &sr); - si = datai[j]; - oldsr = 0.; - oldsi = 0.; - j -= ifp1; - L620: - stmpr = sr; - stmpi = si; - file_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]; - file_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; - file_array_get(datar, i, &valuei); - file_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; - file_array_save(datar, i, sumr + tempr); - //datar[i] = sumr + tempr; - datai[i] = difi + tempi; - file_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]; - file_array_get(datar, i, &valuei); - file_array_save(datar, j, valuei); - datai[j] = -datai[i]; - L755: - i++; - j--; - if (i < imax) - goto L750; - - file_array_get(datar, imin, &valueimin); - file_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]; - file_array_get(datar, i, &valuei); - file_array_save(datar, j, valuei); - datai[j] = datai[i]; - L770: - i--; - j--; - if (i > imin) - goto L765; - //datar[j] = datar[imin] + datai[imin]; - file_array_get(datar, imin, &valueimin); - file_array_save(datar, j, valueimin - datai[imin]); - datai[j] = 0.; - imax = imin; - goto L745; - L780: - file_array_get(datar, 1, &value1); - file_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]; - file_array_get(datar, j, &valuej); - file_array_save(datar, i, valuej); - datai[i] = -datai[j]; - j--; - } - } - j = jmax; - for (i = imin; i <= imax; i += np0) { - //datar[i] = datar[j]; - file_array_get(datar, j, &valuej); - file_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; -} diff --git a/fftma_module/gen/lib_src/gammf.c b/fftma_module/gen/lib_src/gammf.c index 90e5cea..b4ea73e 100644 --- a/fftma_module/gen/lib_src/gammf.c +++ b/fftma_module/gen/lib_src/gammf.c @@ -4,7 +4,7 @@ #include /*gamma covariance function*/ -double gammf(double h, double alpha, int cores) { +double gammf(double h, double alpha) { float delta; double z; diff --git a/fftma_module/gen/lib_src/gasdev.c b/fftma_module/gen/lib_src/gasdev.c index 7388087..f295bd5 100644 --- a/fftma_module/gen/lib_src/gasdev.c +++ b/fftma_module/gen/lib_src/gasdev.c @@ -4,19 +4,19 @@ #define NTAB 32 -double gasdev(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { +double gasdev(long* idum, long* idum2, long* iy, long iv[NTAB]) { /*returns a normally distributed deviate with 0 mean*/ /*and unit variance, using ran2(idum) as the source */ /*of uniform deviates */ - double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores); + double ran2(long* idum, long* idum2, long* iy, long iv[NTAB]); static int iset = 0; static double gset; double fac, rsq, v1, v2; if (iset == 0) { do { - v1 = 2.0 * ran2(idum, idum2, iy, iv, cores) - 1.0; - v2 = 2.0 * ran2(idum, idum2, iy, iv, cores) - 1.0; + v1 = 2.0 * ran2(idum, idum2, iy, iv) - 1.0; + v2 = 2.0 * ran2(idum, idum2, iy, iv) - 1.0; rsq = v1 * v1 + v2 * v2; } while (rsq >= 1.0 || rsq == 0.0); diff --git a/fftma_module/gen/lib_src/generate.c b/fftma_module/gen/lib_src/generate.c index 2cff9e2..8387b87 100644 --- a/fftma_module/gen/lib_src/generate.c +++ b/fftma_module/gen/lib_src/generate.c @@ -10,7 +10,7 @@ /*output: */ /* realization: structure defining the realization*/ -void generate(long* seed, int n, struct realization_mod* realization, int cores) { +void generate(long* seed, int n, struct realization_mod* realization) { int i; long idum2 = 123456789, iy = 0; long* iv; diff --git a/fftma_module/gen/lib_src/geostat.h b/fftma_module/gen/lib_src/geostat.h index 4bfcb71..dc69db2 100644 --- a/fftma_module/gen/lib_src/geostat.h +++ b/fftma_module/gen/lib_src/geostat.h @@ -300,7 +300,7 @@ struct realization_mod { void axes(double* ap, double* scf, int N); /*cardsin covariance value for lag h*/ -double cardsin(double h, int cores); +double cardsin(double h); /*Cholesky decomposition of matrix C */ /* C : symetric positive-definite matrix recorded */ @@ -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(file_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); +void covariance(file_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3]); /*computation of the covariance matrix for the well data*/ /*well coordinates are given as a number of cells */ @@ -358,10 +358,10 @@ void cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well, /*dj: distance along the Y axis */ /*dk: distance along the Z axis */ /* The returned value is the computed covariance value */ -double cov_value(struct vario_mod variogram, double di, double dj, double dk, int cores); +double cov_value(struct vario_mod variogram, double di, double dj, double dk); /*cubic covariance value for lag h*/ -double cubic(double h, int cores); +double cubic(double h); /*truncation of the power spectrum to remove */ /*high frequencies - isotropic case */ @@ -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(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); +void fourt(file_array_t* datar, file_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki); /*calculates F(x) = (1/a)*exp(-x*x/2)*/ double funtrun1(double x); @@ -413,7 +413,7 @@ double funtrun1(double x); float G(float x); /*gamma covariance value for lag h and exponent alpha*/ -double gammf(double h, double alpha, int cores); +double gammf(double h, double alpha); /*returns the value ln(G(x))*/ float gammln(float xx); @@ -425,7 +425,7 @@ float gammp(float a, float x); /*and unit variance, using ran1(idum) as the source */ /*of uniform deviates */ /*idum: seed */ -double gasdev(long* idum, long* idum2, long* iy, long* iv, int cores); +double gasdev(long* idum, long* idum2, long* iy, long* iv); /*gaussian covariance value for lag h*/ double gaussian(double h); @@ -480,7 +480,7 @@ void gradual(struct grad_mod grad, float* Zo, float* Z, float* Zfinal, int n, st /*n: vector with the number of cells along the */ /* 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); +void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3]); /*incomplete gamma function evaluated by its series*/ /*representation as gamser, also returns ln(G(a)) */ @@ -505,7 +505,7 @@ void krig_stat(float* b, int n, struct vario_mod variogram, struct welldata_mod /*i: considered direction */ /*scf: correlation length */ /*ap: normalized anisotropy axes */ -int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores); +int length(int N, int i, double* scf, double* ap, double D, int Nvari); /*calculates L.Z/ /* L : lower triangular matrix recorded */ @@ -519,7 +519,7 @@ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores void LtimeZ(double* L, float* Z, float* b, int n); /*determines the greatest prime factor of an integer*/ -int maxfactor(int n, int cores); +int maxfactor(int n); /*metrop returns a boolean varible that issues a */ /*verdict on whether to accept a reconfiguration */ @@ -545,7 +545,7 @@ double power(double h, double alpha); /*generates uniform deviates between 0 and 1*/ /*idum: seed */ -double ran2(long* idum, long* idum2, long* iy, long* iv, int cores); +double ran2(long* idum, long* idum2, long* iy, long* iv); /*calculates bt.b */ /* b : vector, bi, i = [0...n-1] */ diff --git a/fftma_module/gen/lib_src/length.c b/fftma_module/gen/lib_src/length.c index e47410a..d7f4ce0 100644 --- a/fftma_module/gen/lib_src/length.c +++ b/fftma_module/gen/lib_src/length.c @@ -2,8 +2,8 @@ #include /* compute the length for one dimension*/ -int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores) { - int maxfactor(int n, int cores); +int length(int N, int i, double* scf, double* ap, double D, int Nvari) { + int maxfactor(int n); double temp1, temp2; int n, j, k, nmax; int nlimit = 13; @@ -29,10 +29,10 @@ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores } if ((n % 2) != 0) n = n + 1; - nmax = maxfactor(n, cores); + nmax = maxfactor(n); while (nmax > nlimit) { n += 2; - nmax = maxfactor(n, cores); + nmax = maxfactor(n); } } diff --git a/fftma_module/gen/lib_src/maxfactor.c b/fftma_module/gen/lib_src/maxfactor.c index 95c285c..f950cf2 100644 --- a/fftma_module/gen/lib_src/maxfactor.c +++ b/fftma_module/gen/lib_src/maxfactor.c @@ -1,7 +1,7 @@ #include "genlib.h" /*determines the greatest prime factor of an integer*/ -int maxfactor(int n, int cores) { +int maxfactor(int n) { int test_fact(int* pnum, int fact, int* pmaxfac); int lnum, fact; int maxfac; diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c index d3bf7c0..4aeca86 100644 --- a/fftma_module/gen/lib_src/prebuild_gwn.c +++ b/fftma_module/gen/lib_src/prebuild_gwn.c @@ -21,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, file_array_t* realization, int solver, int cores, long* seed) { +void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, file_array_t* realization, int solver, long* seed) { int i, j, k, maille0, maille1; int ntot; @@ -36,7 +36,7 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin file_array_save(realization, 0, 0.); if (solver == 1) { for (i = 0; i < ntot; i++) { - double value = gasdev(seed, &idum2, &iy, iv, cores); + double value = gasdev(seed, &idum2, &iy, iv); file_array_save(realization, i+1, value); } } else { @@ -45,7 +45,7 @@ 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) { - double value = gasdev(seed, &idum2, &iy, iv, cores); + double value = gasdev(seed, &idum2, &iy, iv); file_array_save(realization, maille1, value); } else { file_array_save(realization, maille1, 0.); diff --git a/fftma_module/gen/lib_src/ran2.c b/fftma_module/gen/lib_src/ran2.c index cd9dcdd..4826971 100644 --- a/fftma_module/gen/lib_src/ran2.c +++ b/fftma_module/gen/lib_src/ran2.c @@ -16,7 +16,7 @@ #define EPS 1.2e-7 #define RNMX (1.0 - EPS) -double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) { +double ran2(long* idum, long* idum2, long* iy, long iv[NTAB]) { int j; long k; double temp; diff --git a/fftma_module/gen/moduleFFTMA.c b/fftma_module/gen/moduleFFTMA.c index fd5b07b..bb50691 100644 --- a/fftma_module/gen/moduleFFTMA.c +++ b/fftma_module/gen/moduleFFTMA.c @@ -31,12 +31,11 @@ static PyObject* genFunc(PyObject* self, PyObject* args) { struct statistic_mod stat; PyObject* out_array; npy_intp out_dims[NPY_MAXDIMS]; - int cores; - if (!Py_getvalues(args, &seed, &grid, &variogram, &stat, &cores)) + if (!Py_getvalues(args, &seed, &grid, &variogram, &stat)) return NULL; - Py_kgeneration(seed, grid, stat, variogram, &Z, &Y, n, cores); + Py_kgeneration(seed, grid, stat, variogram, &Z, &Y, n); out_dims[0] = grid.NZ; out_dims[1] = grid.NY; diff --git a/tests/stages/generation/test.py b/tests/stages/generation/test.py index 77683f3..94118b1 100644 --- a/tests/stages/generation/test.py +++ b/tests/stages/generation/test.py @@ -28,7 +28,7 @@ def generate(N): variance=3.5682389 typ=3 - k = gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8) + k = gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ) print(f"Generation with N = {N} time = {time.time() - start_time}s") return k