remove logs and add tests

improvement-remove_generate_array
Oli 3 years ago
parent bfdc81ee4b
commit cfe4af4e15

@ -1,4 +1,3 @@
install: install:
source ./script_install.sh source ./script_install.sh
@ -11,9 +10,12 @@ fftma:
binaries: binaries:
./script_fortran.sh ./script_fortran.sh
test: binaries test: binaries fftma
cd tests/integration && python3 -m unittest test.py cd tests/integration && python3 -m unittest test.py
perf: binaries perf: binaries
cd tests/performance && python3 generation.py cd tests/performance && python3 generation.py
cd tests/performance && python3 connectivity.py cd tests/performance && python3 connectivity.py
test-gen: fftma
cd tests/stages/generation && python3 -m unittest test.py

@ -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])
@ -28,7 +27,5 @@ mean=15.3245987
variance=3.5682389 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)
np.save(f"out_{N}.npy",k) np.save(f"out_{N}.npy",k)
del k
gc.collect()

@ -11,5 +11,5 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
int Py_getvalues(PyObject*, long*, struct grid_mod*, struct vario_mod*, struct statistic_mod*, 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], int); void Py_kgeneration(long, struct grid_mod, struct statistic_mod, struct vario_mod, struct realization_mod*, struct realization_mod*, int[3]);

@ -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 512
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

@ -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 */ /* n: number of components in the vector */
/*output: */ /*output: */
/* realization: structure defining the realization*/ /* 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); void ikrig(int option, struct statfacies_mod facies, int kk, struct welldata_mod data, struct variotable_mod variogram, struct grid_mod grid, float* krigout);

@ -2,7 +2,6 @@
#include <stddef.h> #include <stddef.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include "chunk_array.h"
#ifndef _GEOSTAT_H #ifndef _GEOSTAT_H
#define _GEOSTAT_H #define _GEOSTAT_H
@ -284,7 +283,6 @@ struct realization_mod {
int n; int n;
int code; int code;
double* vector; double* vector;
chunk_array_t* vector_2;
}; };
/*=====================================================*/ /*=====================================================*/
@ -299,7 +297,7 @@ struct realization_mod {
void axes(double* ap, double* scf, int N); void axes(double* ap, double* scf, int N);
/*cardsin covariance value for lag h*/ /*cardsin covariance value for lag h*/
double cardsin(double h, int cores); double cardsin(double h);
/*Cholesky decomposition of matrix C */ /*Cholesky decomposition of matrix C */
/* C : symetric positive-definite matrix recorded */ /* C : symetric positive-definite matrix recorded */
@ -327,7 +325,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(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3]);
/*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 */
@ -357,10 +355,10 @@ void cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well,
/*dj: distance along the Y axis */ /*dj: distance along the Y axis */
/*dk: distance along the Z axis */ /*dk: distance along the Z axis */
/* The returned value is the computed covariance value */ /* 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*/ /*cubic covariance value for lag h*/
double cubic(double h, int cores); double cubic(double h);
/*truncation of the power spectrum to remove */ /*truncation of the power spectrum to remove */
/*high frequencies - isotropic case */ /*high frequencies - isotropic case */
@ -403,7 +401,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(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki);
/*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);
@ -412,7 +410,7 @@ double funtrun1(double x);
float G(float x); float G(float x);
/*gamma covariance value for lag h and exponent alpha*/ /*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))*/ /*returns the value ln(G(x))*/
float gammln(float xx); float gammln(float xx);
@ -424,7 +422,7 @@ float gammp(float a, float x);
/*and unit variance, using ran1(idum) as the source */ /*and unit variance, using ran1(idum) as the source */
/*of uniform deviates */ /*of uniform deviates */
/*idum: seed */ /*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*/ /*gaussian covariance value for lag h*/
double gaussian(double h); double gaussian(double h);
@ -479,7 +477,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 */ /*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */ /* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */ /* 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*/ /*incomplete gamma function evaluated by its series*/
/*representation as gamser, also returns ln(G(a)) */ /*representation as gamser, also returns ln(G(a)) */
@ -504,7 +502,7 @@ void krig_stat(float* b, int n, struct vario_mod variogram, struct welldata_mod
/*i: considered direction */ /*i: considered direction */
/*scf: correlation length */ /*scf: correlation length */
/*ap: normalized anisotropy axes */ /*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/ /*calculates L.Z/
/* L : lower triangular matrix recorded */ /* L : lower triangular matrix recorded */
@ -518,7 +516,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); void LtimeZ(double* L, float* Z, float* b, int n);
/*determines the greatest prime factor of an integer*/ /*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 */ /*metrop returns a boolean varible that issues a */
/*verdict on whether to accept a reconfiguration */ /*verdict on whether to accept a reconfiguration */
@ -544,7 +542,7 @@ double power(double h, double alpha);
/*generates uniform deviates between 0 and 1*/ /*generates uniform deviates between 0 and 1*/
/*idum: seed */ /*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 */ /*calculates bt.b */
/* b : vector, bi, i = [0...n-1] */ /* b : vector, bi, i = [0...n-1] */

@ -36,7 +36,7 @@
/*realout4: structure defining a yvelocity field */ /*realout4: structure defining a yvelocity field */
/*realout5: structure defining a zvelocity 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 */ /* prebuild_gwn */
/* Produce a first construction in real space of the Gaussian white noise */ /* Produce a first construction in real space of the Gaussian white noise */
@ -51,7 +51,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, double* realization, int solver, int cores, long* seed); void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, long* seed);
/* build_real */ /* build_real */
/* build a realization in the spectral domain */ /* build a realization in the spectral domain */
@ -65,7 +65,7 @@ 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, double* covar, double* realization, double* ireal, int cores); void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal);
void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout); void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout);

