Commit 0f2ce62b authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapbic: replaced critical section with store_link() by store_weightlinks4()

parent fd1ef9e6
......@@ -648,7 +648,7 @@ void remap_set_frac_min(int gridsize, double *array, double missval, remapgrid_t
int timer_remap, timer_remap_init, timer_remap_sort;
int timer_remap_bil, timer_remap_nn, timer_remap_con, timer_remap_con_l1, timer_remap_con_l2;
int timer_remap_bil, timer_remap_bic, timer_remap_dis, timer_remap_con, timer_remap_con_l1, timer_remap_con_l2;
static
void init_remap_timer(void)
......@@ -657,7 +657,8 @@ void init_remap_timer(void)
timer_remap_init = timer_new("remap init");
timer_remap_sort = timer_new("remap sort");
timer_remap_bil = timer_new("remap bil");
timer_remap_nn = timer_new("remap nn");
timer_remap_bic = timer_new("remap bic");
timer_remap_dis = timer_new("remap dis");
timer_remap_con = timer_new("remap con");
timer_remap_con_l1 = timer_new("remap con loop1");
timer_remap_con_l2 = timer_new("remap con loop2");
......
......@@ -170,37 +170,41 @@ void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
{
/* Local variables */
int search_result;
long tgt_grid_size;
long dst_add; /* destination addresss */
long tgt_cell_add; /* destination addresss */
int src_add[4]; /* address for the four source points */
double src_lats[4]; /* latitudes of four bilinear corners */
double src_lons[4]; /* longitudes of four bilinear corners */
double wgts[4][4]; /* bicubic weights for four corners */
double plat, plon; /* lat/lon coords of destination point */
double findex = 0;
extern int timer_remap_bic;
int remap_grid_type = src_grid->remap_grid_type;
if ( cdoVerbose ) cdoPrint("Called %s()", __func__);
progressInit();
if ( cdoTimer ) timer_start(timer_remap_bic);
tgt_grid_size = tgt_grid->size;
progressInit();
/* Compute mappings from source to target grid */
if ( src_grid->rank != 2 )
cdoAbort("Can not do bicubic interpolation when source grid rank != 2");
long tgt_grid_size = tgt_grid->size;
weightlinks4_t *weightlinks = (weightlinks4_t *) malloc(tgt_grid_size*sizeof(weightlinks4_t));
/* Loop over destination grid */
double findex = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
private(dst_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
shared(ompNumThreads, cdoVerbose, weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(dynamic,1)
#endif
/* grid_loop1 */
for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
{
int lprogress = 1;
if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;
......@@ -211,10 +215,12 @@ void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
findex++;
if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
if ( ! tgt_grid->mask[dst_add] ) continue;
weightlinks[tgt_cell_add].nlinks = 0;
if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
plat = tgt_grid->cell_center_lat[dst_add];
plon = tgt_grid->cell_center_lon[dst_add];
plat = tgt_grid->cell_center_lat[tgt_cell_add];
plon = tgt_grid->cell_center_lon[tgt_cell_add];
/* Find nearest square of grid points on source grid */
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
......@@ -239,19 +245,14 @@ void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
{
double iw, jw; /* current guess for bilinear coordinate */
tgt_grid->cell_frac[dst_add] = 1.;
tgt_grid->cell_frac[tgt_cell_add] = 1.;
if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
{
/* Successfully found iw,jw - compute weights */
set_bicubic_weights(iw, jw, wgts);
sort_bicubic_adds(src_add, wgts);
#if defined(_OPENMP)
#pragma omp critical
#endif
store_link_bicub(rv, dst_add, src_add, wgts);
store_weightlinks4(4, src_add, wgts, tgt_cell_add, weightlinks);
}
else
{
......@@ -271,17 +272,18 @@ void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
{
renormalize_weights(src_lats, wgts);
tgt_grid->cell_frac[dst_add] = 1.;
tgt_grid->cell_frac[tgt_cell_add] = 1.;
sort_bicubic_adds(src_add, wgts);
#if defined(_OPENMP)
#pragma omp critical
#endif
store_link_bicub(rv, dst_add, src_add, wgts);
store_weightlinks4(4, src_add, wgts, tgt_cell_add, weightlinks);
}
}
} /* grid_loop1 */
}
if ( cdoTimer ) timer_stop(timer_remap_bic);
weightlinks2remaplinks4(tgt_grid_size, weightlinks, rv);
if ( weightlinks ) free(weightlinks);
} /* scrip_remap_weights_bicubic */
......@@ -296,44 +298,42 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
{
/* Local variables */
int search_result;
long tgt_grid_size;
long dst_add; /* destination addresss */
long tgt_cell_add; /* destination addresss */
int src_add[4]; /* address for the four source points */
double src_lats[4]; /* latitudes of four bilinear corners */
double src_lons[4]; /* longitudes of four bilinear corners */
double wgts[4][4]; /* bicubic weights for four corners */
double plat, plon; /* lat/lon coords of destination point */
double findex = 0;
double *grad1_lat, *grad1_lon, *grad1_latlon;
int remap_grid_type = src_grid->remap_grid_type;
if ( cdoVerbose ) cdoPrint("Called %s()", __func__);
progressInit();
tgt_grid_size = tgt_grid->size;
long tgt_grid_size = tgt_grid->size;
/* Compute mappings from source to target grid */
if ( src_grid->rank != 2 )
cdoAbort("Can not do bicubic interpolation when source grid rank != 2");
grad1_lat = (double*) malloc(src_grid->size*sizeof(double));
grad1_lon = (double*) malloc(src_grid->size*sizeof(double));
grad1_latlon = (double*) malloc(src_grid->size*sizeof(double));
double *grad1_lat = (double*) malloc(src_grid->size*sizeof(double));
double *grad1_lon = (double*) malloc(src_grid->size*sizeof(double));
double *grad1_latlon = (double*) malloc(src_grid->size*sizeof(double));
remap_gradients(*src_grid, src_array, grad1_lat, grad1_lon, grad1_latlon);
/* Loop over destination grid */
double findex = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, grad1_lat, grad1_lon, grad1_latlon, findex) \
private(dst_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
shared(ompNumThreads, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, grad1_lat, grad1_lon, grad1_latlon, findex) \
private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(dynamic,1)
#endif
/* grid_loop1 */
for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
{
int lprogress = 1;
if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;
......@@ -344,12 +344,12 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
findex++;
if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
tgt_array[dst_add] = missval;
tgt_array[tgt_cell_add] = missval;
if ( ! tgt_grid->mask[dst_add] ) continue;
if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
plat = tgt_grid->cell_center_lat[dst_add];
plon = tgt_grid->cell_center_lon[dst_add];
plat = tgt_grid->cell_center_lat[tgt_cell_add];
plon = tgt_grid->cell_center_lon[tgt_cell_add];
/* Find nearest square of grid points on source grid */
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
......@@ -374,7 +374,7 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
{
double iw, jw; /* current guess for bilinear coordinate */
tgt_grid->cell_frac[dst_add] = 1.;
tgt_grid->cell_frac[tgt_cell_add] = 1.;
if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
{
......@@ -383,7 +383,7 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
sort_bicubic_adds(src_add, wgts);
bicubic_remap(&tgt_array[dst_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
bicubic_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
}
else
{
......@@ -403,14 +403,14 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
{
renormalize_weights(src_lats, wgts);
tgt_grid->cell_frac[dst_add] = 1.;
tgt_grid->cell_frac[tgt_cell_add] = 1.;
sort_bicubic_adds(src_add, wgts);
bicubic_remap(&tgt_array[dst_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
bicubic_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
}
}
} /* grid_loop1 */
}
free(grad1_lat);
free(grad1_lon);
......
......@@ -194,7 +194,7 @@ void scrip_remap_weights_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid,
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
shared(ompNumThreads, cdoVerbose, weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(static)
#endif
......@@ -298,7 +298,6 @@ void scrip_remap_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid, const do
double src_lons[4]; /* longitudes of four bilinear corners */
double wgts[4]; /* bilinear weights for four corners */
double plat, plon; /* lat/lon coords of destination point */
double findex = 0;
extern int timer_remap_bil;
int remap_grid_type = src_grid->remap_grid_type;
......@@ -315,15 +314,16 @@ void scrip_remap_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid, const do
if ( src_grid->rank != 2 )
cdoAbort("Can not do bilinear interpolation when source grid rank != 2");
double findex = 0;
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex) \
shared(ompNumThreads, cdoVerbose, cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex) \
private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(static)
#endif
/* grid_loop1 */
for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
{
int lprogress = 1;
......@@ -404,7 +404,7 @@ void scrip_remap_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid, const do
bilinear_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add);
}
}
} /* grid_loop1 */
}
if ( cdoTimer ) timer_stop(timer_remap_bil);
} /* scrip_remap_bilinear */
......@@ -709,7 +709,7 @@ void remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapva
#if defined(_OPENMP)
#pragma omp parallel for default(shared) \
shared(ompNumThreads, cdoTimer, lyac, nbins, num_wts, src_remap_grid_type, tgt_remap_grid_type, src_grid_bound_box, \
shared(ompNumThreads, lyac, nbins, num_wts, src_remap_grid_type, tgt_remap_grid_type, src_grid_bound_box, \
src_edge_type, tgt_edge_type, partial_areas2, partial_weights2, \
rv, cdoVerbose, max_srch_cells2, tgt_num_cell_corners, target_cell_type, \
weightlinks, \
......
......@@ -1380,7 +1380,7 @@ void scrip_remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
#if defined(_OPENMP)
#pragma omp parallel for default(shared) \
shared(ompNumThreads, cdoTimer, nbins, num_wts, src_centroid_lon, src_centroid_lat, \
shared(ompNumThreads, nbins, num_wts, src_centroid_lon, src_centroid_lat, \
remap_store_link_fast, grid_store, link_add1, link_add2, rv, cdoVerbose, max_subseg, \
srch_corner_lat2, srch_corner_lon2, max_srch_cells2, \
src_num_cell_corners, srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, srch_add2, findex) \
......@@ -1605,7 +1605,7 @@ void scrip_remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
#if defined(_OPENMP)
#pragma omp parallel for default(shared) \
shared(ompNumThreads, cdoTimer, nbins, num_wts, tgt_centroid_lon, tgt_centroid_lat, \
shared(ompNumThreads, nbins, num_wts, tgt_centroid_lon, tgt_centroid_lat, \
remap_store_link_fast, grid_store, link_add1, link_add2, rv, cdoVerbose, max_subseg, \
srch_corner_lat2, srch_corner_lon2, max_srch_cells2, \
tgt_num_cell_corners, srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, srch_add2, findex) \
......
......@@ -414,7 +414,7 @@ void scrip_remap_weights_distwgt(int num_neighbors, remapgrid_t *src_grid, remap
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, weightlinks, num_neighbors, remap_grid_type, src_grid, tgt_grid, rv, tgt_grid_size, coslat, coslon, sinlat, sinlon, findex) \
shared(ompNumThreads, weightlinks, num_neighbors, remap_grid_type, src_grid, tgt_grid, rv, tgt_grid_size, coslat, coslon, sinlat, sinlon, findex) \
private(tgt_cell_add, n, nadds, dist_tot, plat, plon)
#endif
for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
......
......@@ -94,10 +94,8 @@ void sort_add_and_wgts4(size_t num_weights, int *src_add, double *wgts[4])
for ( n = 0; n < num_weights; ++n )
{
addweights[n].add = src_add[n];
addweights[n].weight[0] = wgts[0][n];
addweights[n].weight[1] = wgts[1][n];
addweights[n].weight[2] = wgts[2][n];
addweights[n].weight[3] = wgts[3][n];
for ( long k = 0; k < 4; ++k )
addweights[n].weight[k] = wgts[k][n];
}
qsort(addweights, num_weights, sizeof(addweight4_t), cmp_adds4);
......@@ -105,10 +103,8 @@ void sort_add_and_wgts4(size_t num_weights, int *src_add, double *wgts[4])
for ( n = 0; n < num_weights; ++n )
{
src_add[n] = addweights[n].add;
wgts[0][n] = addweights[n].weight[0];
wgts[1][n] = addweights[n].weight[1];
wgts[2][n] = addweights[n].weight[2];
wgts[3][n] = addweights[n].weight[3];
for ( long k = 0; k < 4; ++k )
wgts[k][n] = addweights[n].weight[k];
}
}
......@@ -135,7 +131,7 @@ void store_weightlinks(long num_weights, int *srch_add, double *weights, long ce
}
void store_weightlinks4(long num_weights, int *srch_add, double *weights[4], long cell_add, weightlinks4_t *weightlinks)
void store_weightlinks4(long num_weights, int *srch_add, double weights[4][4], long cell_add, weightlinks4_t *weightlinks)
{
weightlinks[cell_add].nlinks = 0;
weightlinks[cell_add].offset = 0;
......@@ -146,10 +142,8 @@ void store_weightlinks4(long num_weights, int *srch_add, double *weights[4], lon
for ( long n = 0; n < num_weights; ++n )
{
addweights[n].add = srch_add[n];
addweights[n].weight[0] = weights[0][n];
addweights[n].weight[1] = weights[1][n];
addweights[n].weight[2] = weights[2][n];
addweights[n].weight[3] = weights[3][n];
for ( long k = 0; k < 4; ++k )
addweights[n].weight[k] = weights[k][n];
}
sort_addweights4(num_weights, addweights);
......@@ -160,7 +154,7 @@ void store_weightlinks4(long num_weights, int *srch_add, double *weights[4], lon
}
void weightlinks2remaplinks(long tgt_grid_size, weightlinks_t *weightlinks, remapvars_t *rv)
void weightlinks2remaplinks(long tgt_grid_size, weightlinks_t *weightlinks, remapvars_t *rv)
{
long tgt_cell_add;
long nlinks = 0;
......@@ -206,7 +200,7 @@ void weightlinks2remaplinks(long tgt_grid_size, weightlinks_t *weightlinks, rem
}
void weightlinks2remaplinks4(long tgt_grid_size, weightlinks4_t *weightlinks, remapvars_t *rv)
void weightlinks2remaplinks4(long tgt_grid_size, weightlinks4_t *weightlinks, remapvars_t *rv)
{
long tgt_cell_add;
long nlinks = 0;
......@@ -243,10 +237,8 @@ void weightlinks2remaplinks4(long tgt_grid_size, weightlinks4_t *weightlinks, r
{
rv->src_cell_add[offset+ilink] = weightlinks[tgt_cell_add].addweights[ilink].add;
rv->tgt_cell_add[offset+ilink] = tgt_cell_add;
rv->wts[(offset+ilink)*4+0] = weightlinks[tgt_cell_add].addweights[ilink].weight[0];
rv->wts[(offset+ilink)*4+1] = weightlinks[tgt_cell_add].addweights[ilink].weight[1];
rv->wts[(offset+ilink)*4+2] = weightlinks[tgt_cell_add].addweights[ilink].weight[2];
rv->wts[(offset+ilink)*4+3] = weightlinks[tgt_cell_add].addweights[ilink].weight[3];
for ( long k = 0; k < 4; ++k )
rv->wts[(offset+ilink)*4+k] = weightlinks[tgt_cell_add].addweights[ilink].weight[k];
}
free(weightlinks[tgt_cell_add].addweights);
}
......
......@@ -8,26 +8,29 @@ typedef struct
double weight;
} addweight_t;
typedef struct {
int nlinks;
int offset;
addweight_t *addweights;
} weightlinks_t;
typedef struct
{
int add;
double weight[4];
} addweight4_t;
typedef struct {
int nlinks;
int offset;
addweight_t *addweights;
} weightlinks_t;
typedef struct {
int nlinks;
int offset;
addweight4_t *addweights;
} weightlinks4_t;
void store_weightlinks(long num_weights, int *srch_add, double *weights, long cell_add, weightlinks_t *weightlinks);
void weightlinks2remaplinks(long tgt_grid_size, weightlinks_t *weightlinks, remapvars_t *rv);
void store_weightlinks4(long num_weights, int *srch_add, double weights[4][4], long cell_add, weightlinks4_t *weightlinks);
void weightlinks2remaplinks(long tgt_grid_size, weightlinks_t *weightlinks, remapvars_t *rv);
void weightlinks2remaplinks4(long tgt_grid_size, weightlinks4_t *weightlinks, remapvars_t *rv);
void sort_add_and_wgts(size_t num_weights, int *src_add, double *wgts);
void sort_add_and_wgts4(size_t num_weights, int *src_add, double *wgts[4]);
......
......@@ -1000,7 +1000,7 @@ void remap_vars_init(int map_type, long src_grid_size, long tgt_grid_size, remap
if ( map_type == MAP_TYPE_CONSERV ) rv->sort_add = TRUE;
else if ( map_type == MAP_TYPE_CONSERV_YAC ) rv->sort_add = FALSE;
else if ( map_type == MAP_TYPE_BILINEAR ) rv->sort_add = FALSE;
else if ( map_type == MAP_TYPE_BICUBIC ) rv->sort_add = TRUE;
else if ( map_type == MAP_TYPE_BICUBIC ) rv->sort_add = FALSE;
else if ( map_type == MAP_TYPE_DISTWGT ) rv->sort_add = FALSE;
else cdoAbort("Unknown mapping method!");
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment