Skip to content

Commit

Permalink
Revised OpenMP parallelization in cmultiscale code.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lars2015
lars2015 committed Sep 25, 2024
1 parent 9b62292 commit 5495bb9
Show file tree
Hide file tree
Showing 2 changed files with 110 additions and 97 deletions.
204 changes: 107 additions & 97 deletions src/mptrac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -505,112 +505,120 @@ void compress_cms(
else {

/* Init... */
double cr = 0;
cms_module_t *cms_ptr[EP];
cms_sol_t *cms_sol[EP];

/* Loop over levels... */
#pragma omp parallel for default(shared)
for (size_t ip = 0; ip < np; ip++) {

/* Allocate... */
float *tmp_arr;
ALLOC(tmp_arr, float,
nxy);

/* Copy level data... */
for (size_t ix = 0; ix < nx; ++ix)
for (size_t iy = 0; iy < ny; ++iy)
tmp_arr[ARRAY_2D(ix, iy, ny)] = array[ARRAY_3D(ix, iy, ny, ip, np)];

/* Initialize multiscale module... */
cms_ptr[ip] = cms_init(cms_param);

/* Create solution pointer... */
cms_sol[ip] = cms_read_arr(cms_ptr[ip], tmp_arr, lon, lat, nx, ny);

/* Set eps threshold value... */
if (strcasecmp(varname, "Z") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_z);
else if (strcasecmp(varname, "T") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_t);
else if (strcasecmp(varname, "U") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_u);
else if (strcasecmp(varname, "V") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_v);
else if (strcasecmp(varname, "W") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_w);
else if (strcasecmp(varname, "PV") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_pv);
else if (strcasecmp(varname, "H2O") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_h2o);
else if (strcasecmp(varname, "O3") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_o3);
else if (strcasecmp(varname, "LWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_lwc);
else if (strcasecmp(varname, "RWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_rwc);
else if (strcasecmp(varname, "IWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_iwc);
else if (strcasecmp(varname, "SWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_swc);
else if (strcasecmp(varname, "CC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_cc);
else
ERRMSG("Variable name unknown!");

/* Coarsening... */
cms_coarsening(cms_ptr[ip], cms_sol[ip],
(unsigned int) ctl->met_cms_heur);

/* Free... */
free(tmp_arr);
}
/* Loop over batches... */
const size_t dip = (ctl->met_cms_batch <= 0
? (size_t) omp_get_max_threads()
: (size_t) ctl->met_cms_batch);
for (size_t ip0 = 0; ip0 < np; ip0 += dip) {

/* Loop over levels... */
double cr = 0;
for (size_t ip = 0; ip < np; ip++) {

/* Allocate... */
double *tmp_cms, *tmp_org, *tmp_diff;
ALLOC(tmp_cms, double,
nxy);
ALLOC(tmp_org, double,
nxy);
ALLOC(tmp_diff, double,
nxy);

/* Evaluate... */
/* Loop over levels... */
#pragma omp parallel for default(shared)
for (size_t ix = 0; ix < nx; ix++)
for (size_t iy = 0; iy < ny; iy++) {
size_t idx = ARRAY_2D(ix, iy, ny);
double x[] = { lon[ix], lat[iy] };
cms_eval(cms_ptr[ip], cms_sol[ip], x, &tmp_cms[idx]);
tmp_org[idx] = array[ARRAY_3D(ix, iy, ny, ip, np)];
tmp_diff[idx] = tmp_cms[idx] - tmp_org[idx];
}
for (size_t ip = ip0; ip < MIN(ip0 + dip, np); ip++) {

/* Allocate... */
float *tmp_arr;
ALLOC(tmp_arr, float,
nxy);

/* Copy level data... */
for (size_t ix = 0; ix < nx; ++ix)
for (size_t iy = 0; iy < ny; ++iy)
tmp_arr[ARRAY_2D(ix, iy, ny)] =
array[ARRAY_3D(ix, iy, ny, ip, np)];

/* Initialize multiscale module... */
cms_ptr[ip] = cms_init(cms_param);

/* Create solution pointer... */
cms_sol[ip] = cms_read_arr(cms_ptr[ip], tmp_arr, lon, lat, nx, ny);

/* Set eps threshold value... */
if (strcasecmp(varname, "Z") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_z);
else if (strcasecmp(varname, "T") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_t);
else if (strcasecmp(varname, "U") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_u);
else if (strcasecmp(varname, "V") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_v);
else if (strcasecmp(varname, "W") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_w);
else if (strcasecmp(varname, "PV") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_pv);
else if (strcasecmp(varname, "H2O") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_h2o);
else if (strcasecmp(varname, "O3") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_o3);
else if (strcasecmp(varname, "LWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_lwc);
else if (strcasecmp(varname, "RWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_rwc);
else if (strcasecmp(varname, "IWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_iwc);
else if (strcasecmp(varname, "SWC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_swc);
else if (strcasecmp(varname, "CC") == 0)
cms_set_eps(cms_ptr[ip], ctl->met_cms_eps_cc);
else
ERRMSG("Variable name unknown!");

/* Write info... */
LOG(2,
"cmultiscale: var= %s / lev= %lu / ratio= %g / rho= %g"
" / mean= %g / sd= %g / min= %g / max= %g", varname, ip,
cms_compression_rate(cms_ptr[ip], cms_sol[ip]),
gsl_stats_correlation(tmp_cms, 1, tmp_org, 1, nxy),
gsl_stats_mean(tmp_diff, 1, nxy), gsl_stats_sd(tmp_diff, 1, nxy),
gsl_stats_min(tmp_diff, 1, nxy), gsl_stats_max(tmp_diff, 1, nxy));
/* Coarsening... */
cms_coarsening(cms_ptr[ip], cms_sol[ip],
(unsigned int) ctl->met_cms_heur);

/* Calculate mean compression ratio... */
cr += cms_compression_rate(cms_ptr[ip], cms_sol[ip]) / (double) np;
/* Free... */
free(tmp_arr);
}

/* Save binary data... */
cms_save_sol(cms_sol[ip], cms_ptr[ip], inout);
/* Loop over levels... */
for (size_t ip = ip0; ip < MIN(ip0 + dip, np); ip++) {

/* Allocate... */
double *tmp_cms, *tmp_org, *tmp_diff;
ALLOC(tmp_cms, double,
nxy);
ALLOC(tmp_org, double,
nxy);
ALLOC(tmp_diff, double,
nxy);

/* Evaluate... */
#pragma omp parallel for default(shared)
for (size_t ix = 0; ix < nx; ix++)
for (size_t iy = 0; iy < ny; iy++) {
size_t idx = ARRAY_2D(ix, iy, ny);
double x[] = { lon[ix], lat[iy] };
cms_eval(cms_ptr[ip], cms_sol[ip], x, &tmp_cms[idx]);
tmp_org[idx] = array[ARRAY_3D(ix, iy, ny, ip, np)];
tmp_diff[idx] = tmp_cms[idx] - tmp_org[idx];
}

/* Free... */
cms_delete_sol(cms_sol[ip]);
cms_delete_module(cms_ptr[ip]);
free(tmp_cms);
free(tmp_org);
free(tmp_diff);
/* Write info... */
LOG(2,
"cmultiscale: var= %s / lev= %lu / ratio= %g / rho= %g"
" / mean= %g / sd= %g / min= %g / max= %g", varname, ip,
cms_compression_rate(cms_ptr[ip], cms_sol[ip]),
gsl_stats_correlation(tmp_cms, 1, tmp_org, 1, nxy),
gsl_stats_mean(tmp_diff, 1, nxy), gsl_stats_sd(tmp_diff, 1, nxy),
gsl_stats_min(tmp_diff, 1, nxy), gsl_stats_max(tmp_diff, 1, nxy));

/* Calculate mean compression ratio... */
cr += cms_compression_rate(cms_ptr[ip], cms_sol[ip]) / (double) np;

/* Save binary data... */
cms_save_sol(cms_sol[ip], cms_ptr[ip], inout);

/* Free... */
cms_delete_sol(cms_sol[ip]);
cms_delete_module(cms_ptr[ip]);
free(tmp_cms);
free(tmp_org);
free(tmp_diff);
}
}

/* Write info... */
Expand Down Expand Up @@ -5031,6 +5039,8 @@ void read_ctl(
scan_ctl(filename, argc, argv, "MET_ZFP_TOL_T", -1, "5.0", NULL);
ctl->met_zfp_tol_z =
scan_ctl(filename, argc, argv, "MET_ZFP_TOL_Z", -1, "0.5", NULL);
ctl->met_cms_batch =
(int) scan_ctl(filename, argc, argv, "MET_CMS_BATCH", -1, "-1", NULL);
ctl->met_cms_heur =
(int) scan_ctl(filename, argc, argv, "MET_CMS_HEUR", -1, "1", NULL);
ctl->met_cms_eps_z =
Expand Down
3 changes: 3 additions & 0 deletions src/mptrac.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2476,6 +2476,9 @@ typedef struct {
/*! ZFP compression tolerance for geopotential height. */
double met_zfp_tol_z;

/*! cmultiscale batch size. */
int met_cms_batch;

/*! cmultiscale coarsening heuristics
(0=default, 1=mean diff, 2=median diff, 3=max diff). */
int met_cms_heur;
Expand Down

0 comments on commit 5495bb9

Please sign in to comment.