@ -1,6 +1,6 @@
#include "Py_py-api.h" #include "Py_py-api.h"
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "geostat.h" #include "geostat.h"
#include "pressure.h" #include "pressure.h"
#include "simpio.h" #include "simpio.h"
@ -42,14 +42,10 @@
#endif #endif
#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();
log_info("RESULT = in progress");
int i, varioNargs = 12, j = 0; int i, varioNargs = 12, j = 0;
PyObject* listvario; PyObject* listvario;
PyObject* vgr; PyObject* vgr;
//char* gwnfilename;
stat->nblock_mean = 1; stat->nblock_mean = 1;
stat->nblock_var = 1; stat->nblock_var = 1;
@ -59,11 +55,10 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
stat->variance = (double*)malloc(stat->nblock_var * sizeof(double)); stat->variance = (double*)malloc(stat->nblock_var * sizeof(double));
if (stat->variance == NULL) { if (stat->variance == NULL) {
free(stat->mean); free(stat->mean);
log_error("RESULT = failed");
return 0; return 0;
} }
if (!PyArg_ParseTuple(args, "iiidddlO!ddii", /*"iiidddslO!ddi",*/ if (!PyArg_ParseTuple(args, "iiidddlO!ddi", /*"iiidddslO!ddi",*/
&(grid->NX), &(grid->NX),
&(grid->NY), &(grid->NY),
&(grid->NZ), &(grid->NZ),
@ -76,11 +71,9 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
&listvario, &listvario,
stat->mean + 0, stat->mean + 0,
stat->variance + 0, stat->variance + 0,
&(stat->type), &(stat->type))) {
cores)) {
free(stat->mean); free(stat->mean);
free(stat->variance); free(stat->variance);
log_error("RESULT = failed");
return 0; return 0;
} }
@ -90,7 +83,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
if (variogram->var == NULL) { if (variogram->var == NULL) {
free(stat->mean); free(stat->mean);
free(stat->variance); free(stat->variance);
log_error("RESULT = failed");
return 0; return 0;
} }
variogram->vario = (int*)malloc(variogram->Nvario * sizeof(int)); variogram->vario = (int*)malloc(variogram->Nvario * sizeof(int));
@ -98,7 +90,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
free(stat->mean); free(stat->mean);
free(stat->variance); free(stat->variance);
free(variogram->var); free(variogram->var);
log_error("RESULT = failed");
return 0; return 0;
} }
variogram->alpha = (double*)malloc(variogram->Nvario * sizeof(double)); variogram->alpha = (double*)malloc(variogram->Nvario * sizeof(double));
@ -107,7 +98,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
free(stat->variance); free(stat->variance);
free(variogram->var); free(variogram->var);
free(variogram->vario); free(variogram->vario);
log_error("RESULT = failed");
return 0; return 0;
} }
variogram->scf = (double*)malloc(3 * variogram->Nvario * sizeof(double)); variogram->scf = (double*)malloc(3 * variogram->Nvario * sizeof(double));
@ -117,7 +107,6 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
free(variogram->var); free(variogram->var);
free(variogram->vario); free(variogram->vario);
free(variogram->alpha); free(variogram->alpha);
log_error("RESULT = failed");
return 0; return 0;
} }
variogram->ap = (double*)malloc(9 * variogram->Nvario * sizeof(double)); variogram->ap = (double*)malloc(9 * variogram->Nvario * sizeof(double));
@ -128,13 +117,11 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
free(variogram->vario); free(variogram->vario);
free(variogram->alpha); free(variogram->alpha);
free(variogram->scf); free(variogram->scf);
log_error("RESULT = failed");
return 0; return 0;
} }
for (i = 0; i < variogram->Nvario; i++) { for (i = 0; i < variogram->Nvario; i++) {
vgr = PyList_GetItem(listvario, i); vgr = PyList_GetItem(listvario, i);
if (PyTuple_Size(vgr) != 12) { if (PyTuple_Size(vgr) != 12) {
log_error("RESULT = failed");
return 0; return 0;
} }
(variogram->var)[i] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); (variogram->var)[i] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++));
@ -150,10 +137,5 @@ int Py_getvalues(PyObject* args, long* seed, struct grid_mod* grid, struct vario
(variogram->ap)[i * 9 + 4] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); (variogram->ap)[i * 9 + 4] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++));
(variogram->ap)[i * 9 + 5] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++)); (variogram->ap)[i * 9 + 5] = PyFloat_AsDouble(PyTuple_GetItem(vgr, j++));
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
log_info("RESULT = success, ELAPSED = %f seconds", time_taken);
return 1; return 1;
} }

@ -4,8 +4,8 @@
#include "simpio.h" #include "simpio.h"
#include "toolsFFTMA.h" #include "toolsFFTMA.h"
#include "toolsIO.h" #include "toolsIO.h"
#include "log.h"
#include "memory.h"
#include <Python.h> #include <Python.h>
#include <numpy/arrayobject.h> #include <numpy/arrayobject.h>
#include <stdarg.h> #include <stdarg.h>
@ -19,12 +19,7 @@
/* Z is the GWN with 0-mean and 1-variance */ /* Z is the GWN with 0-mean and 1-variance */
/* Y is the realization with mean and variance wanted */ /* 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]) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, N; int i, N;
int typelog; int typelog;
@ -34,47 +29,14 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat,
n[1] = 0; n[1] = 0;
n[2] = 0; n[2] = 0;
log_info("RESULT = in progress, N = %d", N); generate(&seed, N, Z);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
generate(&seed, N, Z, cores);
/*FFTMA*/ /*FFTMA*/
printf("pre fftma2\n"); FFTMA2(variogram, grid, n, Z, Y, &seed);
FFTMA2(variogram, grid, n, Z, Y, cores, &seed);
printf("post fftma2\n");
/* make a log normal realization */ /* make a log normal realization */
if (stat.type == 1 || stat.type == 2) { if (stat.type == 1 || stat.type == 2) {
typelog = stat.type + 2; typelog = stat.type + 2;
nor2log(Y, typelog, Y); nor2log(Y, typelog, Y);
} }
printf("termino nor2log\n");
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
printf("termino pykgeneration\n");
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,6 +1,6 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -21,23 +21,10 @@
/*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, double* covar, double* realization, double* ireal, int cores) { void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, j, k, maille1; int i, j, k, maille1;
double temp; double temp;
log_info("RESULT = in progress, NTOT = %d, covar = %f, n[0] = %d, n[1] = %d, n[2] = %d", NTOT, *covar, n[0], n[1], n[2]);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
/*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++) {
@ -54,23 +41,4 @@ void build_real(int n[3], int NTOT, double* covar, double* realization, double*
} }
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, realization = %f, ireal = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", *realization, *ireal, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,26 +1,12 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
/*cardsin covariance function*/ /*cardsin covariance function*/
double cardsin(double h, int cores) { double cardsin(double h) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, h = %f", h);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
float delta = 20.371; float delta = 20.371;
double z; double z;
@ -30,25 +16,5 @@ double cardsin(double h, int cores) {
} else { } else {
z = 1.; z = 1.;
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", z, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return z; return z;
} }

@ -1,6 +1,6 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
@ -12,23 +12,9 @@
/*n: vector with the number of cells along the */ /*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */ /* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */ /* 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]) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, N; int i, N;
double D; double D;
log_info("RESULT = in progress");
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
if (n == NULL || n[0] == 0 || n[1] == 0 || n[2] == 0) { if (n == NULL || n[0] == 0 || n[1] == 0 || n[2] == 0) {
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
@ -46,31 +32,11 @@ void cgrid(struct vario_mod variogram, struct grid_mod grid, int n[3], int cores
D = grid.DZ; D = grid.DZ;
break; 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 { } else {
if ((n[0] < grid.NX) || (n[1] < grid.NY) || (n[2] < grid.NZ)) { if ((n[0] < grid.NX) || (n[1] < grid.NY) || (n[2] < grid.NZ)) {
log_error("RESULT = failed - Indicated dimensions are inappropriate in cgrid");
exit; exit;
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, n[0] = %d, n[1] = %d, n[2] = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", n[0], n[1], n[2], time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,73 +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) {
printf("Llame a chunk array update\n");
size_t newLen = fread(chunk_array->data, sizeof(double), chunk_array->chunk_size, chunk_array->fp);
chunk_array->init_pos += newLen;
}
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);
}
*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) {
printf("ke");
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;
}

@ -1,6 +1,6 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -9,12 +9,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, double* vectorresult, struct realization_mod* realout, int cores) { void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
double NTOT; double NTOT;
@ -22,19 +17,9 @@ 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*/
log_info("RESULT = in progress, NTOT = %f", NTOT);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
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));
if (realout->vector == NULL) { if (realout->vector == NULL) {
log_error("RESULT = failed - No memory available");
exit; exit;
} }
} }
@ -52,23 +37,4 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
} }
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,26 +1,12 @@
#include "genlib.h" #include "genlib.h"
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <time.h> #include <time.h>
/*selection of model covariance*/ /*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* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, di = %f, dj = %f, dk = %f", di, dj, dk);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
double hx, hy, hz, h; double hx, hy, hz, h;
double cov; double cov;
int k; int k;
@ -45,16 +31,16 @@ double cov_value(struct vario_mod variogram, double di, double dj, double dk, in
cov += variogram.var[k] * spherical(h); cov += variogram.var[k] * spherical(h);
break; break;
case 4: case 4:
cov += variogram.var[k] * cardsin(h, cores); cov += variogram.var[k] * cardsin(h);
break; break;
case 5: case 5:
cov += variogram.var[k] * stable(h, variogram.alpha[k]); cov += variogram.var[k] * stable(h, variogram.alpha[k]);
break; break;
case 6: case 6:
cov += variogram.var[k] * gammf(h, variogram.alpha[k], cores); cov += variogram.var[k] * gammf(h, variogram.alpha[k]);
break; break;
case 7: case 7:
cov += variogram.var[k] * cubic(h, cores); cov += variogram.var[k] * cubic(h);
break; break;
case 8: case 8:
cov += variogram.var[k] * nugget(h); cov += variogram.var[k] * nugget(h);
@ -64,25 +50,5 @@ double cov_value(struct vario_mod variogram, double di, double dj, double dk, in
break; break;
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, cov = %f, hx = %f, hy = %f, hz = %f, h = %f , ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", cov, hx, hy, hz, h, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return cov; return cov;
} }

@ -1,28 +1,14 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.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(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3]) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, j, k, maille, n2[3], symmetric; int i, j, k, maille, n2[3], symmetric;
double di, dj, dk; double di, dj, dk;
log_info("RESULT = in progress");
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
n2[i] = n[i] / 2; n2[i] = n[i] / 2;
@ -35,7 +21,7 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
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;
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); 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]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
/*area 2*/ /*area 2*/
@ -49,7 +35,7 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
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);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); covar[maille] = (double)cov_value(variogram, di, dj, dk);
} }
if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
@ -64,7 +50,7 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
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);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); covar[maille] = (double)cov_value(variogram, di, dj, dk);
} }
if (k > 0 && k < n2[2]) { if (k > 0 && k < n2[2]) {
@ -79,7 +65,7 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
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);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); covar[maille] = (double)cov_value(variogram, di, dj, dk);
} }
if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
@ -90,23 +76,4 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
} }
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, di = %f, dj = %f, dk = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", di, dj, dk, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,26 +1,12 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
/*cubic covariance function*/ /*cubic covariance function*/
double cubic(double h, int cores) { double cubic(double h) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, h = %f", h);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
double z; double z;
if (h >= 1.) { if (h >= 1.) {
@ -28,25 +14,5 @@ double cubic(double h, int cores) {
} else { } else {
z = 1. - 7. * (double)(h * h) + (35. / 4.) * (double)(h * h * h) - 3.5 * (double)(h * h * h * h * h) + .75 * (double)(h * h * h * h * h * h * h); z = 1. - 7. * (double)(h * h) + (35. / 4.) * (double)(h * h * h) - 3.5 * (double)(h * h * h * h * h) + .75 * (double)(h * h * h * h * h * h * h);
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", z, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return z; return z;
} }

@ -1,11 +1,9 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
/*exponential covariance function*/ /*exponential covariance function*/
double exponential(double h) { double exponential(double h) {
log_info("RESULT = in progress, h = %f", h);
return (exp(-3. * (double)h)); return (exp(-3. * (double)h));
} }

@ -1,6 +1,6 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
@ -24,21 +24,7 @@
/*output: */ /*output: */
/*realout: structure defining a realization - */ /*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) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress");
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
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;
@ -48,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); axes(variogram.ap, variogram.scf, variogram.Nvario);
/*pseudo-grid definition*/ /*pseudo-grid definition*/
cgrid(variogram, grid, n, cores); cgrid(variogram, grid, n);
/*constant definition*/ /*constant definition*/
NTOT = n[0] * n[1] * n[2]; NTOT = n[0] * n[1] * n[2];
@ -82,51 +68,30 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
testmemory(worki); testmemory(worki);
/*covariance function creation*/ /*covariance function creation*/
covariance(covar, variogram, grid, n, cores); covariance(covar, variogram, grid, n);
/*power spectrum*/ /*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*/ /*organization of the input Gaussian white noise*/
printf("pre prebuild_gwn\n");
solver = 0; solver = 0;
prebuild_gwn(grid, n, realin, realization, solver, cores, seed); prebuild_gwn(grid, n, realin, realization, solver, seed);
printf("post prebuild_gwn\n");
/*forward fourier transform of the GWN*/ /*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 realization in spectral domain */
build_real(n, NTOT, covar, realization, ireal, cores); build_real(n, NTOT, covar, realization, ireal);
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
free(covar); free(covar);
/*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);
free(ireal); free(ireal);
free(workr); free(workr);
free(worki); free(worki);
/*output realization*/ /*output realization*/
clean_real(realin, n, grid, realization, realout, cores); clean_real(realin, n, grid, realization, realout);
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
log_info("RESULT = success, NTOT = %d, NMAX = %d, NDIM = %d, ntot = %d, nmax = %d, NXYZ = %d, nxyz = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", NTOT, NMAX, NDIM, ntot, nmax, NXYZ, nxyz, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,5 +1,5 @@
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
@ -92,21 +92,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(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress");
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
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;
@ -599,22 +585,5 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
nprev = n; nprev = n;
} }
L920: L920:
t = clock() - t; return;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
log_info("RESULT = success, ELAPSED = %f, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,50 +1,16 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
/*gamma covariance function*/ /*gamma covariance function*/
double gammf(double h, double alpha, int cores) { double gammf(double h, double alpha) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, h = %f, alpha = %f", h, alpha);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
float delta; float delta;
double z; double z;
delta = pow(20., 1. / alpha) - 1.; delta = pow(20., 1. / alpha) - 1.;
z = 1. / (double)(pow(1. + h * delta, alpha)); z = 1. / (double)(pow(1. + h * delta, alpha));
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
log_info("RESULT = success, delta = %f, z = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", delta, z, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return z; return z;
} }

@ -1,87 +1,33 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <time.h> #include <time.h>
#define NTAB 32 #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*/ /*returns a normally distributed deviate with 0 mean*/
/*and unit variance, using ran2(idum) as the source */ /*and unit variance, using ran2(idum) as the source */
/*of uniform deviates */ /*of uniform deviates */
double* used_ram_t0 = malloc(sizeof(double)); double ran2(long* idum, long* idum2, long* iy, long iv[NTAB]);
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, idum = %f, idum2 = %f, iy = %f", *idum, *idum2, *iy);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores);
static int iset = 0; static int iset = 0;
static double gset; static double gset;
double fac, rsq, v1, v2; double fac, rsq, v1, v2;
if (iset == 0) { if (iset == 0) {
do { do {
v1 = 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, cores) - 1.0; v2 = 2.0 * ran2(idum, idum2, iy, iv) - 1.0;
rsq = v1 * v1 + v2 * v2; rsq = v1 * v1 + v2 * v2;
} while (rsq >= 1.0 || rsq == 0.0); } while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log(rsq) / rsq); fac = sqrt(-2.0 * log(rsq) / rsq);
gset = v1 * fac; gset = v1 * fac;
iset = 1; iset = 1;
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, gset = %f, fac = %f, v1 = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", gset, fac, v1, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return (v2 * fac); return (v2 * fac);
} else { } else {
iset = 0; iset = 0;
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, gset = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", gset, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return gset; return gset;
} }
} }

