ssantisi
/
simulacion-permeabilidad
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

clean code base

Browse Source
improvement-file_array
Oli 3 years ago
parent 7d68fc075e
commit 40f5c5ba74
18 changed files with 270 additions and 478 deletions
Show Stats Download Patch File Download Diff File

36
fftma_module/gen/include/chunk_array.h
Unescape Escape View File

@ -1,36 +0,0 @@
#ifndef _CHUNKARRAY_H
#define _CHUNKARRAY_H
#include <stdarg.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#define MAX_CHUNK_SIZE 1500
typedef struct chunk_array {
size_t init_pos;
size_t chunk_size;
size_t total_size;
FILE* fp;
double* data;
}chunk_array_t;
chunk_array_t* chunk_array_create(char* filename, size_t total_size, size_t chunk_size);
void chunk_array_read(chunk_array_t* chunk_array);
//void chunk_array_write(chunk_array_t* chunk_array, char* filename);
void chunk_array_free(chunk_array_t* chunk_array);
bool chunk_array_get(chunk_array_t* chunk_array, size_t pos, double *value_ptr);
bool chunk_array_save(chunk_array_t* chunk_array, size_t pos, double valor);
void chunk_array_flush(chunk_array_t* chunk_array);
size_t chunk_array_size(chunk_array_t* chunk_array);
#endif

30
fftma_module/gen/include/file_array.h
Unescape Escape View File

@ -0,0 +1,30 @@
#ifndef _CHUNKARRAY_H
#define _CHUNKARRAY_H
#include <stdarg.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <stdbool.h>
typedef struct file_array {
size_t init_pos;
size_t total_size;
FILE* fp;
}file_array_t;
file_array_t* file_array_create(char* filename, size_t total_size);
void file_array_read(file_array_t* file_array);
void file_array_free(file_array_t* file_array);
bool file_array_get(file_array_t* file_array, size_t pos, double *value_ptr);
bool file_array_save(file_array_t* file_array, size_t pos, double valor);
void file_array_flush(file_array_t* file_array);
size_t file_array_size(file_array_t* file_array);
#endif

8
fftma_module/gen/include/geostat.h
Unescape Escape View File

@ -2,7 +2,7 @@
#include <stddef.h> #include <stddef.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include "chunk_array.h" #include "file_array.h"
#ifndef _GEOSTAT_H #ifndef _GEOSTAT_H
#define _GEOSTAT_H #define _GEOSTAT_H
@ -284,7 +284,7 @@ struct realization_mod {
int n; int n;
int code; int code;
double* vector; double* vector;
chunk_array_t* vector_2; file_array_t* vector_2;
}; };
/*=====================================================*/ /*=====================================================*/
@ -327,7 +327,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
/*variogram: structure defined above */ /*variogram: structure defined above */
/*grid: structure defined above */ /*grid: structure defined above */
/*n: number of gridblocks along X,Y and Z*/ /*n: number of gridblocks along X,Y and Z*/
void covariance(chunk_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], int cores);
/*computation of the covariance matrix for the well data*/ /*computation of the covariance matrix for the well data*/
/*well coordinates are given as a number of cells */ /*well coordinates are given as a number of cells */
@ -403,7 +403,7 @@ double exponential(double h);
/*workr: utility real part vector for storage */ /*workr: utility real part vector for storage */
/*worki: utility imaginary part vector for storage */ /*worki: utility imaginary part vector for storage */
/*The transformed data are returned in datar and datai*/ /*The transformed data are returned in datar and datai*/
void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); 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);
/*calculates F(x) = (1/a)*exp(-x*x/2)*/ /*calculates F(x) = (1/a)*exp(-x*x/2)*/
double funtrun1(double x); double funtrun1(double x);

8
fftma_module/gen/include/toolsFFTMA.h
Unescape Escape View File

@ -1,6 +1,6 @@
#include "genlib.h" #include "genlib.h"
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h" #include "file_array.h"
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
@ -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 */ /* must be a Gaussian white noise */
/*realization: structure defining a realization*/ /*realization: structure defining a realization*/
void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_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, int cores, long* seed);
/* build_real */ /* build_real */
/* build a realization in the spectral domain */ /* 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 */ /*realization: vector defining the real part */
/*ireal: vector defining the i-part */ /*ireal: vector defining the i-part */
void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, chunk_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 cores);
void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores); void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, file_array_t* vectorresult, struct realization_mod* realout, int cores);
#endif // define _TOOLSFFTMA_H #endif // define _TOOLSFFTMA_H

20
fftma_module/gen/lib_src/build_real.c
Unescape Escape View File

@ -6,7 +6,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <time.h> #include <time.h>
#include "chunk_array.h" #include "file_array.h"
/* build_real */ /* build_real */
/* build a realization in the spectral domain */ /* build a realization in the spectral domain */
@ -20,20 +20,20 @@
/*realization: vector defining the real part */ /*realization: vector defining the real part */
/*ireal: vector defining the i-part */ /*ireal: vector defining the i-part */
void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, chunk_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 cores) {
int i, j, k, maille1; int i, j, k, maille1;
double temp; double temp;
chunk_array_read(realization); file_array_read(realization);
chunk_array_read(ireal); file_array_read(ireal);
chunk_array_read(covar); file_array_read(covar);
/*decomposition and multiplication in the spectral domain*/ /*decomposition and multiplication in the spectral domain*/
for (k = 1; k <= n[2]; k++) { for (k = 1; k <= n[2]; k++) {
for (j = 1; j <= n[1]; j++) { for (j = 1; j <= n[1]; j++) {
for (i = 1; i <= n[0]; i++) { for (i = 1; i <= n[0]; i++) {
maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0];
chunk_array_get(covar, maille1, &temp); file_array_get(covar, maille1, &temp);
if (temp > 0.) { if (temp > 0.) {
temp = sqrt(temp) / (double)NTOT; temp = sqrt(temp) / (double)NTOT;
} else if (temp < 0.) { } else if (temp < 0.) {
@ -41,11 +41,11 @@ void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realiza
} }
double valuerealizationmaille1, valueirealmaille1; double valuerealizationmaille1, valueirealmaille1;
chunk_array_get(realization, maille1, &valuerealizationmaille1); file_array_get(realization, maille1, &valuerealizationmaille1);
chunk_array_get(ireal, maille1, &valueirealmaille1); file_array_get(ireal, maille1, &valueirealmaille1);
chunk_array_save(realization, maille1, temp * valuerealizationmaille1); file_array_save(realization, maille1, temp * valuerealizationmaille1);
chunk_array_save(ireal, maille1, temp * valueirealmaille1); file_array_save(ireal, maille1, temp * valueirealmaille1);
} }
} }
} }

92
fftma_module/gen/lib_src/chunk_array.c
Unescape Escape View File

@ -1,92 +0,0 @@
#include "chunk_array.h"
#include "stdbool.h"
void chunk_array_free(chunk_array_t* chunk_array) {
fclose(chunk_array->fp);
free(chunk_array->data);
free(chunk_array);
}
bool chunk_array_update_read(chunk_array_t* chunk_array, size_t pos) {
int init_pos = pos/chunk_array->chunk_size;
fseek(chunk_array->fp, init_pos * chunk_array->chunk_size * sizeof(double), SEEK_SET);
size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp);
chunk_array->init_pos += newLen;
}
/*
bool chunk_array_get(chunk_array_t* chunk_array, size_t pos, double *valor) {
if (pos>((chunk_array->init_pos + chunk_array->chunk_size)-1)) {
chunk_array_update_read(chunk_array, pos);
}
*valor=chunk_array->data[pos%chunk_array->chunk_size];
return true;
}
bool chunk_array_save(chunk_array_t* chunk_array, size_t pos, double valor) {
if (pos>((chunk_array->init_pos + chunk_array->chunk_size)-1)) {
chunk_array_flush(chunk_array);
}
chunk_array->data[pos%chunk_array->chunk_size]=valor;
return true;
}
*/
chunk_array_t* chunk_array_create(char* filename, size_t total_size, size_t chunk_size) {
chunk_array_t* chunk_array = (chunk_array_t*)malloc(sizeof(chunk_array_t));
chunk_array->fp = fopen(filename, "w+");
if (chunk_array == NULL || chunk_array->fp == NULL) {
return NULL;
}
chunk_array->data = malloc(chunk_size * sizeof(double));
if (chunk_size > 0 && chunk_array->data == NULL) {
free(chunk_array);
return NULL;
}
chunk_array->init_pos = 0;
chunk_array->chunk_size = chunk_size;
chunk_array->total_size = total_size;
return chunk_array;
}
void chunk_array_read(chunk_array_t* chunk_array) {
rewind(chunk_array->fp);
chunk_array->init_pos = 0;
size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp);
}
/*
void chunk_array_write(chunk_array_t* chunk_array, char* filename) {
chunk_array->fp = fopen(filename, "w");
if (chunk_array->fp == NULL) {
fclose(chunk_array->fp);
chunk_array->fp = fopen(filename, "w");
}
chunk_array->init_pos = 0;
}*/
void chunk_array_flush(chunk_array_t* chunk_array) {
size_t newLen = fwrite(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp);
chunk_array->init_pos += newLen;
}
bool chunk_array_get(chunk_array_t* chunk_array, size_t pos, double *valor) {
fseek(chunk_array->fp, pos * sizeof(double), SEEK_SET);
fread(valor, sizeof(double), 1, chunk_array->fp);
return true;
}
bool chunk_array_save(chunk_array_t* chunk_array, size_t pos, double valor) {
fseek(chunk_array->fp, pos * sizeof(double), SEEK_SET);
fwrite(&valor, sizeof(double), 1, chunk_array->fp);
return true;
}

8
fftma_module/gen/lib_src/clean_real.c
Unescape Escape View File

@ -1,5 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h" #include "file_array.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -8,7 +8,7 @@
#include <string.h> #include <string.h>
#include <time.h> #include <time.h>
void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_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 cores) {
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
double NTOT; double NTOT;
@ -16,7 +16,7 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
/*is the output realization already allocated?*/ /*is the output realization already allocated?*/
/*if not, memory allocation*/ /*if not, memory allocation*/
chunk_array_read(vectorresult); file_array_read(vectorresult);
if (realout->vector == NULL || realout->n != realin->n) { if (realout->vector == NULL || realout->n != realin->n) {
realout->vector = (double*)malloc(realin->n * sizeof(double)); realout->vector = (double*)malloc(realin->n * sizeof(double));
@ -35,7 +35,7 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
/* Modif du 18 juin 2003 */ /* Modif du 18 juin 2003 */
/*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/ /*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/
double valuemaille1; double valuemaille1;
chunk_array_get(vectorresult, maille1, &valuemaille1); file_array_get(vectorresult, maille1, &valuemaille1);
realout->vector[maille0] = valuemaille1; realout->vector[maille0] = valuemaille1;
} }
} }

30
fftma_module/gen/lib_src/covariance.c
Unescape Escape View File

@ -1,17 +1,17 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h" #include "file_array.h"
#include <time.h> #include <time.h>
/*builds the sampled covariance function*/ /*builds the sampled covariance function*/
/*dimensions are even*/ /*dimensions are even*/
void covariance(chunk_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 cores) {
int i, j, k, maille, n2[3], symmetric; int i, j, k, maille, n2[3], symmetric;
double di, dj, dk; double di, dj, dk;
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
n2[i] = n[i] / 2; n2[i] = n[i] / 2;
chunk_array_read(covar); file_array_read(covar);
for (i = 0; i <= n2[0]; i++) { for (i = 0; i <= n2[0]; i++) {
for (j = 0; j <= n2[1]; j++) { for (j = 0; j <= n2[1]; j++) {
@ -22,14 +22,14 @@ void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mo
di = (double)i * mesh.DX; di = (double)i * mesh.DX;
dj = (double)j * mesh.DY; dj = (double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
/*area 2*/ /*area 2*/
symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k)); symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k));
double value; double value;
chunk_array_get(covar, maille, &value); file_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value); file_array_save(covar, symmetric, value);
} }
if (i > 0 && i < n2[0]) { if (i > 0 && i < n2[0]) {
@ -38,15 +38,15 @@ void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mo
dj = (double)j * mesh.DY; dj = (double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k); maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k);
chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
/*area 8*/ /*area 8*/
symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k)); symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k));
double value; double value;
chunk_array_get(covar, maille, &value); file_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value); file_array_save(covar, symmetric, value);
} }
if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) { if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) {
@ -55,15 +55,15 @@ void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mo
dj = -(double)j * mesh.DY; dj = -(double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k); maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k);
chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2]) { if (k > 0 && k < n2[2]) {
/*area 6*/ /*area 6*/
symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k)); symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k));
double value; double value;
chunk_array_get(covar, maille, &value); file_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value); file_array_save(covar, symmetric, value);
} }
if (j > 0 && j < n2[1]) { if (j > 0 && j < n2[1]) {
@ -72,15 +72,15 @@ void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mo
dj = -(double)j * mesh.DY; dj = -(double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + i + n[0] * (n[1] - j + n[1] * k); maille = 1 + i + n[0] * (n[1] - j + n[1] * k);
chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores)); file_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
/*area 7*/ /*area 7*/
symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k)); symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k));
double value; double value;
chunk_array_get(covar, maille, &value); file_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value); file_array_save(covar, symmetric, value);
} }
} }
} }

