27 #include <grass/interpf.h>
41 double zmin,
double zmax,
42 double *zminac,
double *zmaxac,
43 double *gmin,
double *gmax,
44 double *c1min,
double *c1max,
45 double *c2min,
double *c2max,
49 double dnorm,
int threads)
51 int some_thread_failed = 0;
58 double ***matrix =
NULL;
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);
74 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
82 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
89 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
96 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
98 (params->
KMAX2 + 2) * (params->
KMAX2 + 2) + 1))) {
105 cut_tree(tree, all_leafs, &i);
108 #pragma omp parallel firstprivate( \
109 tid, i, j, zmin, zmax, tree, totsegm, offset1, dnorm, smseg, ertot, \
110 params, info, all_leafs, bitmask, b, indx, matrix, data_local, A) \
111 shared(cursegm, threads, some_thread_failed, zminac, zmaxac, gmin, gmax, \
112 c1min, c1max, c2min, c2max) default(none)
114 #pragma omp for schedule(dynamic)
115 for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
118 tid = omp_get_thread_num();
121 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
124 double ew_res, ns_res;
129 int m_skip, skip_index, k, segtest;
141 if (all_leafs[i_cnt] ==
NULL) {
142 some_thread_failed = -1;
145 if (all_leafs[i_cnt]->data ==
NULL) {
146 some_thread_failed = -1;
162 xmx = ((
struct quaddata *)(all_leafs[i_cnt]->data))->
xmax;
164 ymx = ((
struct quaddata *)(all_leafs[i_cnt]->data))->
ymax;
174 pr = pow(2., (xmx - xmn) / smseg - 1.);
175 MINPTS = params->
kmin *
176 (pr / (1 + params->
kmin * pr / params->
KMAX2));
182 xmn - distx, ymn - disty, xmx + distx, ymx + disty, 0, 0, 0,
187 while ((npt < MINPTS) || (npt > params->
KMAX2)) {
189 G_warning(
_(
"Taking too long to find points for "
191 "please change the region to area where "
193 "Continuing calculations..."));
197 if (npt > params->
KMAX2)
203 distx = distxp - fabs(distx - temp1) * 0.5;
206 disty = distyp - fabs(disty - temp2) * 0.5;
215 disty = fabs(disty - temp1) * 0.5 + distyp;
216 distx = fabs(distx - temp2) * 0.5 + distxp;
224 data_local[tid]->
x_orig = xmn - distx;
225 data_local[tid]->
y_orig = ymn - disty;
226 data_local[tid]->
xmax = xmx + distx;
227 data_local[tid]->
ymax = ymx + disty;
233 if (totsegm != 0 && tid == 0) {
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;
262 some_thread_failed = -1;
271 for (i = 0; i < data_local[tid]->
n_points; i++) {
273 (data_local[tid]->
points[i].
x -
274 data_local[tid]->
x_orig) /
277 (data_local[tid]->
points[i].
y -
278 data_local[tid]->
y_orig) /
281 point[i].
x = data_local[tid]->
points[i].
x;
282 point[i].
y = data_local[tid]->
points[i].
y;
283 point[i].
z = data_local[tid]->
points[i].
z;
307 for (skip_index = 0; skip_index < m_skip; skip_index++) {
311 xx = point[skip_index].
x * dnorm +
313 yy = point[skip_index].
y * dnorm +
317 xx <= data_local[tid]->
xmax + params->
x_orig &&
319 yy <= data_local[tid]->
ymax + params->
y_orig) {
321 skip_point.
x = point[skip_index].
x;
322 skip_point.
y = point[skip_index].
y;
323 skip_point.
z = point[skip_index].
z;
324 for (k = 0; k < m_skip; k++) {
325 if (k != skip_index && params->
cv) {
326 data_local[tid]->
points[j].
x = point[k].
x;
327 data_local[tid]->
points[j].
y = point[k].
y;
328 data_local[tid]->
points[j].
z = point[k].
z;
337 params, data_local[tid]->
points,
338 data_local[tid]->
n_points, matrix[tid],
339 indx[tid], A[tid]) < 0) {
340 some_thread_failed = -1;
344 else if (segtest == 1) {
347 params, data_local[tid]->
points,
348 data_local[tid]->
n_points - 1, matrix[tid],
349 indx[tid], A[tid]) < 0) {
350 some_thread_failed = -1;
355 for (i = 0; i < data_local[tid]->
n_points; i++) {
356 b[tid][i + 1] = data_local[tid]->
points[i].
z;
363 ertot, zmin, dnorm, skip_point);
365 else if (segtest == 1) {
366 for (i = 0; i < data_local[tid]->
n_points - 1; i++) {
367 b[tid][i + 1] = data_local[tid]->
points[i].
z;
373 ertot, zmin, dnorm, skip_point);
386 if (params->
grid_calc(params, data_local[tid],
387 bitmask, zmin, zmax, zminac,
388 zmaxac, gmin, gmax, c1min,
389 c1max, c2min, c2max, ertot,
390 b[tid], offset1, dnorm) < 0) {
391 some_thread_failed = -1;
400 if (totsegm < cursegm) {
401 G_debug(1,
"%d %d", totsegm, cursegm);
404 if (totsegm != 0 && tid == 0) {
418 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
431 if (some_thread_failed != 0) {
449 for (i = 0; i < 4; i++) {
450 cut_tree(tree->
leafs[i], cut_leafs, where_to_add);
455 cut_leafs[*where_to_add] = tree;
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
void G_percent(long, long, int)
Print percent complete messages.
void G_free(void *)
Free allocated memory.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
void G_message(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int * G_alloc_ivector(size_t)
Vector matrix memory allocation.
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
int IL_matrix_create_alloc(struct interp_params *, struct triple *, int, double **, int *, double *)
Creates system of linear equations from interpolated points.
double smallest_segment(struct multtree *, int)
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
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