@ -1,10 +1,9 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
/*gaussian covariance function*/ /*gaussian covariance function*/
double gaussian(double h) { double gaussian(double h) {
log_info("RESULT = in progress, h = %f", h);
return (exp(-3. * (double)(h * h))); return (exp(-3. * (double)(h * h)));
} }

@ -2,8 +2,8 @@
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
#include "log.h"
#include "memory.h"
/* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */ /* GENERATION OF A GAUSSIAN WHITE NOISE VECTOR */
/*input: */ /*input: */
@ -12,28 +12,11 @@
/*output: */ /*output: */
/* realization: structure defining the realization*/ /* 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) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
log_info("RESULT = in progress, n = %d", n);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
clock_t t = clock();
int i; int i;
long idum2 = 123456789, iy = 0; long idum2 = 123456789, iy = 0;
long* iv;
int iset = 0; int iset = 0;
iv = (long*)malloc(NTAB * sizeof(long));
/*negative seed*/ /*negative seed*/
if (*seed > 0.0) if (*seed > 0.0)
*seed = -(*seed); *seed = -(*seed);
@ -41,38 +24,4 @@ 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);*/
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
free(iv);
log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -2,7 +2,6 @@
#include <stddef.h> #include <stddef.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include "chunk_array.h"
#ifndef _GEOSTAT_H #ifndef _GEOSTAT_H
#define _GEOSTAT_H #define _GEOSTAT_H
@ -285,7 +284,6 @@ struct realization_mod {
int n; int n;
int code; int code;
double* vector; double* vector;
chunk_array_t* vector_2;
}; };
/*=====================================================*/ /*=====================================================*/
@ -300,7 +298,7 @@ struct realization_mod {
void axes(double* ap, double* scf, int N); void axes(double* ap, double* scf, int N);
/*cardsin covariance value for lag h*/ /*cardsin covariance value for lag h*/
double cardsin(double h, int cores); double cardsin(double h);
/*Cholesky decomposition of matrix C */ /*Cholesky decomposition of matrix C */
/* C : symetric positive-definite matrix recorded */ /* C : symetric positive-definite matrix recorded */
@ -328,7 +326,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(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3]);
/*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 */
@ -358,10 +356,10 @@ void cov_matrix(double* C, struct vario_mod variogram, struct welldata_mod well,
/*dj: distance along the Y axis */ /*dj: distance along the Y axis */
/*dk: distance along the Z axis */ /*dk: distance along the Z axis */
/* The returned value is the computed covariance value */ /* 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*/ /*cubic covariance value for lag h*/
double cubic(double h, int cores); double cubic(double h);
/*truncation of the power spectrum to remove */ /*truncation of the power spectrum to remove */
/*high frequencies - isotropic case */ /*high frequencies - isotropic case */
@ -404,7 +402,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(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki);
/*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);
@ -413,7 +411,7 @@ double funtrun1(double x);
float G(float x); float G(float x);
/*gamma covariance value for lag h and exponent alpha*/ /*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))*/ /*returns the value ln(G(x))*/
float gammln(float xx); float gammln(float xx);
@ -425,7 +423,7 @@ float gammp(float a, float x);
/*and unit variance, using ran1(idum) as the source */ /*and unit variance, using ran1(idum) as the source */
/*of uniform deviates */ /*of uniform deviates */
/*idum: seed */ /*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*/ /*gaussian covariance value for lag h*/
double gaussian(double h); double gaussian(double h);
@ -480,7 +478,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 */ /*n: vector with the number of cells along the */
/* X, Y and Z axes for the underlying grid */ /* X, Y and Z axes for the underlying grid */
/* i = [0 1 2] */ /* 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*/ /*incomplete gamma function evaluated by its series*/
/*representation as gamser, also returns ln(G(a)) */ /*representation as gamser, also returns ln(G(a)) */
@ -505,7 +503,7 @@ void krig_stat(float* b, int n, struct vario_mod variogram, struct welldata_mod
/*i: considered direction */ /*i: considered direction */
/*scf: correlation length */ /*scf: correlation length */
/*ap: normalized anisotropy axes */ /*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/ /*calculates L.Z/
/* L : lower triangular matrix recorded */ /* L : lower triangular matrix recorded */
@ -519,7 +517,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); void LtimeZ(double* L, float* Z, float* b, int n);
/*determines the greatest prime factor of an integer*/ /*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 */ /*metrop returns a boolean varible that issues a */
/*verdict on whether to accept a reconfiguration */ /*verdict on whether to accept a reconfiguration */
@ -545,7 +543,7 @@ double power(double h, double alpha);
/*generates uniform deviates between 0 and 1*/ /*generates uniform deviates between 0 and 1*/
/*idum: seed */ /*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 */ /*calculates bt.b */
/* b : vector, bi, i = [0...n-1] */ /* b : vector, bi, i = [0...n-1] */

@ -1,25 +1,11 @@
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <time.h> #include <time.h>
/* compute the length for one dimension*/ /* compute the length for one dimension*/
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) {
double* used_ram_t0 = malloc(sizeof(double)); int maxfactor(int n);
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, N = %d, i = %d, D = %f, Nvari = %d", N, i, D, Nvari);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
int maxfactor(int n, int cores);
double temp1, temp2; double temp1, temp2;
int n, j, k, nmax; int n, j, k, nmax;
int nlimit = 13; int nlimit = 13;
@ -45,31 +31,11 @@ int length(int N, int i, double* scf, double* ap, double D, int Nvari, int cores
} }
if ((n % 2) != 0) if ((n % 2) != 0)
n = n + 1; n = n + 1;
nmax = maxfactor(n, cores); nmax = maxfactor(n);
while (nmax > nlimit) { while (nmax > nlimit) {
n += 2; n += 2;
nmax = maxfactor(n, cores); nmax = maxfactor(n);
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, n = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", n, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return n; return n;
} }