14
fftma_module/gen/lib_src/fftma2.c
Unescape Escape View File

@ -1,5 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h" #include "file_array.h"
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
@ -27,7 +27,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz; int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz;
int solver; int solver;
double temp; double temp;
chunk_array_t *covar, *ireal, *realization; file_array_t *covar, *ireal, *realization;
double *workr, *worki; double *workr, *worki;
/*covariance axis normalization*/ /*covariance axis normalization*/
@ -52,9 +52,9 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
nxyz = NXYZ + 1; nxyz = NXYZ + 1;
/*array initialization*/ /*array initialization*/
covar = chunk_array_create("covar.txt", ntot, 1500); covar = file_array_create("covar.txt", ntot);
ireal = chunk_array_create("ireal.txt", ntot, 1500); ireal = file_array_create("ireal.txt", ntot);
realization = chunk_array_create("realization.txt", ntot, 1500); realization = file_array_create("realization.txt", ntot);
workr = (double*)malloc(nmax * sizeof(double)); workr = (double*)malloc(nmax * sizeof(double));
testmemory(workr); testmemory(workr);
@ -78,13 +78,13 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
/* build realization in spectral domain */ /* build realization in spectral domain */
build_real(n, NTOT, covar, realization, ireal, cores); build_real(n, NTOT, covar, realization, ireal, cores);
chunk_array_free(covar); file_array_free(covar);
remove("covar.txt"); remove("covar.txt");
/*backward fourier transform*/ /*backward fourier transform*/
fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores);
chunk_array_free(ireal); file_array_free(ireal);
remove("ireal.txt"); remove("ireal.txt");
free(workr); free(workr);

38
fftma_module/gen/lib_src/file_array.c
Unescape Escape View File

@ -0,0 +1,38 @@
#include "file_array.h"
#include "stdbool.h"
void file_array_free(file_array_t* file_array) {
fclose(file_array->fp);
free(file_array);
}
file_array_t* file_array_create(char* filename, size_t total_size) {
file_array_t* file_array = (file_array_t*)malloc(sizeof(file_array_t));
file_array->fp = fopen(filename, "w+");
if (file_array == NULL || file_array->fp == NULL) {
return NULL;
}
file_array->init_pos = 0;
file_array->total_size = total_size;
return file_array;
}
void file_array_read(file_array_t* file_array) {
rewind(file_array->fp);
file_array->init_pos = 0;
}
bool file_array_get(file_array_t* file_array, size_t pos, double *valor) {
fseek(file_array->fp, pos * sizeof(double), SEEK_SET);
fread(valor, sizeof(double), 1, file_array->fp);
return true;
}
bool file_array_save(file_array_t* file_array, size_t pos, double valor) {
fseek(file_array->fp, pos * sizeof(double), SEEK_SET);
fwrite(&valor, sizeof(double), 1, file_array->fp);
return true;
}

324
fftma_module/gen/lib_src/fourt.c
Unescape Escape View File

