GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
segmen2d_parallel.c
Go to the documentation of this file.
1/*!
2 * \file segmen2d.c
3 *
4 * \author H. Mitasova, I. Kosinovsky, D. Gerdes (single core)
5 * \author Stanislav Zubal, Michal Lacko (OpenMP version)
6 * \author Anna Petrasova (OpenMP version GRASS integration)
7 *
8 * \copyright
9 * (C) 1993-2017 by Helena Mitasova and the GRASS Development Team
10 *
11 * \copyright
12 * This program is free software under the
13 * GNU General Public License (>=v2).
14 * Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 */
18
19
20#include <stdio.h>
21#include <stdlib.h>
22#include <math.h>
23#if defined(_OPENMP)
24#include <omp.h>
25#endif
26#include <grass/gis.h>
27#include <grass/glocale.h>
28#include <grass/interpf.h>
29#include <grass/gmath.h>
30
31static int cut_tree(struct multtree *, struct multtree **, int *);
32
33
34/*!
35 * See documentation for IL_interp_segments_2d.
36 * This is a parallel processing implementation.
37 */
38int IL_interp_segments_2d_parallel(struct interp_params *params, struct tree_info *info, /*!< info for the quad tree */
39 struct multtree *tree, /*!< current leaf of the quad tree */
40 struct BM *bitmask, /*!< bitmask */
41 double zmin, double zmax, /*!< min and max input z-values */
42 double *zminac, double *zmaxac, /*!< min and max interp. z-values */
43 double *gmin, double *gmax, /*!< min and max inperp. slope val. */
44 double *c1min, double *c1max, /*!< min and max interp. curv. val. */
45 double *c2min, double *c2max, /*!< min and max interp. curv. val. */
46 double *ertot, /*!< total interplating func. error */
47 int totsegm, /*!< total number of segments */
48 off_t offset1, /*!< offset for temp file writing */
49 double dnorm, int threads)
50{
51 int some_thread_failed = 0;
52 int tid;
53 int i = 0;
54 int j = 0;
55 int i_cnt;
56 int cursegm = 0;
57 double smseg;
58 double ***matrix = NULL;
59 int **indx = NULL;
60 double **b = NULL;
61 double **A = NULL;
62 struct quaddata **data_local;
63 struct multtree **all_leafs;
64
65 all_leafs =
66 (struct multtree **)G_malloc(sizeof(struct multtree *) * totsegm);
67 data_local =
68 (struct quaddata **)G_malloc(sizeof(struct quaddata *) * threads);
69 matrix = (double ***)G_malloc(sizeof(double **) * threads);
70 indx = (int **)G_malloc(sizeof(int *) * threads);
71 b = (double **)G_malloc(sizeof(double *) * threads);
72 A = (double **)G_malloc(sizeof(double *) * threads);
73
74 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
75 if (!
76 (matrix[i_cnt] =
77 G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
78 G_fatal_error(_("Out of memory"));
79 return -1;
80 }
81 }
82
83 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
84 if (!(indx[i_cnt] = G_alloc_ivector(params->KMAX2 + 1))) {
85 G_fatal_error(_("Out of memory"));
86 return -1;
87 }
88 }
89
90 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
91 if (!(b[i_cnt] = G_alloc_vector(params->KMAX2 + 3))) {
92 G_fatal_error(_("Out of memory"));
93 return -1;
94 }
95 }
96
97 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
98 if (!
99 (A[i_cnt] =
100 G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
101 G_fatal_error(_("Out of memory"));
102 return -1;
103 }
104 }
105
106 smseg = smallest_segment(tree, 4);
107 cut_tree(tree, all_leafs, &i);
108
109 G_message(_("Starting parallel work"));
110#pragma omp parallel firstprivate(tid, i, j, zmin, zmax, tree, totsegm, offset1, dnorm, smseg, ertot, params, info, all_leafs, bitmask, b, indx, matrix, data_local, A) shared(cursegm, threads, some_thread_failed, zminac, zmaxac, gmin, gmax, c1min, c1max, c2min, c2max) default(none)
111 {
112#pragma omp for schedule(dynamic)
113 for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
114 /* Obtain thread id */
115#if defined(_OPENMP)
116 tid = omp_get_thread_num();
117#endif
118
119 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
120 temp2;
121 int npt, nptprev, MAXENC;
122 double ew_res, ns_res;
123 int MINPTS;
124 double pr;
125 struct triple *point;
126 struct triple skip_point;
127 int m_skip, skip_index, k, segtest;
128 double xx, yy, zz;
129
130
131 //struct quaddata *data_local;
132
133 ns_res = (((struct quaddata *)(tree->data))->ymax -
134 ((struct quaddata *)(tree->data))->y_orig) /
135 params->nsizr;
136 ew_res =
137 (((struct quaddata *)(tree->data))->xmax -
138 ((struct quaddata *)(tree->data))->x_orig) / params->nsizc;
139
140 if (all_leafs[i_cnt] == NULL) {
141 some_thread_failed = -1;
142 continue;
143 }
144 if (all_leafs[i_cnt]->data == NULL) {
145 some_thread_failed = -1;
146 continue;
147 }
148 if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
149 continue;
150 }
151 else {
152 distx =
153 (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
154 ew_res) * 0.1;
155 disty =
156 (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
157 ns_res) * 0.1;
158 distxp = 0;
159 distyp = 0;
160 xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
161 xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
162 ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
163 ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
164 i = 0;
165 MAXENC = 0;
166 /* data is a window with zero points; some fields don't make sense in this case
167 so they are zero (like resolution,dimensions */
168 /* CHANGE */
169 /* Calcutaing kmin for surrent segment (depends on the size) */
170
171 /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
172 pr = pow(2., (xmx - xmn) / smseg - 1.);
173 MINPTS =
174 params->kmin * (pr /
175 (1 + params->kmin * pr / params->KMAX2));
176 /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
177
178 data_local[tid] =
179 (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
180 xmx + distx, ymx + disty,
181 0, 0, 0, params->KMAX2);
182 npt =
183 MT_region_data(info, tree, data_local[tid], params->KMAX2,
184 4);
185
186 while ((npt < MINPTS) || (npt > params->KMAX2)) {
187 if (i >= 70) {
188 G_warning(_("Taking too long to find points for interpolation - "
189 "please change the region to area where your points are. "
190 "Continuing calculations..."));
191 break;
192 }
193 i++;
194 if (npt > params->KMAX2)
195 /* decrease window */
196 {
197 MAXENC = 1;
198 nptprev = npt;
199 temp1 = distxp;
200 distxp = distx;
201 distx = distxp - fabs(distx - temp1) * 0.5;
202 temp2 = distyp;
203 distyp = disty;
204 disty = distyp - fabs(disty - temp2) * 0.5;
205 /* decrease by 50% of a previous change in window */
206 }
207 else {
208 nptprev = npt;
209 temp1 = distyp;
210 distyp = disty;
211 temp2 = distxp;
212 distxp = distx;
213 if (MAXENC) {
214 disty = fabs(disty - temp1) * 0.5 + distyp;
215 distx = fabs(distx - temp2) * 0.5 + distxp;
216 }
217 else {
218 distx += distx;
219 disty += disty;
220 }
221 /* decrease by 50% of extra distance */
222 }
223 data_local[tid]->x_orig = xmn - distx; /* update window */
224 data_local[tid]->y_orig = ymn - disty;
225 data_local[tid]->xmax = xmx + distx;
226 data_local[tid]->ymax = ymx + disty;
227 data_local[tid]->n_points = 0;
228 npt =
229 MT_region_data(info, tree, data_local[tid],
230 params->KMAX2, 4);
231 }
232
233 if (totsegm != 0 && tid == 0) {
234 G_percent(cursegm, totsegm, 1);
235 }
236 data_local[tid]->n_rows =
237 ((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows;
238 data_local[tid]->n_cols =
239 ((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols;
240
241 /* for printing out overlapping segments */
242 ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig =
243 xmn - distx;
244 ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig =
245 ymn - disty;
246 ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax =
247 xmx + distx;
248 ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax =
249 ymx + disty;
250
251 data_local[tid]->x_orig = xmn;
252 data_local[tid]->y_orig = ymn;
253 data_local[tid]->xmax = xmx;
254 data_local[tid]->ymax = ymx;
255
256 /* allocate memory for CV points only if cv is performed */
257 if (params->cv) {
258 if (!
259 (point =
260 (struct triple *)G_malloc(sizeof(struct triple) *
261 data_local[tid]->
262 n_points))) {
263 G_warning(_("Out of memory"));
264 some_thread_failed = -1;
265 continue;
266 }
267 }
268
269 /*normalize the data so that the side of average segment is about 1m */
270 /* put data_points into point only if CV is performed */
271
272 for (i = 0; i < data_local[tid]->n_points; i++) {
273 data_local[tid]->points[i].x =
274 (data_local[tid]->points[i].x -
275 data_local[tid]->x_orig) / dnorm;
276 data_local[tid]->points[i].y =
277 (data_local[tid]->points[i].y -
278 data_local[tid]->y_orig) / dnorm;
279 if (params->cv) {
280 point[i].x = data_local[tid]->points[i].x; /*cv stuff */
281 point[i].y = data_local[tid]->points[i].y; /*cv stuff */
282 point[i].z = data_local[tid]->points[i].z; /*cv stuff */
283 }
284
285 /* commented out by Helena january 1997 as this is not necessary
286 although it may be useful to put normalization of z back?
287 data->points[i].z = data->points[i].z / dnorm;
288 this made smoothing self-adjusting based on dnorm
289 if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
290 */
291 }
292
293 /* cv stuff */
294 if (params->cv) {
295 m_skip = data_local[tid]->n_points;
296 }
297 else {
298 m_skip = 1;
299 }
300 /* remove after cleanup - this is just for testing */
301 skip_point.x = 0.;
302 skip_point.y = 0.;
303 skip_point.z = 0.;
304
305 for (skip_index = 0; skip_index < m_skip; skip_index++) {
306 if (params->cv) {
307 segtest = 0;
308 j = 0;
309 xx = point[skip_index].x * dnorm +
310 data_local[tid]->x_orig + params->x_orig;
311 yy = point[skip_index].y * dnorm +
312 data_local[tid]->y_orig + params->y_orig;
313 zz = point[skip_index].z;
314 if (xx >= data_local[tid]->x_orig + params->x_orig &&
315 xx <= data_local[tid]->xmax + params->x_orig &&
316 yy >= data_local[tid]->y_orig + params->y_orig &&
317 yy <= data_local[tid]->ymax + params->y_orig) {
318 segtest = 1;
319 skip_point.x = point[skip_index].x;
320 skip_point.y = point[skip_index].y;
321 skip_point.z = point[skip_index].z;
322 for (k = 0; k < m_skip; k++) {
323 if (k != skip_index && params->cv) {
324 data_local[tid]->points[j].x = point[k].x;
325 data_local[tid]->points[j].y = point[k].y;
326 data_local[tid]->points[j].z = point[k].z;
327 j++;
328 }
329 }
330 } /* segment area test */
331 }
332 if (!params->cv) {
333 if ( /* params */
335 data_local[tid]->points,
336 data_local[tid]->
337 n_points, matrix[tid],
338 indx[tid],
339 A[tid]) < 0) {
340 some_thread_failed = -1;
341 continue;
342 }
343 }
344 else if (segtest == 1) {
345 if ( /* params */
347 data_local[tid]->points,
348 data_local[tid]->
349 n_points - 1,
350 matrix[tid], indx[tid],
351 A[tid]) < 0) {
352 some_thread_failed = -1;
353 continue;
354 }
355 }
356 if (!params->cv) {
357 for (i = 0; i < data_local[tid]->n_points; i++) {
358 b[tid][i + 1] = data_local[tid]->points[i].z;
359 }
360 b[tid][0] = 0.;
361 G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
362 indx[tid], b[tid]);
363 /* put here condition to skip error if not needed */
364 params->check_points(params, data_local[tid], b[tid],
365 ertot, zmin, dnorm, skip_point);
366 }
367 else if (segtest == 1) {
368 for (i = 0; i < data_local[tid]->n_points - 1; i++) {
369 b[tid][i + 1] = data_local[tid]->points[i].z;
370 }
371 b[tid][0] = 0.;
372 G_lubksb(matrix[tid], data_local[tid]->n_points,
373 indx[tid], b[tid]);
374 params->check_points(params, data_local[tid], b[tid],
375 ertot, zmin, dnorm, skip_point);
376 }
377 } /*end of cv loop */
378
379
380 if (!params->cv) {
381 if ((params->Tmp_fd_z != NULL) ||
382 (params->Tmp_fd_dx != NULL) ||
383 (params->Tmp_fd_dy != NULL) ||
384 (params->Tmp_fd_xx != NULL) ||
385 (params->Tmp_fd_yy != NULL) ||
386 (params->Tmp_fd_xy != NULL)) {
387#pragma omp critical
388 {
389 if (params->grid_calc
390 (params, data_local[tid], bitmask, zmin, zmax,
391 zminac, zmaxac, gmin, gmax, c1min, c1max,
392 c2min, c2max, ertot, b[tid], offset1,
393 dnorm) < 0) {
394 some_thread_failed = -1;
395 }
396 }
397 }
398 }
399
400 /* show after to catch 100% */
401#pragma omp atomic
402 cursegm++;
403 if (totsegm < cursegm) {
404 G_debug(1, "%d %d", totsegm, cursegm);
405 }
406
407 if (totsegm != 0 && tid == 0) {
408 G_percent(cursegm, totsegm, 1);
409 }
410 /*
411 G_free_matrix(matrix);
412 G_free_ivector(indx);
413 G_free_vector(b);
414 */
415 G_free(data_local[tid]->points);
416 G_free(data_local[tid]);
417 }
418 }
419 } /* All threads join master thread and terminate */
420
421 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
422 G_free(matrix[i_cnt]);
423 G_free(indx[i_cnt]);
424 G_free(b[i_cnt]);
425 G_free(A[i_cnt]);
426 }
427 G_free(all_leafs);
428 G_free(data_local);
429 G_free(matrix);
430 G_free(indx);
431 G_free(b);
432 G_free(A);
433
434 if (some_thread_failed != 0) {
435 return -1;
436 }
437 return 1;
438}
439
440
441/* cut given tree into separate leafs */
442int cut_tree(struct multtree *tree, /* tree we want to cut */
443 struct multtree **cut_leafs, /* array of leafs */
444 int *where_to_add /* index of leaf which will be next */ )
445{
446 if (tree == NULL)
447 return -1;
448 if (tree->data == NULL)
449 return -1;
450 if (((struct quaddata *)(tree->data))->points == NULL) {
451 int i;
452
453 for (i = 0; i < 4; i++) {
454 cut_tree(tree->leafs[i], cut_leafs, where_to_add);
455 }
456 return 1;
457 }
458 else {
459 cut_leafs[*where_to_add] = tree;
460 (*where_to_add)++;
461 return 1;
462 }
463}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#define NULL
Definition: ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition: dataquad.c:62
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double b
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition: ialloc.c:41
int IL_matrix_create_alloc(struct interp_params *, struct triple *, int, double **, int *, double *)
Creates system of linear equations from interpolated points.
Definition: matrix.c:72
double smallest_segment(struct multtree *, int)
Definition: segmen2d.c:346
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:104
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition: qtree.c:194
int IL_interp_segments_2d_parallel(struct interp_params *params, struct tree_info *info, struct multtree *tree, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot, int totsegm, off_t offset1, double dnorm, int threads)
check_points_fn * check_points
Definition: interpf.h:99
FILE * Tmp_fd_xx
Definition: interpf.h:93
FILE * Tmp_fd_xy
Definition: interpf.h:94
FILE * Tmp_fd_yy
Definition: interpf.h:94
grid_calc_fn * grid_calc
Definition: interpf.h:97
double x_orig
Definition: interpf.h:87
FILE * Tmp_fd_dx
Definition: interpf.h:92
double y_orig
Definition: interpf.h:87
FILE * Tmp_fd_z
Definition: interpf.h:92
FILE * Tmp_fd_dy
Definition: interpf.h:93
Definition: qtree.h:57
struct multtree ** leafs
Definition: qtree.h:59
struct quaddata * data
Definition: qtree.h:58
double ymax
Definition: dataquad.h:53
double y_orig
Definition: dataquad.h:51
double x_orig
Definition: dataquad.h:50
struct triple * points
Definition: dataquad.h:57
int n_points
Definition: dataquad.h:56
double xmax
Definition: dataquad.h:52
int n_cols
Definition: dataquad.h:55
int n_rows
Definition: dataquad.h:54
double z
Definition: dataquad.h:44
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43