@ -1,177 +0,0 @@
/*
* Copyright (c) 2020 rxi
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal in the Software without restriction, including without limitation the
* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
* sell copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
* IN THE SOFTWARE.
*/
#include <stdlib.h>
#include "log.h"
#include <string.h>
#define MAX_CALLBACKS 32
typedef struct {
log_LogFn fn;
void *udata;
int level;
} Callback;
static struct {
void *udata;
log_LockFn lock;
int level;
bool quiet;
Callback callbacks[MAX_CALLBACKS];
} L;
static const char *level_strings[] = {
"TRACE", "DEBUG", "INFO", "WARN", "ERROR", "FATAL"
};
#ifdef LOG_USE_COLOR
static const char *level_colors[] = {
"\x1b[94m", "\x1b[36m", "\x1b[32m", "\x1b[33m", "\x1b[31m", "\x1b[35m"
};
#endif
static void stdout_callback(log_Event *ev) {
char buf[16];
buf[strftime(buf, sizeof(buf), "%H:%M:%S", ev->time)] = '\0';
#ifdef LOG_USE_COLOR
fprintf(
ev->udata, "%s %s%-5s\x1b[0m \x1b[0m%s:%d:\x1b[0m ",
buf, level_colors[ev->level], level_strings[ev->level],
ev->file, ev->line);
#else
fprintf(
ev->udata, "%s %-5s %s:%d: ",
buf, level_strings[ev->level], ev->file, ev->line);
#endif
vfprintf(ev->udata, ev->fmt, ev->ap);
fprintf(ev->udata, "\n");
fflush(ev->udata);
}
static void file_callback(log_Event *ev) {
char buf[64];
buf[strftime(buf, sizeof(buf), "%Y-%m-%d %H:%M:%S", ev->time)] = '\0';
fprintf(
ev->udata, "%s %-5s %s:%d: ",
buf, level_strings[ev->level], ev->file, ev->line);
vfprintf(ev->udata, ev->fmt, ev->ap);
fprintf(ev->udata, "\n");
fflush(ev->udata);
}
static void lock(void) {
if (L.lock) { L.lock(true, L.udata); }
}
static void unlock(void) {
if (L.lock) { L.lock(false, L.udata); }
}
const char* log_level_string(int level) {
return level_strings[level];
}
void log_set_lock(log_LockFn fn, void *udata) {
L.lock = fn;
L.udata = udata;
}
void log_set_level(int level) {
L.level = level;
}
void log_set_quiet(bool enable) {
L.quiet = enable;
}
int log_add_callback(log_LogFn fn, void *udata, int level) {
for (int i = 0; i < MAX_CALLBACKS; i++) {
if (!L.callbacks[i].fn) {
L.callbacks[i] = (Callback) { fn, udata, level };
return 0;
}
}
return -1;
}
int log_add_fp(FILE *fp, int level) {
return log_add_callback(file_callback, fp, level);
}
static void init_event(log_Event *ev, void *udata) {
if (!ev->time) {
time_t t = time(NULL);
ev->time = localtime(&t);
}
ev->udata = udata;
}
void log_log(int level, const char *file, int line, const char *fmt, ...) {
log_Event ev = {
.fmt = fmt,
.file = file,
.line = line,
.level = level,
};
char* env_var = getenv("ENV");
if (env_var != NULL && strcmp("false", env_var) == 0) return;
char* substr_mem = strstr(fmt, "MEM");
char* substr_cpu = strstr(fmt, "CPU");
if (env_var != NULL && strcmp("analysis", env_var) == 0 && substr_mem == NULL && substr_cpu == NULL) return;
lock();
if (!L.quiet && level >= L.level) {
init_event(&ev, stderr);
va_start(ev.ap, fmt);
stdout_callback(&ev);
va_end(ev.ap);
}
for (int i = 0; i < MAX_CALLBACKS && L.callbacks[i].fn; i++) {
Callback *cb = &L.callbacks[i];
if (level >= cb->level) {
init_event(&ev, cb->udata);
va_start(ev.ap, fmt);
cb->fn(&ev);
va_end(ev.ap);
}
}
unlock();
}