@ -1,7 +1,7 @@
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
#include "chunk_array.h" #include "file_array.h"
/*fast fourier transform */ /*fast fourier transform */
/* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */
@ -91,7 +91,7 @@
/* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */
/* 10-06-2000, MLR */ /* 10-06-2000, MLR */
void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { 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) {
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; 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 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; double valueri, valueri1, valueri3, valueii3, valuerj3, valueij3, valuerj, valueij, valueii, valuerimin, valueiimin, valuei1, valuer1;
@ -101,8 +101,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
ntot *= nn[idim]; ntot *= nn[idim];
} }
chunk_array_read(datar); file_array_read(datar);
chunk_array_read(datai); file_array_read(datai);
/*main loop for each dimension*/ /*main loop for each dimension*/
np1 = 1; np1 = 1;
@ -189,14 +189,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
ntot /= 2; ntot /= 2;
i = 1; i = 1;
for (j = 1; j <= ntot; j++) { for (j = 1; j <= ntot; j++) {
chunk_array_get(datar, i, &valueri); file_array_get(datar, i, &valueri);
////printf("[1] datar[%d] = %f\n", i, valueri); file_array_get(datar, i+1, &valueri1);
chunk_array_get(datar, i+1, &valueri1); file_array_save(datar, j, valueri);
////printf("[2] datar[%d] = %f\n", i+1, valueri1); file_array_save(datai, j, 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; i += 2;
} }
@ -214,28 +210,18 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (i3 = i1; i3 <= ntot; i3 += np2) { for (i3 = i1; i3 <= ntot; i3 += np2) {
j3 = j + i3 - i2; j3 = j + i3 - i2;
chunk_array_get(datar, i3, &valueri3); file_array_get(datar, i3, &valueri3);
chunk_array_get(datai, i3, &valueii3); file_array_get(datai, i3, &valueii3);
chunk_array_get(datar, j3, &valuerj3); file_array_get(datar, j3, &valuerj3);
chunk_array_get(datai, j3, &valueij3); file_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);
tempr = valueri3; tempr = valueri3;
tempi = valueii3; tempi = valueii3;
//printf("[50] Saving in datar the value %f in pos %d\n", valuerj3, i3); file_array_save(datar, i3, valuerj3);
//printf("[51] Saving in datai the value %f in pos %d\n", valueij3, i3); file_array_save(datai, i3, valueij3);
//printf("[52] Saving in datar the value %f in pos %d\n", tempr, j3); file_array_save(datar, j3, tempr);
//printf("[53] Saving in datai the value %f in pos %d\n", tempi, j3); file_array_save(datai, j3, tempi);
chunk_array_save(datar, i3, valuerj3);
chunk_array_save(datai, i3, valueij3);
chunk_array_save(datar, j3, tempr);
chunk_array_save(datai, j3, tempi);
} }
} }
@ -260,11 +246,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (i3 = i1; i3 <= ntot; i3 += np2) { for (i3 = i1; i3 <= ntot; i3 += np2) {
j = i3; j = i3;
for (i = 1; i <= n; i++) { for (i = 1; i <= n; i++) {
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, j, &valueij);
//printf("[7] datar[%d] = %f\n", j, valuerj);
//printf("[8] datai[%d] = %f\n", j, valueij);
if (icase != 3) { if (icase != 3) {
workr[i] = valuerj; workr[i] = valuerj;
@ -289,11 +272,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
i2max = i3 + np2 - np1; i2max = i3 + np2 - np1;
i = 1; i = 1;
for (i2 = i3; i2 <= i2max; i2 += np1) { for (i2 = i3; i2 <= i2max; i2 += np1) {
//printf("[54] Saving in datar the value %f in pos %d\n", workr[i], i2); file_array_save(datar, i2, workr[i]);
//printf("[55] Saving in datai the value %f in pos %d\n", worki[i], i2); file_array_save(datai, i2, worki[i]);
chunk_array_save(datar, i2, workr[i]);
chunk_array_save(datai, i2, worki[i]);
i++; i++;
} }
} }
@ -314,25 +294,15 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L310: L310:
j = i1; j = i1;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
chunk_array_get(datar, i, &tempr); file_array_get(datar, i, &tempr);
chunk_array_get(datai, i, &tempi); file_array_get(datai, i, &tempi);
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, j, &valueij);
//printf("[9] tempr = %f\n", i, tempr); file_array_save(datar, i, valuerj - tempr);
//printf("[10] tempi = %f\n", i, tempi); file_array_save(datai, i, valueij - tempi);
//printf("[11] datar[%d] = %f\n", j, valuerj); file_array_save(datar, j, valuerj + tempr);
//printf("[12] datai[%d] = %f\n", j, valueij); file_array_save(datai, j, valueij + tempi);
chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi);
chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi);
//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; j += istep;
} }
@ -352,44 +322,26 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
if (ifrwd != 0) { if (ifrwd != 0) {
chunk_array_get(datai, i, &tempr); file_array_get(datai, i, &tempr);
chunk_array_get(datar, i, &tempi); file_array_get(datar, i, &tempi);
//printf("[13] datai[%d] = %f\n", i, tempr);
//printf("[14] datar[%d] = %f\n", i, tempi);
tempi = -tempi; tempi = -tempi;
} else { } else {
chunk_array_get(datai, i, &tempr); file_array_get(datai, i, &tempr);
//printf("[15] datai[%d] = %f\n", i, tempr);
tempr = -tempr; tempr = -tempr;
chunk_array_get(datar, i, &tempi); file_array_get(datar, i, &tempi);
//printf("[16] datar[%d] = %f\n", i, tempi);
} }
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, 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);
chunk_array_get(datar, j, &valuerj); file_array_save(datar, i, valuerj - tempr);
chunk_array_get(datai, j, &valueij); file_array_save(datai, i, valueij - tempi);
//printf("[19] datar[%d] = %f\n", j, valuerj); file_array_get(datar, j, &valuerj);
//printf("[20] datai[%d] = %f\n", j, valueij); file_array_get(datai, j, &valueij);
chunk_array_save(datar, j, valuerj + tempr); file_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi); file_array_save(datai, j, valueij + tempi);
//printf("[62] Saving in datar the value %f in pos %d\n", valuerj + tempr, j);
//printf("[63] Saving in datai the value %f in pos %d\n", valueij + tempi, j);
j += istep; j += istep;
} }
@ -426,29 +378,18 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L510: L510:
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
chunk_array_get(datar, i, &valueri); file_array_get(datar, i, &valueri);
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, i, &valueii); file_array_get(datai, i, &valueii);
chunk_array_get(datai, j, &valueij); file_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);
tempr = valueri * wr - valueii * wi; tempr = valueri * wr - valueii * wi;
tempi = valueri * wi + valueii * wr; tempi = valueri * wi + valueii * wr;
chunk_array_save(datar, i, valuerj - tempr); file_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi); file_array_save(datai, i, valueij - tempi);
chunk_array_save(datar, j, valuerj + tempr); file_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi); file_array_save(datai, j, valueij + tempi);
//printf("[64] Saving in datar the value %f in pos %d\n", valuerj - tempr, i);
//printf("[65] Saving in datai the value %f in pos %d\n", valueij - tempi, i);
//printf("[66] Saving in datar the value %f in pos %d\n", valuerj + tempr, j);
//printf("[67] Saving in datai the value %f in pos %d\n", valueij + tempi, j);
j += istep; j += istep;
} }
imin = 2 * imin - i1; imin = 2 * imin - i1;
@ -501,11 +442,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
j3max = j2 + np2 - ifp2; j3max = j2 + np2 - ifp2;
for (j3 = j2; j3 <= j3max; j3 += ifp2) { for (j3 = j2; j3 <= j3max; j3 += ifp2) {
j = jmin + ifp2 - ifp1; j = jmin + ifp2 - ifp1;
chunk_array_get(datar, j, &sr); file_array_get(datar, j, &sr);
chunk_array_get(datai, j, &si); file_array_get(datai, j, &si);
//printf("[25] datar[%d] = %f\n", j, sr);
//printf("[26] datai[%d] = %f\n", j, si);
oldsr = 0.; oldsr = 0.;
oldsi = 0.; oldsi = 0.;
@ -514,11 +452,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
stmpr = sr; stmpr = sr;
stmpi = si; stmpi = si;
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, j, &valueij);
//printf("[27] datar[%d] = %f\n", j, valuerj);
//printf("[28] datai[%d] = %f\n", j, valueij);
sr = twowr * sr - oldsr + valuerj; sr = twowr * sr - oldsr + valuerj;
si = twowr * si - oldsi + valueij; si = twowr * si - oldsi + valueij;
@ -528,18 +463,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
if (j > jmin) if (j > jmin)
goto L620; goto L620;
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, j, &valueij);
workr[i] = wr * sr - wi * si - oldsr + valuerj; workr[i] = wr * sr - wi * si - oldsr + valuerj;
worki[i] = wi * sr + wr * si - oldsi + valueij; worki[i] = wi * sr + wr * si - oldsi + valueij;
//printf("[85] wr = %f, sr = %f, wi = %f, si = %f, oldsr = %f, datar[j] = %f\n", wr, sr, wi, si, oldsr, valuerj);
//printf("[86] wi = %f, sr = %f, wr = %f, si = %f, oldsi = %f, datai[j] = %f\n", wi, sr, wr, si, oldsi, valueij);
//printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i);
//printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i);
jmin += ifp2; jmin += ifp2;
i++; i++;
} }
@ -551,11 +480,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (j2 = i3; j2 <= j2max; j2 += ifp1) { for (j2 = i3; j2 <= j2max; j2 += ifp1) {
j3max = j2 + np2 - ifp2; j3max = j2 + np2 - ifp2;
for (j3 = j2; j3 <= j3max; j3 += 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); file_array_save(datar, j3, workr[i]);
//printf("[69] Saving in datai the value %f in pos %d, i = %d\n", worki[i], j3, i); file_array_save(datai, j3, worki[i]);
chunk_array_save(datar, j3, workr[i]);
chunk_array_save(datai, j3, worki[i]);
i++; i++;
} }
} }
@ -599,15 +525,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L710: L710:
j = jmin; j = jmin;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
chunk_array_get(datar, i, &valueri); file_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); file_array_get(datai, i, &valueii);
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_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);
sumr = (valueri + valuerj) / 2.; sumr = (valueri + valuerj) / 2.;
sumi = (valueii + valueij) / 2.; sumi = (valueii + valueij) / 2.;
@ -616,15 +537,10 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
tempr = wr * sumi + wi * difr; tempr = wr * sumi + wi * difr;
tempi = wi * sumi - wr * difr; tempi = wi * sumi - wr * difr;
chunk_array_save(datar, i, sumr + tempr); file_array_save(datar, i, sumr + tempr);
chunk_array_save(datai, i, difi + tempi); file_array_save(datai, i, difi + tempi);
chunk_array_save(datar, j, sumr - tempr); file_array_save(datar, j, sumr - tempr);
chunk_array_save(datai, j, tempi - difi); file_array_save(datai, j, tempi - difi);
//printf("[70] Saving in datar the value %f in pos %d\n", sumr + tempr, i);
//printf("[71] Saving in datai the value %f in pos %d\n", difi + tempi, i);
//printf("[72] Saving in datar the value %f in pos %d\n", sumr - tempr, j);
//printf("[73] Saving in datai the value %f in pos %d\n", tempi - difi, j);
j += np2; j += np2;
} }
@ -642,10 +558,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
if (ifrwd == 0) if (ifrwd == 0)
goto L740; goto L740;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
chunk_array_get(datai, i, &valueii); file_array_get(datai, i, &valueii);
//printf("[33] datai[%d] = %f\n", i, valueii); file_array_save(datai, i, -valueii);
chunk_array_save(datai, i, -valueii);
//printf("[73] Saving in datai the value %f in pos %d\n", -valueii, i);
} }
L740: L740:
np2 *= 2; np2 *= 2;
@ -657,86 +571,52 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
i = imin; i = imin;
goto L755; goto L755;
L750: ; L750: ;
chunk_array_get(datar, i, &valueri); file_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); file_array_get(datai, 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);
file_array_save(datar, j, valueri);
file_array_save(datai, j, -valueii);
L755: L755:
i++; i++;
j--; j--;
if (i < imax) if (i < imax)
goto L750; goto L750;
chunk_array_get(datar, imin, &valuerimin); file_array_get(datar, imin, &valuerimin);
chunk_array_get(datai, imin, &valueiimin); file_array_get(datai, imin, &valueiimin);
//printf("[36] datar[%d] = %f\n", imin, valuerimin); file_array_save(datar, j, valuerimin - valueiimin);
//printf("[37] datai[%d] = %f\n", imin, valueiimin); file_array_save(datai, j, 0.);
chunk_array_save(datar, j, valuerimin - valueiimin);
chunk_array_save(datai, j, 0.);
//printf("[75] Saving in datar the value %f in pos %d\n", valuerimin - valueiimin, j);
//printf("[76] Saving in datai the value %f in pos %d\n", 0., j);
if (i >= j) { if (i >= j) {
goto L780; goto L780;
} else { } else {
goto L770; goto L770;
} }
L765: L765:
chunk_array_get(datar, i, &valueri); file_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); file_array_get(datai, i, &valueii);
//printf("[38] datar[%d] = %f\n", i, valueri); file_array_save(datar, j, valueri);
//printf("[39] datai[%d] = %f\n", i, valueii); file_array_save(datai, j, 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);
L770: L770:
i--; i--;
j--; j--;
if (i > imin) if (i > imin)
goto L765; goto L765;
chunk_array_get(datar, imin, &valuerimin); file_array_get(datar, imin, &valuerimin);
chunk_array_get(datai, imin, &valueiimin); file_array_get(datai, 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);
file_array_save(datar, j, valuerimin + valueiimin);
file_array_save(datai, j, 0.);
imax = imin; imax = imin;
goto L745; goto L745;
L780: ; L780: ;
chunk_array_get(datai, 1, &valuei1); file_array_get(datai, 1, &valuei1);
chunk_array_get(datar, 1, &valuer1); file_array_get(datar, 1, &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);
file_array_save(datar, 1, valuei1 + valuer1);
file_array_save(datai, 1, 0.);
goto L900; goto L900;
/*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/
@ -754,35 +634,21 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
if (idim > 2) { if (idim > 2) {
j = jmax + np0; j = jmax + np0;
for (i = imin; i <= imax; i++) { for (i = imin; i <= imax; i++) {
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, 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);
file_array_save(datar, i, valuerj);
file_array_save(datai, i, -valueij);
j--; j--;
} }
} }
j = jmax; j = jmax;
for (i = imin; i <= imax; i += np0) { for (i = imin; i <= imax; i += np0) {
chunk_array_get(datar, j, &valuerj); file_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); file_array_get(datai, 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);
file_array_save(datar, i, valuerj);
file_array_save(datai, i, -valueij);
j -= np0; j -= np0;
} }
} }

94
fftma_module/gen/lib_src/fourt_covar.c
Unescape Escape View File

@ -1,7 +1,7 @@
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
#include "chunk_array.h" #include "file_array.h"
/*fast fourier transform */ /*fast fourier transform */
/* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */
@ -91,7 +91,7 @@
/* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */
/* 10-06-2000, MLR */ /* 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) { 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; 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 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; double value1, valuei, valuej, valuei1, valueimin, valuei3, valuej3;
@ -101,7 +101,7 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
ntot *= nn[idim]; ntot *= nn[idim];
} }
chunk_array_read(datar); file_array_read(datar);
/*main loop for each dimension*/ /*main loop for each dimension*/
np1 = 1; np1 = 1;
@ -188,9 +188,9 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
ntot /= 2; ntot /= 2;
i = 1; i = 1;
for (j = 1; j <= ntot; j++) { for (j = 1; j <= ntot; j++) {
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_get(datar, i, &valuei1); file_array_get(datar, i, &valuei1);
chunk_array_save(datar, j, valuei); file_array_save(datar, j, valuei);
//datar[j] = datar[i]; //datar[j] = datar[i];
datai[j] = valuei1; datai[j] = valuei1;
@ -217,10 +217,10 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
//datar[j3] = tempr; //datar[j3] = tempr;
datai[j3] = tempi; datai[j3] = tempi;
chunk_array_get(datar, i3, &valuei3); file_array_get(datar, i3, &valuei3);
chunk_array_get(datar, j3, &valuej3); file_array_get(datar, j3, &valuej3);
chunk_array_save(datar, i3, valuej3); file_array_save(datar, i3, valuej3);
chunk_array_save(datar, j3, valuei3); file_array_save(datar, j3, valuei3);
} }
} }
@ -247,10 +247,10 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
for (i = 1; i <= n; i++) { for (i = 1; i <= n; i++) {
if (icase != 3) { if (icase != 3) {
//workr[i] = datar[j]; //workr[i] = datar[j];
chunk_array_get(datar, j, &workr[i]); file_array_get(datar, j, &workr[i]);
worki[i] = datai[j]; worki[i] = datai[j];
} else { } else {
chunk_array_get(datar, j, &workr[i]); file_array_get(datar, j, &workr[i]);
//workr[i] = datar[j]; //workr[i] = datar[j];
worki[i] = 0.; worki[i] = 0.;
} }
@ -270,7 +270,7 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
i2max = i3 + np2 - np1; i2max = i3 + np2 - np1;
i = 1; i = 1;
for (i2 = i3; i2 <= i2max; i2 += np1) { for (i2 = i3; i2 <= i2max; i2 += np1) {
chunk_array_save(datar, i2, workr[i]); file_array_save(datar, i2, workr[i]);
//datar[i2] = workr[i]; //datar[i2] = workr[i];
datai[i2] = worki[i]; datai[i2] = worki[i];
i++; i++;
@ -300,11 +300,11 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
//datar[j] = datar[j] + tempr; //datar[j] = datar[j] + tempr;
datai[j] = datai[j] + tempi; datai[j] = datai[j] + tempi;
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
chunk_array_save(datar, i, valuej - valuei); file_array_save(datar, i, valuej - valuei);
chunk_array_save(datar, j, valuej + valuei); file_array_save(datar, j, valuej + valuei);
j += istep; j += istep;
} }
@ -326,17 +326,17 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
if (ifrwd != 0) { if (ifrwd != 0) {
tempr = datai[i]; tempr = datai[i];
//tempi = -datar[i]; //tempi = -datar[i];
chunk_array_get(datar, i, &tempi); file_array_get(datar, i, &tempi);
tempi = -tempi; tempi = -tempi;
} else { } else {
tempr = -datai[i]; tempr = -datai[i];
//tempi = datar[i]; //tempi = datar[i];
chunk_array_get(datar, i, &tempi); file_array_get(datar, i, &tempi);
} }
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
chunk_array_save(datar, i, valuej - tempr); file_array_save(datar, i, valuej - tempr);
chunk_array_save(datar, j, valuej - tempr); file_array_save(datar, j, valuej - tempr);
//datar[i] = datar[j] - tempr; //datar[i] = datar[j] - tempr;
datai[i] = datai[j] - tempi; datai[i] = datai[j] - tempi;
@ -379,14 +379,14 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
double valuei, valuej; double valuei, valuej;
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
tempr = valuei * wr - datai[i] * wi; tempr = valuei * wr - datai[i] * wi;
tempi = valuei * wi + datai[i] * wr; tempi = valuei * wi + datai[i] * wr;
chunk_array_save(datar, i, valuej - tempr); file_array_save(datar, i, valuej - tempr);
//datar[i] = valuej - tempr; //datar[i] = valuej - tempr;
datai[i] = datai[j] - tempi; datai[i] = datai[j] - tempi;
chunk_array_save(datar, i, valuej + tempr); file_array_save(datar, i, valuej + tempr);
//datar[j] += tempr; //datar[j] += tempr;
datai[j] += tempi; datai[j] += tempi;
j += istep; j += istep;
@ -442,7 +442,7 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
for (j3 = j2; j3 <= j3max; j3 += ifp2) { for (j3 = j2; j3 <= j3max; j3 += ifp2) {
j = jmin + ifp2 - ifp1; j = jmin + ifp2 - ifp1;
//sr = datar[j]; //sr = datar[j];
chunk_array_get(datar, j, &sr); file_array_get(datar, j, &sr);
si = datai[j]; si = datai[j];
oldsr = 0.; oldsr = 0.;
oldsi = 0.; oldsi = 0.;
@ -450,7 +450,7 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
L620: L620:
stmpr = sr; stmpr = sr;
stmpi = si; stmpi = si;
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
sr = twowr * sr - oldsr + valuej; sr = twowr * sr - oldsr + valuej;
si = twowr * si - oldsi + datai[j]; si = twowr * si - oldsi + datai[j];
oldsr = stmpr; oldsr = stmpr;
@ -472,7 +472,7 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
j3max = j2 + np2 - ifp2; j3max = j2 + np2 - ifp2;
for (j3 = j2; j3 <= j3max; j3 += ifp2) { for (j3 = j2; j3 <= j3max; j3 += ifp2) {
//datar[j3] = workr[i]; //datar[j3] = workr[i];
chunk_array_save(datar, j3, workr[i]); file_array_save(datar, j3, workr[i]);
datai[j3] = worki[i]; datai[j3] = worki[i];
i++; i++;
} }
@ -518,18 +518,18 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
j = jmin; j = jmin;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
double valuei, valuej; double valuei, valuej;
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
sumr = (valuei + valuej) / 2.; sumr = (valuei + valuej) / 2.;
sumi = (datai[i] + datai[j]) / 2.; sumi = (datai[i] + datai[j]) / 2.;
difr = (valuei - valuej) / 2.; difr = (valuei - valuej) / 2.;
difi = (datai[i] - datai[j]) / 2.; difi = (datai[i] - datai[j]) / 2.;
tempr = wr * sumi + wi * difr; tempr = wr * sumi + wi * difr;
tempi = wi * sumi - wr * difr; tempi = wi * sumi - wr * difr;
chunk_array_save(datar, i, sumr + tempr); file_array_save(datar, i, sumr + tempr);
//datar[i] = sumr + tempr; //datar[i] = sumr + tempr;
datai[i] = difi + tempi; datai[i] = difi + tempi;
chunk_array_save(datar, j, sumr - tempr); file_array_save(datar, j, sumr - tempr);
//datar[j] = sumr - tempr; //datar[j] = sumr - tempr;
datai[j] = tempi - difi; datai[j] = tempi - difi;
j += np2; j += np2;
@ -561,8 +561,8 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
goto L755; goto L755;
L750: L750:
//datar[j] = datar[i]; //datar[j] = datar[i];
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_save(datar, j, valuei); file_array_save(datar, j, valuei);
datai[j] = -datai[i]; datai[j] = -datai[i];
L755: L755:
i++; i++;
@ -570,8 +570,8 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
if (i < imax) if (i < imax)
goto L750; goto L750;
chunk_array_get(datar, imin, &valueimin); file_array_get(datar, imin, &valueimin);
chunk_array_save(datar, j, valueimin - datai[imin]); file_array_save(datar, j, valueimin - datai[imin]);
//datar[j] = datar[imin] - datai[imin]; //datar[j] = datar[imin] - datai[imin];
datai[j] = 0.; datai[j] = 0.;
if (i >= j) { if (i >= j) {
@ -581,8 +581,8 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
} }
L765: L765:
//datar[j] = datar[i]; //datar[j] = datar[i];
chunk_array_get(datar, i, &valuei); file_array_get(datar, i, &valuei);
chunk_array_save(datar, j, valuei); file_array_save(datar, j, valuei);
datai[j] = datai[i]; datai[j] = datai[i];
L770: L770:
i--; i--;
@ -590,14 +590,14 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
if (i > imin) if (i > imin)
goto L765; goto L765;
//datar[j] = datar[imin] + datai[imin]; //datar[j] = datar[imin] + datai[imin];
chunk_array_get(datar, imin, &valueimin); file_array_get(datar, imin, &valueimin);
chunk_array_save(datar, j, valueimin - datai[imin]); file_array_save(datar, j, valueimin - datai[imin]);
datai[j] = 0.; datai[j] = 0.;
imax = imin; imax = imin;
goto L745; goto L745;
L780: L780:
chunk_array_get(datar, 1, &value1); file_array_get(datar, 1, &value1);
chunk_array_save(datar, 1, value1 + datai[1]); file_array_save(datar, 1, value1 + datai[1]);
//datar[1] += datai[1]; //datar[1] += datai[1];
datai[1] = 0.; datai[1] = 0.;
goto L900; goto L900;
@ -618,8 +618,8 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
j = jmax + np0; j = jmax + np0;
for (i = imin; i <= imax; i++) { for (i = imin; i <= imax; i++) {
//datar[i] = datar[j]; //datar[i] = datar[j];
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
chunk_array_save(datar, i, valuej); file_array_save(datar, i, valuej);
datai[i] = -datai[j]; datai[i] = -datai[j];
j--; j--;
} }
@ -627,8 +627,8 @@ void fourt_covar(chunk_array_t* datar, double* datai, int nn[3], int ndim, int i
j = jmax; j = jmax;
for (i = imin; i <= imax; i += np0) { for (i = imin; i <= imax; i += np0) {
//datar[i] = datar[j]; //datar[i] = datar[j];
chunk_array_get(datar, j, &valuej); file_array_get(datar, j, &valuej);
chunk_array_save(datar, i, valuej); file_array_save(datar, i, valuej);
datai[i] = -datai[j]; datai[i] = -datai[j];
j -= np0; j -= np0;
} }

13
fftma_module/gen/lib_src/generate.c
Unescape Escape View File

@ -25,19 +25,6 @@ void generate(long* seed, int n, struct realization_mod* realization, int cores)
/*realization definition*/ /*realization definition*/
(*realization).n = n; (*realization).n = n;
(*realization).code = 0; (*realization).code = 0;
/*(*realization).vector_2 = chunk_array_create("realization1.txt", n, 512);
if ((*realization).vector_2 == NULL) {
log_error("RESULT = failed - No memory available in generate");
exit(1);
}*/
/*Gaussian white noise generation*/
/*for (i = 0; i < n; i++) {
double value = gasdev(seed, &idum2, &iy, iv, cores);
chunk_array_save((*realization).vector_2, i, value);
}
chunk_array_flush((*realization).vector_2);*/
free(iv); free(iv);
} }

8
fftma_module/gen/lib_src/geostat.h
Unescape Escape View File

@ -2,7 +2,7 @@
#include <stddef.h> #include <stddef.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include "chunk_array.h" #include "file_array.h"
#ifndef _GEOSTAT_H #ifndef _GEOSTAT_H
#define _GEOSTAT_H #define _GEOSTAT_H
@ -285,7 +285,7 @@ struct realization_mod {
int n; int n;
int code; int code;
double* vector; double* vector;
chunk_array_t* vector_2; file_array_t* vector_2;
}; };
/*=====================================================*/ /*=====================================================*/
@ -328,7 +328,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
/*variogram: structure defined above */ /*variogram: structure defined above */
/*grid: structure defined above */ /*grid: structure defined above */
/*n: number of gridblocks along X,Y and Z*/ /*n: number of gridblocks along X,Y and Z*/
void covariance(chunk_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], int cores);
/*computation of the covariance matrix for the well data*/ /*computation of the covariance matrix for the well data*/
/*well coordinates are given as a number of cells */ /*well coordinates are given as a number of cells */
@ -404,7 +404,7 @@ double exponential(double h);
/*workr: utility real part vector for storage */ /*workr: utility real part vector for storage */
/*worki: utility imaginary part vector for storage */ /*worki: utility imaginary part vector for storage */
/*The transformed data are returned in datar and datai*/ /*The transformed data are returned in datar and datai*/
void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); 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);
/*calculates F(x) = (1/a)*exp(-x*x/2)*/ /*calculates F(x) = (1/a)*exp(-x*x/2)*/
double funtrun1(double x); double funtrun1(double x);

