[GRASS-SVN] r70293 - in grass/trunk: lib/rst/interp_float vector/v.surf.rst
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jan 7 12:37:16 PST 2017
Author: annakrat
Date: 2017-01-07 12:37:16 -0800 (Sat, 07 Jan 2017)
New Revision: 70293
Added:
grass/trunk/lib/rst/interp_float/segmen2d_parallel.c
Modified:
grass/trunk/lib/rst/interp_float/Makefile
grass/trunk/lib/rst/interp_float/interpf.h
grass/trunk/lib/rst/interp_float/matrix.c
grass/trunk/lib/rst/interp_float/segmen2d.c
grass/trunk/vector/v.surf.rst/Makefile
grass/trunk/vector/v.surf.rst/main.c
grass/trunk/vector/v.surf.rst/surf.h
Log:
v.surf.rst: added new parallel implementation by S. Zubal and J. Hofierka
Modified: grass/trunk/lib/rst/interp_float/Makefile
===================================================================
--- grass/trunk/lib/rst/interp_float/Makefile 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/lib/rst/interp_float/Makefile 2017-01-07 20:37:16 UTC (rev 70293)
@@ -7,8 +7,10 @@
include $(MODULE_TOPDIR)/include/Make/Lib.make
+LIBES = $(VECTORLIB) $(INTERPDATALIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(QTREELIB) $(BITMAPLIB)
EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(OMPCFLAGS)
+EXTRA_LIBS = $(OMPLIB)
default: $(ARCH_INCDIR)/interpf.h
$(MAKE) lib
Modified: grass/trunk/lib/rst/interp_float/interpf.h
===================================================================
--- grass/trunk/lib/rst/interp_float/interpf.h 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/lib/rst/interp_float/interpf.h 2017-01-07 20:37:16 UTC (rev 70293)
@@ -138,6 +138,8 @@
/* matrix.c */
int IL_matrix_create(struct interp_params *, struct triple *, int, double **,
int *);
+int IL_matrix_create_alloc(struct interp_params *, struct triple *, int, double **,
+ int *, double *);
/* minmax.c */
int min1(int, int);
int max1(int, int);
@@ -176,11 +178,18 @@
double *, double *, double *, double *, double *,
double *, int, int);
/* segmen2d.c */
+double smallest_segment(struct multtree *, int);
int IL_interp_segments_2d(struct interp_params *, struct tree_info *,
struct multtree *, struct BM *, double, double,
double *, double *, double *, double *, double *,
double *, double *, double *, double *, int, off_t,
double);
+/* segmen2d_parallel.c */
+int IL_interp_segments_2d_parallel(struct interp_params *, struct tree_info *,
+ struct multtree *, struct BM *, double, double,
+ double *, double *, double *, double *, double *,
+ double *, double *, double *, double *, int, off_t,
+ double, int);
/* vinput2d.c */
int IL_vector_input_data_2d(struct interp_params *, struct Map_info *, int,
char *, char *, struct tree_info *, double *,
Modified: grass/trunk/lib/rst/interp_float/matrix.c
===================================================================
--- grass/trunk/lib/rst/interp_float/matrix.c 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/lib/rst/interp_float/matrix.c 2017-01-07 20:37:16 UTC (rev 70293)
@@ -36,6 +36,25 @@
#include <grass/gmath.h>
+int IL_matrix_create(struct interp_params *params,
+ struct triple *points, /* points for interpolation */
+ int n_points, /* number of points */
+ double **matrix, /* matrix */
+ int *indx)
+{
+ static double *A = NULL;
+
+ if (!A) {
+ if (!
+ (A =
+ G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
+ fprintf(stderr, "Cannot allocate memory for A\n");
+ return -1;
+ }
+ }
+ return IL_matrix_create_alloc(params, points, n_points, matrix, indx, A);
+}
+
/*!
* \brief Creates system of linear equations from interpolated points
*
@@ -50,18 +69,18 @@
*
* \return -1 on failure, 1 on success
*/
-int IL_matrix_create(struct interp_params *params,
- struct triple *points, /* points for interpolation */
- int n_points, /* number of points */
- double **matrix, /* matrix */
- int *indx)
+int IL_matrix_create_alloc(struct interp_params *params,
+ struct triple *points, /* points for interpolation */
+ int n_points, /* number of points */
+ double **matrix, /* matrix */
+ int *indx,
+ double *A /* temporary matrix unique for all threads */)
{
double xx, yy;
double rfsta2, r;
double d;
int n1, k1, k2, k, i1, l, m, i, j;
double fstar2 = params->fi * params->fi / 4.;
- static double *A = NULL;
double RO, amaxa;
double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
double xxr, yyr;
@@ -77,15 +96,6 @@
n1 = n_points + 1;
- if (!A) {
- if (!
- (A =
- G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
- fprintf(stderr, "Cannot allocate memory for A\n");
- return -1;
- }
- }
-
/*
C GENERATION OF MATRIX
C FIRST COLUMN
Modified: grass/trunk/lib/rst/interp_float/segmen2d.c
===================================================================
--- grass/trunk/lib/rst/interp_float/segmen2d.c 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/lib/rst/interp_float/segmen2d.c 2017-01-07 20:37:16 UTC (rev 70293)
@@ -23,9 +23,7 @@
#include <grass/interpf.h>
#include <grass/gmath.h>
-static double smallest_segment(struct multtree *, int);
-
/*!
* Interpolate recursively a tree of segments
*
@@ -345,7 +343,7 @@
return 1;
}
-static double smallest_segment(struct multtree *tree, int n_leafs)
+double smallest_segment(struct multtree *tree, int n_leafs)
{
static int first_time = 1;
int ii;
Added: grass/trunk/lib/rst/interp_float/segmen2d_parallel.c
===================================================================
--- grass/trunk/lib/rst/interp_float/segmen2d_parallel.c (rev 0)
+++ grass/trunk/lib/rst/interp_float/segmen2d_parallel.c 2017-01-07 20:37:16 UTC (rev 70293)
@@ -0,0 +1,459 @@
+/*!
+ * \file segmen2d.c
+ *
+ * \author H. Mitasova, I. Kosinovsky, D. Gerdes
+ *
+ * \copyright
+ * (C) 1993 by Helena Mitasova and the GRASS Development Team
+ *
+ * \copyright
+ * This program is free software under the
+ * GNU General Public License (>=v2).
+ * Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <omp.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/interpf.h>
+#include <grass/gmath.h>
+
+static int cut_tree(struct multtree *, struct multtree **, int *);
+
+
+/*!
+ * See documentation for IL_interp_segments_2d.
+ * This is a parallel processing implementation.
+ */
+int IL_interp_segments_2d_parallel(struct interp_params *params, struct tree_info *info, /*!< info for the quad tree */
+ struct multtree *tree, /*!< current leaf of the quad tree */
+ struct BM *bitmask, /*!< bitmask */
+ double zmin, double zmax, /*!< min and max input z-values */
+ double *zminac, double *zmaxac, /*!< min and max interp. z-values */
+ double *gmin, double *gmax, /*!< min and max inperp. slope val. */
+ double *c1min, double *c1max, /*!< min and max interp. curv. val. */
+ double *c2min, double *c2max, /*!< min and max interp. curv. val. */
+ double *ertot, /*!< total interplating func. error */
+ int totsegm, /*!< total number of segments */
+ off_t offset1, /*!< offset for temp file writing */
+ double dnorm, int threads)
+{
+ int some_thread_failed = 0;
+ int tid;
+ int i = 0;
+ int j = 0;
+ int i_cnt;
+ int cursegm = 0;
+ double smseg;
+ double ***matrix = NULL;
+ int **indx = NULL;
+ double **b = NULL;
+ double **A = NULL;
+ struct quaddata **data_local;
+ struct multtree **all_leafs;
+
+ all_leafs =
+ (struct multtree **)G_malloc(sizeof(struct multtree *) * totsegm);
+ data_local =
+ (struct quaddata **)G_malloc(sizeof(struct quaddata *) * threads);
+ matrix = (double ***)G_malloc(sizeof(double **) * threads);
+ indx = (int **)G_malloc(sizeof(int *) * threads);
+ b = (double **)G_malloc(sizeof(double *) * threads);
+ A = (double **)G_malloc(sizeof(double *) * threads);
+
+ for (i_cnt = 0; i_cnt < threads; i_cnt++) {
+ if (!
+ (matrix[i_cnt] =
+ G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ }
+
+ for (i_cnt = 0; i_cnt < threads; i_cnt++) {
+ if (!(indx[i_cnt] = G_alloc_ivector(params->KMAX2 + 1))) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ }
+
+ for (i_cnt = 0; i_cnt < threads; i_cnt++) {
+ if (!(b[i_cnt] = G_alloc_vector(params->KMAX2 + 3))) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ }
+
+ for (i_cnt = 0; i_cnt < threads; i_cnt++) {
+ if (!
+ (A[i_cnt] =
+ G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ }
+
+ smseg = smallest_segment(tree, 4);
+ cut_tree(tree, all_leafs, &i);
+
+ G_message(_("Starting parallel work"));
+#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)
+ {
+#pragma omp for schedule(dynamic)
+ for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
+ /* Obtain thread id */
+#if defined(_OPENMP)
+ tid = omp_get_thread_num();
+#endif
+
+ double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
+ temp2;
+ int npt, nptprev, MAXENC;
+ double ew_res, ns_res;
+ int MINPTS;
+ double pr;
+ struct triple *point;
+ struct triple skip_point;
+ int m_skip, skip_index, k, segtest;
+ double xx, yy, zz;
+
+
+ //struct quaddata *data_local;
+
+ ns_res = (((struct quaddata *)(tree->data))->ymax -
+ ((struct quaddata *)(tree->data))->y_orig) /
+ params->nsizr;
+ ew_res =
+ (((struct quaddata *)(tree->data))->xmax -
+ ((struct quaddata *)(tree->data))->x_orig) / params->nsizc;
+
+ if (all_leafs[i_cnt] == NULL) {
+ some_thread_failed = -1;
+ continue;
+ }
+ if (all_leafs[i_cnt]->data == NULL) {
+ some_thread_failed = -1;
+ continue;
+ }
+ if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
+ continue;
+ }
+ else {
+ distx =
+ (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
+ ew_res) * 0.1;
+ disty =
+ (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
+ ns_res) * 0.1;
+ distxp = 0;
+ distyp = 0;
+ xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
+ xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
+ ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
+ ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
+ i = 0;
+ MAXENC = 0;
+ /* data is a window with zero points; some fields don't make sence in this case
+ so they are zero (like resolution,dimentions */
+ /* CHANGE */
+ /* Calcutaing kmin for surrent segment (depends on the size) */
+
+ /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
+ pr = pow(2., (xmx - xmn) / smseg - 1.);
+ MINPTS =
+ params->kmin * (pr /
+ (1 + params->kmin * pr / params->KMAX2));
+ /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
+
+ data_local[tid] =
+ (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
+ xmx + distx, ymx + disty,
+ 0, 0, 0, params->KMAX2);
+ npt =
+ MT_region_data(info, tree, data_local[tid], params->KMAX2,
+ 4);
+
+ while ((npt < MINPTS) || (npt > params->KMAX2)) {
+ if (i >= 70) {
+ G_warning(_("Taking too long to find points for interpolation - "
+ "please change the region to area where your points are. "
+ "Continuing calculations..."));
+ break;
+ }
+ i++;
+ if (npt > params->KMAX2)
+ /* decrease window */
+ {
+ MAXENC = 1;
+ nptprev = npt;
+ temp1 = distxp;
+ distxp = distx;
+ distx = distxp - fabs(distx - temp1) * 0.5;
+ temp2 = distyp;
+ distyp = disty;
+ disty = distyp - fabs(disty - temp2) * 0.5;
+ /* decrease by 50% of a previous change in window */
+ }
+ else {
+ nptprev = npt;
+ temp1 = distyp;
+ distyp = disty;
+ temp2 = distxp;
+ distxp = distx;
+ if (MAXENC) {
+ disty = fabs(disty - temp1) * 0.5 + distyp;
+ distx = fabs(distx - temp2) * 0.5 + distxp;
+ }
+ else {
+ distx += distx;
+ disty += disty;
+ }
+ /* decrease by 50% of extra distance */
+ }
+ data_local[tid]->x_orig = xmn - distx; /* update window */
+ data_local[tid]->y_orig = ymn - disty;
+ data_local[tid]->xmax = xmx + distx;
+ data_local[tid]->ymax = ymx + disty;
+ data_local[tid]->n_points = 0;
+ npt =
+ MT_region_data(info, tree, data_local[tid],
+ params->KMAX2, 4);
+ }
+
+ if (totsegm != 0 && tid == 0) {
+ G_percent(cursegm, totsegm, 1);
+ }
+ data_local[tid]->n_rows =
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows;
+ data_local[tid]->n_cols =
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols;
+
+ /* for printing out overlapping segments */
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig =
+ xmn - distx;
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig =
+ ymn - disty;
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax =
+ xmx + distx;
+ ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax =
+ ymx + disty;
+
+ data_local[tid]->x_orig = xmn;
+ data_local[tid]->y_orig = ymn;
+ data_local[tid]->xmax = xmx;
+ data_local[tid]->ymax = ymx;
+
+ /* allocate memory for CV points only if cv is performed */
+ if (params->cv) {
+ if (!
+ (point =
+ (struct triple *)G_malloc(sizeof(struct triple) *
+ data_local[tid]->
+ n_points))) {
+ G_warning(_("Out of memory"));
+ some_thread_failed = -1;
+ continue;
+ }
+ }
+
+ /*normalize the data so that the side of average segment is about 1m */
+ /* put data_points into point only if CV is performed */
+
+ for (i = 0; i < data_local[tid]->n_points; i++) {
+ data_local[tid]->points[i].x =
+ (data_local[tid]->points[i].x -
+ data_local[tid]->x_orig) / dnorm;
+ data_local[tid]->points[i].y =
+ (data_local[tid]->points[i].y -
+ data_local[tid]->y_orig) / dnorm;
+ if (params->cv) {
+ point[i].x = data_local[tid]->points[i].x; /*cv stuff */
+ point[i].y = data_local[tid]->points[i].y; /*cv stuff */
+ point[i].z = data_local[tid]->points[i].z; /*cv stuff */
+ }
+
+ /* commented out by Helena january 1997 as this is not necessary
+ although it may be useful to put normalization of z back?
+ data->points[i].z = data->points[i].z / dnorm;
+ this made smoothing self-adjusting based on dnorm
+ if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
+ */
+ }
+
+ /* cv stuff */
+ if (params->cv) {
+ m_skip = data_local[tid]->n_points;
+ }
+ else {
+ m_skip = 1;
+ }
+ /* remove after cleanup - this is just for testing */
+ skip_point.x = 0.;
+ skip_point.y = 0.;
+ skip_point.z = 0.;
+
+ for (skip_index = 0; skip_index < m_skip; skip_index++) {
+ if (params->cv) {
+ segtest = 0;
+ j = 0;
+ xx = point[skip_index].x * dnorm +
+ data_local[tid]->x_orig + params->x_orig;
+ yy = point[skip_index].y * dnorm +
+ data_local[tid]->y_orig + params->y_orig;
+ zz = point[skip_index].z;
+ if (xx >= data_local[tid]->x_orig + params->x_orig &&
+ xx <= data_local[tid]->xmax + params->x_orig &&
+ yy >= data_local[tid]->y_orig + params->y_orig &&
+ yy <= data_local[tid]->ymax + params->y_orig) {
+ segtest = 1;
+ skip_point.x = point[skip_index].x;
+ skip_point.y = point[skip_index].y;
+ skip_point.z = point[skip_index].z;
+ for (k = 0; k < m_skip; k++) {
+ if (k != skip_index && params->cv) {
+ data_local[tid]->points[j].x = point[k].x;
+ data_local[tid]->points[j].y = point[k].y;
+ data_local[tid]->points[j].z = point[k].z;
+ j++;
+ }
+ }
+ } /* segment area test */
+ }
+ if (!params->cv) {
+ if ( /* params */
+ IL_matrix_create_alloc(params,
+ data_local[tid]->points,
+ data_local[tid]->
+ n_points, matrix[tid],
+ indx[tid],
+ A[tid]) < 0) {
+ some_thread_failed = -1;
+ continue;
+ }
+ }
+ else if (segtest == 1) {
+ if ( /* params */
+ IL_matrix_create_alloc(params,
+ data_local[tid]->points,
+ data_local[tid]->
+ n_points - 1,
+ matrix[tid], indx[tid],
+ A[tid]) < 0) {
+ some_thread_failed = -1;
+ continue;
+ }
+ }
+ if (!params->cv) {
+ for (i = 0; i < data_local[tid]->n_points; i++) {
+ b[tid][i + 1] = data_local[tid]->points[i].z;
+ }
+ b[tid][0] = 0.;
+ G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
+ indx[tid], b[tid]);
+ /* put here condition to skip error if not needed */
+ params->check_points(params, data_local[tid], b[tid],
+ ertot, zmin, dnorm, skip_point);
+ }
+ else if (segtest == 1) {
+ for (i = 0; i < data_local[tid]->n_points - 1; i++) {
+ b[tid][i + 1] = data_local[tid]->points[i].z;
+ }
+ b[tid][0] = 0.;
+ G_lubksb(matrix[tid], data_local[tid]->n_points,
+ indx[tid], b[tid]);
+ params->check_points(params, data_local[tid], b[tid],
+ ertot, zmin, dnorm, skip_point);
+ }
+ } /*end of cv loop */
+
+
+ if (!params->cv) {
+ if ((params->Tmp_fd_z != NULL) ||
+ (params->Tmp_fd_dx != NULL) ||
+ (params->Tmp_fd_dy != NULL) ||
+ (params->Tmp_fd_xx != NULL) ||
+ (params->Tmp_fd_yy != NULL) ||
+ (params->Tmp_fd_xy != NULL)) {
+#pragma omp critical
+ {
+ if (params->grid_calc
+ (params, data_local[tid], bitmask, zmin, zmax,
+ zminac, zmaxac, gmin, gmax, c1min, c1max,
+ c2min, c2max, ertot, b[tid], offset1,
+ dnorm) < 0) {
+ some_thread_failed = -1;
+ }
+ }
+ }
+ }
+
+ /* show after to catch 100% */
+#pragma omp atomic
+ cursegm++;
+ if (totsegm < cursegm) {
+ G_debug(1, "%d %d", totsegm, cursegm);
+ }
+
+ if (totsegm != 0 && tid == 0) {
+ G_percent(cursegm, totsegm, 1);
+ }
+ /*
+ G_free_matrix(matrix);
+ G_free_ivector(indx);
+ G_free_vector(b);
+ */
+ G_free(data_local[tid]->points);
+ G_free(data_local[tid]);
+ }
+ }
+ } /* All threads join master thread and terminate */
+
+ for (i_cnt = 0; i_cnt < threads; i_cnt++) {
+ G_free(matrix[i_cnt]);
+ G_free(indx[i_cnt]);
+ G_free(b[i_cnt]);
+ G_free(A[i_cnt]);
+ }
+ G_free(all_leafs);
+ G_free(data_local);
+ G_free(matrix);
+ G_free(indx);
+ G_free(b);
+ G_free(A);
+
+ if (some_thread_failed != 0) {
+ return -1;
+ }
+ return 1;
+}
+
+
+/* cut given tree into separate leafs */
+int cut_tree(struct multtree *tree, /* tree we want to cut */
+ struct multtree **cut_leafs, /* array of leafs */
+ int *where_to_add /* index of leaf which will be next */ )
+{
+ if (tree == NULL)
+ return -1;
+ if (tree->data == NULL)
+ return -1;
+ if (((struct quaddata *)(tree->data))->points == NULL) {
+ int i;
+
+ for (i = 0; i < 4; i++) {
+ cut_tree(tree->leafs[i], cut_leafs, where_to_add);
+ }
+ return 1;
+ }
+ else {
+ cut_leafs[*where_to_add] = tree;
+ (*where_to_add)++;
+ return 1;
+ }
+}
Modified: grass/trunk/vector/v.surf.rst/Makefile
===================================================================
--- grass/trunk/vector/v.surf.rst/Makefile 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/vector/v.surf.rst/Makefile 2017-01-07 20:37:16 UTC (rev 70293)
@@ -4,10 +4,10 @@
EXTRA_CLEAN_DIRS=doxygenhtml
-LIBES = $(INTERPFLLIB) $(QTREELIB) $(INTERPDATALIB) $(GMATHLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+LIBES = $(INTERPFLLIB) $(QTREELIB) $(INTERPDATALIB) $(GMATHLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(OMPLIB)
DEPENDENCIES = $(INTERPFLDEP) $(QTREEDEP) $(INTERPDATADEP) $(GMATHDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
-EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(OMPCFLAGS)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass/trunk/vector/v.surf.rst/main.c
===================================================================
--- grass/trunk/vector/v.surf.rst/main.c 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/vector/v.surf.rst/main.c 2017-01-07 20:37:16 UTC (rev 70293)
@@ -11,6 +11,7 @@
* modified by Mitasova in November 1999 (dmax, timestamp update)
* dnorm independent tension - -t flag
* cross-validation -v flag by Jaro Hofierka 2004
+ * Parallel version: S. Zubal 2015
*
* PURPOSE: Surface interpolation from vector point data by splines
* COPYRIGHT: (C) 2003-2009, 2013 by the GRASS Development Team
@@ -21,6 +22,7 @@
*
*****************************************************************************/
+#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
@@ -120,6 +122,9 @@
struct multtree *tree;
int open_check, with_z;
char buf[1024];
+#if defined(_OPENMP)
+ int threads;
+#endif
struct GModule *module;
struct
@@ -127,7 +132,7 @@
struct Option *input, *field, *zcol, *wheresql, *scol, *elev, *slope,
*aspect, *pcurv, *tcurv, *mcurv, *treefile, *overfile, *maskmap,
*dmin, *dmax, *zmult, *fi, *rsm, *segmax, *npmin, *cvdev, *devi,
- *theta, *scalex;
+ *theta, *scalex, *threads;
} parm;
struct
{
@@ -246,6 +251,17 @@
_("Name for output vector map showing overlapping windows");
parm.overfile->guisection = _("Outputs");
+#if defined(_OPENMP)
+ parm.threads = G_define_option();
+ parm.threads->key = "nprocs";
+ parm.threads->type = TYPE_INTEGER;
+ parm.threads->answer = NUM_THREADS;
+ parm.threads->required = NO;
+ parm.threads->description =
+ _("Number of threads for parallel computing");
+ parm.threads->guisection = _("Parameters");
+#endif
+
parm.maskmap = G_define_standard_option(G_OPT_R_INPUT);
parm.maskmap->key = "mask";
parm.maskmap->required = NO;
@@ -383,6 +399,21 @@
treefile = parm.treefile->answer;
overfile = parm.overfile->answer;
+#if defined(_OPENMP)
+ sscanf(parm.threads->answer, "%d", &threads);
+ if (threads < 1)
+ {
+ G_warning(_("<%d> is not valid number of threads. Number of threads will be set on <%d>"),
+ threads, abs(threads));
+ threads = abs(threads);
+ }
+ if (parm.devi->answer && threads > 1) {
+ G_warning(_("Parallel computation disabled when deviation output is required"));
+ threads = 1;
+ }
+ omp_set_num_threads(threads);
+#endif
+
if (devi) {
if (Vect_legal_filename(devi) == -1)
G_fatal_error(_("Output vector map name <%s> is not valid map name"),
@@ -679,6 +710,16 @@
}
ertot = 0.;
+#if defined(_OPENMP)
+ G_message(_("Processing segments in parallel..."));
+ if (IL_interp_segments_2d_parallel(¶ms, info, info->root, bitmask,
+ zmin, zmax, &zminac, &zmaxac, &gmin, &gmax,
+ &c1min, &c1max, &c2min, &c2max, &ertot, totsegm,
+ n_cols, dnorm, threads) < 0) {
+ clean();
+ G_fatal_error(_("Interp_segmets failed"));
+ }
+#else
G_message(_("Processing segments..."));
if (IL_interp_segments_2d(¶ms, info, info->root, bitmask,
zmin, zmax, &zminac, &zmaxac, &gmin, &gmax,
@@ -687,7 +728,7 @@
clean();
G_fatal_error(_("Interp_segmets failed"));
}
-
+#endif
G_free_vector(az);
if (cond1) {
G_free_vector(adx);
Modified: grass/trunk/vector/v.surf.rst/surf.h
===================================================================
--- grass/trunk/vector/v.surf.rst/surf.h 2017-01-07 20:30:24 UTC (rev 70292)
+++ grass/trunk/vector/v.surf.rst/surf.h 2017-01-07 20:37:16 UTC (rev 70293)
@@ -31,5 +31,6 @@
#define TOPPARAM "0"
#define PREPROCESS "0"
#define ZMULT "1.0"
+#define NUM_THREADS "1"
int print_tree(struct multtree *, double, double, struct Map_info *);
More information about the grass-commit
mailing list