@ -1,49 +0,0 @@
/**
* Copyright (c) 2020 rxi
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the MIT license. See `log.c` for details.
*/
#ifndef LOG_H
#define LOG_H
#include <stdio.h>
#include <stdarg.h>
#include <stdbool.h>
#include <time.h>
#define LOG_VERSION "0.1.0"
typedef struct {
va_list ap;
const char *fmt;
const char *file;
struct tm *time;
void *udata;
int line;
int level;
} log_Event;
typedef void (*log_LogFn)(log_Event *ev);
typedef void (*log_LockFn)(bool lock, void *udata);
enum { LOG_TRACE, LOG_DEBUG, LOG_INFO, LOG_WARN, LOG_ERROR, LOG_FATAL };
#define log_trace(...) log_log(LOG_TRACE, __FILE__, __LINE__, __VA_ARGS__)
#define log_debug(...) log_log(LOG_DEBUG, __FILE__, __LINE__, __VA_ARGS__)
#define log_info(...) log_log(LOG_INFO, __FILE__, __LINE__, __VA_ARGS__)
#define log_warn(...) log_log(LOG_WARN, __FILE__, __LINE__, __VA_ARGS__)
#define log_error(...) log_log(LOG_ERROR, __FILE__, __LINE__, __VA_ARGS__)
#define log_fatal(...) log_log(LOG_FATAL, __FILE__, __LINE__, __VA_ARGS__)
const char* log_level_string(int level);
void log_set_lock(log_LockFn fn, void *udata);
void log_set_level(int level);
void log_set_quiet(bool enable);
int log_add_callback(log_LogFn fn, void *udata, int level);
int log_add_fp(FILE *fp, int level);
void log_log(int level, const char *file, int line, const char *fmt, ...);
#endif

@ -1,23 +1,9 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
/*determines the greatest prime factor of an integer*/
int maxfactor(int n, int cores) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress, n = %d", n);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
/*determines the greatest prime factor of an integer*/
int maxfactor(int n) {
int test_fact(int* pnum, int fact, int* pmaxfac); int test_fact(int* pnum, int fact, int* pmaxfac);
int lnum, fact; int lnum, fact;
int maxfac; int maxfac;
@ -46,25 +32,5 @@ int maxfactor(int n, int cores) {
maxfac = lnum; maxfac = lnum;
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, maxfac = %d, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", maxfac, time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
return maxfac; return maxfac;
} }