12
fftma_module/gen/lib_src/prebuild_gwn.c
Unescape Escape View File

@ -1,5 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h" #include "file_array.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -21,7 +21,7 @@
/* must be a Gaussian white noise */ /* must be a Gaussian white noise */
/*realization: structure defining a realization*/ /*realization: structure defining a realization*/
void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_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, int cores, long* seed) {
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
int ntot; int ntot;
@ -33,11 +33,11 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
*seed = -(*seed); *seed = -(*seed);
ntot = n[0] * n[1] * n[2]; ntot = n[0] * n[1] * n[2];
chunk_array_save(realization, 0, 0.); file_array_save(realization, 0, 0.);
if (solver == 1) { if (solver == 1) {
for (i = 0; i < ntot; i++) { for (i = 0; i < ntot; i++) {
double value = gasdev(seed, &idum2, &iy, iv, cores); double value = gasdev(seed, &idum2, &iy, iv, cores);
chunk_array_save(realization, i+1, value); file_array_save(realization, i+1, value);
} }
} else { } else {
for (k = 1; k <= n[2]; k++) { for (k = 1; k <= n[2]; k++) {
@ -46,9 +46,9 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0];
if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) { 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, cores);
chunk_array_save(realization, maille1, value); file_array_save(realization, maille1, value);
} else { } else {
chunk_array_save(realization, maille1, 0.); file_array_save(realization, maille1, 0.);
} }
} }
} }

4
fftma_module/gen/moduleFFTMA.c
Unescape Escape View File

@ -7,7 +7,7 @@
#include "toolsFFTMA.h" #include "toolsFFTMA.h"
#include "toolsFFTPSIM.h" #include "toolsFFTPSIM.h"
#include "toolsIO.h" #include "toolsIO.h"
#include "chunk_array.h" #include "file_array.h"
#include <Python.h> #include <Python.h>
#include <math.h> #include <math.h>
#include <numpy/arrayobject.h> #include <numpy/arrayobject.h>
@ -56,7 +56,7 @@ static PyObject* genFunc(PyObject* self, PyObject* args) {
free(variogram.ap); free(variogram.ap);
free(variogram.vario); free(variogram.vario);
/*chunk_array_free(Z.vector_2); /*file_array_free(Z.vector_2);
remove("realization1.txt");*/ remove("realization1.txt");*/
return out_array; return out_array;

2
fftma_module/gen/setup.py
Unescape Escape View File

@ -42,7 +42,7 @@ module_FFTMA = Extension(
"./lib_src/clean_real.c", "./lib_src/clean_real.c",
"./lib_src/testmemory.c", "./lib_src/testmemory.c",
"./lib_src/genlib.c", "./lib_src/genlib.c",
"./lib_src/chunk_array.c" "./lib_src/file_array.c"
], ],
) )

7
fftma_module/gen/test.py
Unescape Escape View File

@ -1,6 +1,5 @@
from FFTMA import gen from FFTMA import gen
import numpy as np import numpy as np
import gc
import sys import sys
N=int(sys.argv[1]) N=int(sys.argv[1])
@ -29,8 +28,8 @@ variance=3.5682389
typ=3 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, 8)
np.save(f"out_{N}.npy",k) #np.save(f"out_{N}.npy",k)
k2 = np.load(f"out_{N}_ceci.npy")
del k print(k==k2)
gc.collect()
Write Preview
Loading…
Cancel
Save