#! /bin/sh /usr/share/dpatch/dpatch-run ## omp-repair.dpatch by ## ## All lines beginning with `## DP:' are a description of the patch. ## DP: Specifies shared and private data in some "#pragma omp" statements that caused problems. ## DP: Also includes some minor improvements in function headers: const declaration for pointers. @DPATCH@ diff -urNad libgpiv-0.6.1~/include/gpiv/gpiv-piv_utils.h libgpiv-0.6.1/include/gpiv/gpiv-piv_utils.h --- libgpiv-0.6.1~/include/gpiv/gpiv-piv_utils.h 2009-10-09 10:11:19.000000000 +0200 +++ libgpiv-0.6.1/include/gpiv/gpiv-piv_utils.h 2009-10-09 10:11:29.000000000 +0200 @@ -112,7 +112,7 @@ * @return NULL on success or error message on failure */ gchar * -gpiv_0_pivdata (GpivPivData *piv_data +gpiv_0_pivdata (const GpivPivData *piv_data ); @@ -141,7 +141,7 @@ */ gchar * gpiv_ovwrt_pivdata (const GpivPivData *piv_data_in, - GpivPivData *piv_data_out + const GpivPivData *piv_data_out ); diff -urNad libgpiv-0.6.1~/lib/imgproc_deform.c libgpiv-0.6.1/lib/imgproc_deform.c --- libgpiv-0.6.1~/lib/imgproc_deform.c 2009-10-09 10:11:19.000000000 +0200 +++ libgpiv-0.6.1/lib/imgproc_deform.c 2009-10-09 10:11:29.000000000 +0200 @@ -185,9 +185,9 @@ * and second frame with half shift backwards. */ /* OK: results approved */ -#pragma omp parallel sections num_threads(2) +#pragma omp parallel sections { -#pragma omp section +#pragma omp section { if ((err_msg = imgproc_deform_frame_from_pivdata(image->frame1, image->header, @@ -200,7 +200,7 @@ #endif } } -#pragma omp section +#pragma omp section { if ((err_msg = imgproc_deform_frame_from_pivdata(image->frame2, image->header, @@ -454,7 +454,7 @@ /* BUGFIX SamplesToCoefficients: OMP causes error */ -/* #pragma omp parallel for */ +/* KEEP OUT! #pragma omp parallel for */ for (y = 0L; y < Height; y++) { GetRow(Image, y, Line, Width); ConvertToInterpolationCoefficients(Line, Width, Pole, NbPoles, @@ -473,7 +473,7 @@ /* BUGFIX SamplesToCoefficients: OMP causes error */ -/* #pragma omp parallel for */ +/* KEEP OUT! #pragma omp parallel for */ for (x = 0L; x < Width; x++) { GetColumn(Image, Width, x, Line, Height); ConvertToInterpolationCoefficients(Line, Height, Pole, NbPoles, DBL_EPSILON); @@ -712,7 +712,7 @@ pd = gpiv_alloc_pivdata (image_par->ncolumns, image_par->nrows); /* OK: results approved */ -#pragma omp parallel for +#pragma omp parallel for shared(pd) private(i, j) for (i = 0; i < pd->ny; i++) { for (j = 0; j < pd->nx; j++) { pd->point_x[i][j] = (float) j; @@ -745,7 +745,7 @@ #endif /* SMOOTH */ /* OK: results approved */ -#pragma omp parallel for +#pragma omp parallel for shared(pd) private(i, j) for (i = 0; i < pd->ny; i++) { for (j = 0; j < pd->nx; j++) { pd->dx[i][j] = magnitude * pd->dx[i][j]; @@ -811,7 +811,7 @@ p = ImageRasterArray; /* OK: results approved */ -#pragma omp parallel for +#pragma omp parallel for shared(image_par, img, p) private(i, j) for (i = 0; i < image_par->nrows; i++) { for (j = 0; j < image_par->ncolumns; j++) { p[i*image_par->ncolumns + j] = (float) img[i][j]; @@ -839,7 +839,7 @@ /* BUGFIX imgproc_deform_frame: OMP causes differences in results */ /* BUGFIX imgproc_deform_frame: When disabled OMP in imgproc_deform */ -/* #pragma omp parallel for private(x1, y1) */ +/* KEEP OUT! #pragma omp parallel for private(x1, y1) */ for (i = 0; i < image_par->nrows; i++) { for (j = 0; j < image_par->ncolumns; j++) { x1 = (double) j - (double) piv_data->dx[i][j]; @@ -875,7 +875,7 @@ p = OutputImage; /* OK: results approved */ -#pragma omp parallel for private(rounded) +#pragma omp parallel for private(i, j, rounded) for (i = 0; i < image_par->nrows; i++) { for (j = 0; j < image_par->ncolumns; j++) { rounded = (long)floor ((double) p[i*image_par->ncolumns + j] + 0.5); diff -urNad libgpiv-0.6.1~/lib/piv_utils.c libgpiv-0.6.1/lib/piv_utils.c --- libgpiv-0.6.1~/lib/piv_utils.c 2009-10-09 10:11:19.000000000 +0200 +++ libgpiv-0.6.1/lib/piv_utils.c 2009-10-09 10:11:29.000000000 +0200 @@ -115,7 +115,7 @@ /* * Initializing values for structure members */ -#pragma omp parallel for +#pragma omp parallel for shared(piv_data) private(i, j) for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { piv_data->point_x[i][j] = 0.0; @@ -253,7 +253,7 @@ return NULL; } -#pragma omp parallel for +#pragma omp parallel for shared(piv_data_out, piv_data) private(i, j) for (i = 0; i < piv_data_out->ny; i++) { for (j = 0; j < piv_data_out->nx; j++) { piv_data_out->point_x[i][j] = piv_data->point_x[i][j]; @@ -264,7 +264,8 @@ piv_data_out->peak_no[i][j] = piv_data->peak_no[i][j]; } } - + + #ifdef DEBUG2 gpiv_write_pivdata (NULL, piv_data_out, FALSE); @@ -276,7 +277,7 @@ gchar * gpiv_ovwrt_pivdata (const GpivPivData *piv_data_in, - GpivPivData *piv_data_out + const GpivPivData *piv_data_out ) /*----------------------------------------------------------------------------- * DESCRIPTION: @@ -321,7 +322,7 @@ /* * Overwriting */ -#pragma omp parallel for +#pragma omp parallel for shared(piv_data_out, piv_data_in) private(i, j) for (i = 0; i < piv_data_in->ny; i++) { for (j = 0; j < piv_data_in->nx; j++) { piv_data_out->point_x[i][j] = piv_data_in->point_x[i][j]; @@ -332,7 +333,6 @@ piv_data_out->peak_no[i][j] = piv_data_in->peak_no[i][j]; } } - #ifdef DEBUG2 gpiv_write_pivdata (NULL, piv_data_out, FALSE); @@ -343,7 +343,7 @@ gchar * -gpiv_0_pivdata (GpivPivData *piv_data +gpiv_0_pivdata (const GpivPivData *piv_data ) /*----------------------------------------------------------------------------- * DESCRIPTION: @@ -402,11 +402,9 @@ } -#pragma omp parallel for +#pragma omp parallel for shared(piv_data) private(i, j) for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { - /* piv_data->point_x[i][j] = 0.0; */ - /* piv_data->point_y[i][j] = 0.0; */ piv_data->dx[i][j] = 0.0; piv_data->dy[i][j] = 0.0; piv_data->snr[i][j] = 0.0; @@ -450,7 +448,7 @@ return err_msg; } -#pragma omp parallel for +#pragma omp parallel for shared(piv_data_out, piv_data_in) private(i, j) for (i = 0; i < piv_data_out->ny; i++) { for (j = 0; j < piv_data_out->nx; j++) { piv_data_out->point_x[i][j] = piv_data_in->point_x[i][j]; @@ -488,7 +486,6 @@ } -#pragma omp parallel for for (i = 0; i < gpd->ny; i++) { for (j = 0; j < gpd->nx; j++) { *sum += gpd->dx[i][j]; diff -urNad libgpiv-0.6.1~/lib/valid.c libgpiv-0.6.1/lib/valid.c --- libgpiv-0.6.1~/lib/valid.c 2009-10-09 10:11:19.000000000 +0200 +++ libgpiv-0.6.1/lib/valid.c 2009-10-09 10:12:22.000000000 +0200 @@ -124,9 +124,11 @@ float *dx_m, *dy_m, *dx_org, *dy_org; float *r_m; float residu; + /* * Obtain length of array */ +#pragma parallel for shared(data, incl_point, vlength) private(k, l) for (k = i - N; k <= i + N; k++) { if ((k >= 0) && (k < data->ny)) { for (l = j - N; l <= j + N; l++) { @@ -179,8 +181,18 @@ /* * Sorting dx_m and dy_m arrays and searching median index and value */ - qsort(dx_m, vlength, sizeof(float), compare_float); - qsort(dy_m, vlength, sizeof(float), compare_float); +#pragma omp sections + { +#pragma omp section + { + qsort(dx_m, vlength, sizeof(float), compare_float); + } +#pragma omp section + { + qsort(dy_m, vlength, sizeof(float), compare_float); + } + } + if (incl_point) { median_index = (int) (vlength - 1) / 2; } else { @@ -282,7 +294,6 @@ { gchar *err_msg = NULL; gint i, j; -/* guint i_mx, j_mx, i_my, j_my; */ /* out_data = gpiv_cp_pivdata (piv_data); */ @@ -292,7 +303,6 @@ } else if (valid_par->residu_type == GPIV_VALID_RESIDUTYPE__MEDIAN) { -/* #pragma omp parallel for */ for (i = 1; i < piv_data->ny - 1; i++) { for (j = 1; j < piv_data->nx - 1; j++) { guint i_mx, j_mx, i_my, j_my; @@ -309,8 +319,6 @@ } else if (valid_par->residu_type == GPIV_VALID_RESIDUTYPE__NORMMEDIAN) { gfloat residu_from_median_dxdy = 1.111, residu_norm = 5.555; -/* BUGFIX: gpiv_valid_residu: different results with OMP */ -/* #pragma omp parallel for private(residu_from_median_dxdy, residu_norm) */ for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { guint i_mx, j_mx, i_my, j_my; @@ -473,6 +481,10 @@ gfloat dx, dy, snr; FILE *fp; + guint count = 0; + gfloat dx_sum = 0, dy_sum = 0; + + /* * Create a local PIV data set, so input data will not be modified halfway * during substitution. Else, this might affect PIV displacements @@ -489,10 +501,9 @@ * no substitution, only sets peak_no to 0 */ -#pragma omp parallel for +#pragma omp parallel for shared(piv_data, l_data, outlier_found) private(i, j) for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { - guint i_mx, j_mx, i_my, j_my; if (piv_data->peak_no[i][j] != -1 && piv_data->snr[i][j] > valid_par->residu_max) { @@ -514,12 +525,13 @@ } else if (valid_par->subst_type == GPIV_VALID_SUBSTYPE__L_MEAN) { -#pragma omp parallel for +#pragma omp parallel for shared(piv_data) private(i, j, k, l, count, \ + dx_sum, dy_sum) for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { - guint i_mx, j_mx, i_my, j_my; - guint count = 0; - gfloat dx_sum = 0, dy_sum = 0; + count = 0; + dx_sum = 0.0; + dy_sum = 0.0; if (piv_data->peak_no[i][j] != -1 && piv_data->snr[i][j] > valid_par->residu_max) { @@ -555,11 +567,9 @@ /* * Substitution with median particle displacements */ - -#pragma omp parallel for for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { - guint i_mx, j_mx, i_my, j_my; + guint i_mx, j_mx, i_my, j_my; if (piv_data->peak_no[i][j] != -1 && piv_data->snr[i][j] > valid_par->residu_max) { @@ -601,9 +611,6 @@ } } - - - /* @@ -753,7 +760,8 @@ /* * Subdividing fractional particle displacements in bins */ -/* #pragma omp parallel for private(fract) */ +#pragma omp parallel for shared(piv_data, bound, count) \ + private(i, j, k, fract) for (i = 0; i < piv_data->ny; i++) { for (j = 0; j < piv_data->nx; j++) { fract = fabs(piv_data->dx[i][j] - (int) piv_data->dx[i][j]); @@ -897,13 +905,6 @@ l_piv_par = gpiv_piv_cp_parameters (piv_par); l_valid_par = gpiv_valid_cp_parameters (valid_par); -/* l_data = gpiv_alloc_pivdata (piv_data->nx, piv_data->ny); */ -/* for (i = 0; i < piv_data->ny; i++) { */ -/* for (j = 0; j < piv_data->nx; j++) { */ -/* l_data->point_x[i][j] = piv_data->point_x[i][j]; */ -/* l_data->point_y[i][j] = piv_data->point_y[i][j]; */ -/* } */ -/* } */ while (outlier_found && count < GPIV_VALID_MAX_SWEEP) { @@ -936,25 +937,14 @@ /* * piv_data are updated with corrected values - * l_data will be freed for an eventually next loop + * l_data will be made free for an eventually next loop + * BUGFIX: possible memleak? Comment: when outside if {} causes crashing. */ if (outlier_found) { count++; - -#pragma omp parallel for - for (i = 0; i < piv_data->ny; i++) { - for (j = 0; j < piv_data->nx; j++) { - piv_data->point_x[i][j] = l_data->point_x[i][j]; - piv_data->point_y[i][j] = l_data->point_y[i][j]; - piv_data->dx[i][j] = l_data->dx[i][j]; - piv_data->dy[i][j] = l_data->dy[i][j]; - piv_data->snr[i][j] = l_data->snr[i][j]; - piv_data->peak_no[i][j] = l_data->peak_no[i][j]; - } - } - + gpiv_ovwrt_pivdata (l_data, piv_data); + } gpiv_free_pivdata (l_data); - } }