@ -1,94 +0,0 @@
#include <sys/types.h>
#include <sys/sysinfo.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/times.h>
#include <sys/vtimes.h>
#include <unistd.h>
#include "log.h"
#include "memory.h"
static unsigned long long lastTotalUser, lastTotalUserLow, lastTotalSys, lastTotalIdle;
static clock_t lastCPU, lastSysCPU, lastUserCPU;
static int numProcessors;
void getTotalVirtualMem(double* total_ram) {
const double megabyte = 1024 * 1024;
struct sysinfo si;
sysinfo(&si);
*total_ram = si.totalram / megabyte;
}
void getVirtualMemUsed(double* used_ram) {
const double megabyte = 1024 * 1024;
struct sysinfo si;
sysinfo(&si);
*used_ram = (si.totalram - si.freeram) / megabyte;
}
int parseLine(char* line) {
// This assumes that a digit will be found and the line ends in " Kb".
int i = strlen(line);
const char* p = line;
while (*p <'0' || *p > '9') p++;
line[i-3] = '\0';
i = atoi(p);
return i;
}
int getVirtualMemUsedByCurrentProcess() {
FILE* file = fopen("/proc/self/status", "r");
int result = -1;
char line[128];
while (fgets(line, 128, file) != NULL){
if (strncmp(line, "VmSize:", 7) == 0){
result = parseLine(line);
break;
}
}
fclose(file);
return result / 1024;
}
void skip_lines(FILE *fp, int numlines) {
int cnt = 0;
char ch;
while ((cnt < numlines) && ((ch = getc(fp)) != EOF)) {
if (ch == '\n') cnt++;
}
}
void get_stats(struct cpustat *st, int cpunum) {
FILE *fp = fopen("/proc/stat", "r");
int lskip = cpunum+1;
skip_lines(fp, lskip);
char cpun[255];
fscanf(fp, "%s %d %d %d %d %d %d %d", cpun, &(st->t_user), &(st->t_nice),
&(st->t_system), &(st->t_idle), &(st->t_iowait), &(st->t_irq),
&(st->t_softirq));
fclose(fp);
return;
}
double calculate_load(struct cpustat *prev, struct cpustat *cur) {
int idle_prev = (prev->t_idle) + (prev->t_iowait);
int idle_cur = (cur->t_idle) + (cur->t_iowait);
int nidle_prev = (prev->t_user) + (prev->t_nice) + (prev->t_system) + (prev->t_irq) + (prev->t_softirq);
int nidle_cur = (cur->t_user) + (cur->t_nice) + (cur->t_system) + (cur->t_irq) + (cur->t_softirq);
int total_prev = idle_prev + nidle_prev;
int total_cur = idle_cur + nidle_cur;
double totald = (double) total_cur - (double) total_prev;
double idled = (double) idle_cur - (double) idle_prev;
if (totald == 0 && idled == 0) return 0;
double cpu_perc = (1000 * (totald - idled) / totald + 1) / 10;
return cpu_perc;
}

@ -1,24 +0,0 @@
#include "sys/types.h"
#include "sys/sysinfo.h"
#include "stdlib.h"
#include "stdio.h"
#include "string.h"
#include "sys/times.h"
#include "sys/vtimes.h"
void getTotalVirtualMem();
void getVirtualMemUsed();
int getVirtualMemUsedByCurrentProcess();
struct cpustat {
unsigned long t_user;
unsigned long t_nice;
unsigned long t_system;
unsigned long t_idle;
unsigned long t_iowait;
unsigned long t_irq;
unsigned long t_softirq;
};
void get_stats(struct cpustat *st, int cpunum);
double calculate_load(struct cpustat *prev, struct cpustat *cur);

@ -1,5 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
@ -15,12 +15,6 @@
/* lognormal numbers */ /* lognormal numbers */
void nor2log(struct realization_mod* realin, int typelog, struct realization_mod* realout) { void nor2log(struct realization_mod* realin, int typelog, struct realization_mod* realout) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
log_info("RESULT = in progress");
int i; int i;
double coeff; double coeff;
@ -69,20 +63,8 @@ void nor2log(struct realization_mod* realin, int typelog, struct realization_mod
(*realout).vector[i] = exp((*realin).vector[i] * coeff); (*realout).vector[i] = exp((*realin).vector[i] * coeff);
break; break;
default: default:
log_error("RESULT = failed - Unexpected case in nor2log");
return; return;
break; break;
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,12 +1,10 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
/*nugget covariance function*/ /*nugget covariance function*/
double nugget(double h) { double nugget(double h) {
log_info("RESULT = in progress, h = %f", h);
if (h == 0) { if (h == 0) {
return (1.); return (1.);
} else { } else {

@ -1,10 +1,9 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
/*power covariance function*/ /*power covariance function*/
double power(double h, double alpha) { double power(double h, double alpha) {
log_info("RESULT = in progress, h = %f, alpha = %f", h, alpha);
return pow(h, alpha); return pow(h, alpha);
} }

@ -1,6 +1,6 @@
#include "geostat.h" #include "geostat.h"
#include "log.h"
#include "memory.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -22,12 +22,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, double* realization, int solver, int cores, long* seed) { void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, long* seed) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
int ntot; int ntot;
@ -38,25 +33,12 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
if (*seed > 0.0) if (*seed > 0.0)
*seed = -(*seed); *seed = -(*seed);
log_info("RESULT = in progress, n[0] = %d, n[1] = %d, n[2] = %d, solver = %d", n[0], n[1], n[2], solver);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
ntot = n[0] * n[1] * n[2]; ntot = n[0] * n[1] * n[2];
realization[0] = 0.; realization[0] = 0.;
/*printf("Antes de llamar a chunkarray read\n");
chunk_array_read((*realin).vector_2);
printf("Despues de llamar a chunkarray read\n");*/
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);
realization[i+1] = value; realization[i+1] = value;
//chunk_array_get((*realin).vector_2, i, &realization[i + 1]);
} }
} else { } else {
for (k = 1; k <= n[2]; k++) { for (k = 1; k <= n[2]; k++) {
@ -64,11 +46,8 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
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];
if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) { if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) {
//maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; double value = gasdev(seed, &idum2, &iy, iv);
//printf("Maille0 es %d", maille0);
double value = gasdev(seed, &idum2, &iy, iv, cores);
realization[maille1] = value; realization[maille1] = value;
//chunk_array_get((*realin).vector_2, maille0, &realization[maille1]);
} else { } else {
realization[maille1] = 0.; realization[maille1] = 0.;
} }
@ -76,23 +55,4 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
} }
} }
} }
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
log_info("RESULT = success, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", time_taken, *used_ram_tf - *used_ram_t0);
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,7 +1,7 @@
#include <time.h> #include <time.h>
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include "memory.h"
#define IM1 2147483563 #define IM1 2147483563
#define IM2 2147483399 #define IM2 2147483399
@ -18,25 +18,11 @@
#define EPS 1.2e-7 #define EPS 1.2e-7
#define RNMX (1.0 - EPS) #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]) {
double* used_ram_t0 = malloc(sizeof(double));
getVirtualMemUsed(used_ram_t0);
clock_t t = clock();
int j; int j;
long k; long k;
double temp; double temp;
log_info("RESULT = in progress, idum = %f, idum2 = %f, iy = %f", *idum, *idum2, *iy);
struct cpustat initial[cores];
struct cpustat final[cores];
for (int i = 0; i < cores; i++) {
get_stats(&initial[i], i - 1);
}
if (*idum <= 0) { if (*idum <= 0) {
if (-(*idum) < 1) if (-(*idum) < 1)
*idum = 1; *idum = 1;
@ -68,29 +54,10 @@ double ran2(long* idum, long* idum2, long* iy, long iv[NTAB], int cores) {
if (*iy < 1) if (*iy < 1)
(*iy) += IMM1; (*iy) += IMM1;
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time
for (int i = 0; i < cores; i++) {
get_stats(&final[i], i - 1);
}
for (int i = 0; i < cores; i++) {
log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
}
double* used_ram_tf = malloc(sizeof(double));
getVirtualMemUsed(used_ram_tf);
if ((temp = AM * (*iy)) > RNMX) { if ((temp = AM * (*iy)) > RNMX) {
log_info("RESULT = success, RNMX = %f, ELAPSED = %f, DIF USED VIRTUAL MEM = %5.1f MB", RNMX, time_taken, *used_ram_tf - *used_ram_t0);
return (RNMX); return (RNMX);
} }
else { else {
log_info("RESULT = success, temp = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", temp, time_taken, *used_ram_tf - *used_ram_t0);
return temp; return temp;
} }
free(used_ram_t0);
free(used_ram_tf);
} }

@ -1,20 +1,15 @@
#include "genlib.h" #include "genlib.h"
#include "log.h"
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
/*spherical covariance function*/ /*spherical covariance function*/
double spherical(double h) { double spherical(double h) {
double z; double z;
log_info("RESULT = in progress, h = %f", h);
if (h >= 1.) { if (h >= 1.) {
z = 0.; z = 0.;
} else { } else {
z = 1. - 1.5 * (double)h + 0.5 * (double)(h * h * h); z = 1. - 1.5 * (double)h + 0.5 * (double)(h * h * h);
} }
log_info("RESULT = success, z = %f", z);
return z; return z;
} }

@ -7,8 +7,6 @@
#include "toolsFFTMA.h" #include "toolsFFTMA.h"
#include "toolsFFTPSIM.h" #include "toolsFFTPSIM.h"
#include "toolsIO.h" #include "toolsIO.h"
#include "lib_src/log.h"
#include "chunk_array.h"
#include <Python.h> #include <Python.h>
#include <math.h> #include <math.h>
#include <numpy/arrayobject.h> #include <numpy/arrayobject.h>
@ -23,7 +21,6 @@
/* Y is the realization with mean and variance wanted */ /* Y is the realization with mean and variance wanted */
static PyObject* genFunc(PyObject* self, PyObject* args) { static PyObject* genFunc(PyObject* self, PyObject* args) {
log_info("RESULT = in progress");
int n[3]; int n[3];
struct realization_mod Z, Y; struct realization_mod Z, Y;
struct grid_mod grid; struct grid_mod grid;
@ -32,13 +29,11 @@ static PyObject* genFunc(PyObject* self, PyObject* args) {
struct statistic_mod stat; struct statistic_mod stat;
PyObject* out_array; PyObject* out_array;
npy_intp out_dims[NPY_MAXDIMS]; 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; return NULL;
printf("deberia imprimir\n"); Py_kgeneration(seed, grid, stat, variogram, &Z, &Y, n);
Py_kgeneration(seed, grid, stat, variogram, &Z, &Y, n, cores);
out_dims[0] = grid.NZ; out_dims[0] = grid.NZ;
out_dims[1] = grid.NY; out_dims[1] = grid.NY;
@ -50,34 +45,13 @@ static PyObject* genFunc(PyObject* self, PyObject* args) {
PyArray_ENABLEFLAGS(out_array, NPY_ARRAY_OWNDATA); PyArray_ENABLEFLAGS(out_array, NPY_ARRAY_OWNDATA);
printf("me localizo\n");
free(stat.mean); free(stat.mean);
free(stat.variance); free(stat.variance);
printf("termino stat\n");
free(variogram.var); free(variogram.var);
printf("1\n");
printf("2\n");
free(variogram.alpha); free(variogram.alpha);
printf("3\n");
free(variogram.scf); free(variogram.scf);
printf("4\n");
free(variogram.ap); free(variogram.ap);
//free(variogram.vario);
printf("Termino variogram\n");
printf("aca no deberia llegar\n");
/*chunk_array_free(Z.vector_2);
remove("realization1.txt");*/
log_info("RESULT = success");
return out_array; return out_array;
} }

@ -42,9 +42,6 @@ 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/log.c",
"./lib_src/memory.c",
"./lib_src/chunk_array.c"
], ],
) )

Binary file not shown.

@ -0,0 +1,54 @@
from FFTMA import gen
import numpy as np
import unittest
import time
def generate(N):
start_time = time.time()
nx, ny, nz = N,N,N
dx, dy, dz = 1.0, 1.0, 1.0
seed= 1548762 #rdi(10000,99999)
var=1
vario=2
alpha=1
lcx=2
lcy=4
lcz=16
ap1x=1
ap1y=0
ap1z=0
ap2x=0
ap2y=1
ap2z=0
v1 = (var, vario, alpha, lcx, lcy, lcz, ap1x, ap1y, ap1z, ap2x, ap2y, ap2z)
variograms = [v1]
mean=15.3245987
variance=3.5682389
typ=3
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
def test(N):
k = generate(N)
k_correct = np.load(f"out_{N}.npy")
comparison = k == k_correct
return comparison.all()
class TestGeneration(unittest.TestCase):
def test_8(self):
self.assertTrue(test(8))
def test_16(self):
self.assertTrue(test(16))
def test_32(self):
self.assertTrue(test(32))
def test_64(self):
self.assertTrue(test(64))
Loading…
Cancel
Save