r.proj parallelization with openmp

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@51141 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.proj/Makefile

@@ -4,6 +4,9 @@ PGM = r.proj
 
 LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB)
 DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
+
+EXTRA_LIBS = $(OMPLIB)
+EXTRA_CFLAGS = $(OMPCFLAGS)
 EXTRA_INC = $(PROJINC)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make

raster/r.proj/bilinear.c

@@ -18,8 +18,8 @@
 void p_bilinear(struct cache *ibuffer,	/* input buffer                  */
 		void *obufptr,	/* ptr in output buffer          */
 		int cell_type,	/* raster map type of obufptr    */
-		double *col_idx,	/* column index          */
-		double *row_idx,	/* row index                     */
+		double col_idx,	/* column index          */
+		double row_idx,	/* row index                     */
 		struct Cell_head *cellhd	/* information of output map     */
     )
 {
@@ -31,8 +31,8 @@ void p_bilinear(struct cache *ibuffer,	/* input buffer                  */
     FCELL c[2][2];
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx - 0.5);
-    col = (int)floor(*col_idx - 0.5);
+    row = (int)floor(row_idx - 0.5);
+    col = (int)floor(col_idx - 0.5);
 
     /* check for out of bounds - if out of bounds set NULL value and return */
     if (row < 0 || row + 1 >= cellhd->rows || col < 0 || col + 1 >= cellhd->cols) {
@@ -42,17 +42,17 @@ void p_bilinear(struct cache *ibuffer,	/* input buffer                  */
 
     for (i = 0; i < 2; i++)
 	for (j = 0; j < 2; j++) {
-	    const FCELL *cellp = CPTR(ibuffer, row + i, col + j);
-	    if (Rast_is_f_null_value(cellp)) {
+	    const FCELL cell = CVAL(ibuffer, row + i, col + j);
+	    if (Rast_is_f_null_value(&cell)) {
 		Rast_set_null_value(obufptr, 1, cell_type);
 		return;
 	    }
-	    c[i][j] = *cellp;
+	    c[i][j] = cell;
 	}
 
     /* do the interpolation  */
-    t = *col_idx - 0.5 - col;
-    u = *row_idx - 0.5 - row;
+    t = col_idx - 0.5 - col;
+    u = row_idx - 0.5 - row;
 
     result = Rast_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]);
 

raster/r.proj/bilinear_f.c

@@ -15,18 +15,18 @@
 void p_bilinear_f(struct cache *ibuffer,	/* input buffer                  */
 		void *obufptr,	/* ptr in output buffer          */
 		int cell_type,	/* raster map type of obufptr    */
-		double *col_idx,	/* column index          */
-		double *row_idx,	/* row index                     */
+		double col_idx,	/* column index          */
+		double row_idx,	/* row index                     */
 	    struct Cell_head *cellhd	/* cell header of input layer    */
     )
 {
     /* start nearest neighbor to do some basic tests */
     int row, col;		/* row/col of nearest neighbor   */
-    FCELL *cellp, cell;
+    FCELL cell;
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx);
-    col = (int)floor(*col_idx);
+    row = (int)floor(row_idx);
+    col = (int)floor(col_idx);
 
     /* check for out of bounds - if out of bounds set NULL value     */
     if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
@@ -34,13 +34,12 @@ void p_bilinear_f(struct cache *ibuffer,	/* input buffer                  */
         return;
     }
 
-    cellp = CPTR(ibuffer, row, col);
+    cell = CVAL(ibuffer, row, col);
     /* if nearest is null, all the other interps will be null */
-    if (Rast_is_f_null_value(cellp)) {
+    if (Rast_is_f_null_value(&cell)) {
         Rast_set_null_value(obufptr, 1, cell_type);
         return;
     }
-    cell = *cellp;
 
     p_bilinear(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
     /* fallback to nearest if bilinear is null */

raster/r.proj/cubic.c

@@ -21,8 +21,8 @@
 void p_cubic(struct cache *ibuffer,	/* input buffer                  */
 	     void *obufptr,	/* ptr in output buffer          */
 	     int cell_type,	/* raster map type of obufptr    */
-	     double *col_idx,	/* column index (decimal)        */
-	     double *row_idx,	/* row index (decimal)           */
+	     double col_idx,	/* column index (decimal)        */
+	     double row_idx,	/* row index (decimal)           */
 	     struct Cell_head *cellhd	/* information of output map     */
     )
 {
@@ -32,11 +32,11 @@ void p_cubic(struct cache *ibuffer,	/* input buffer                  */
     FCELL t, u;			/* intermediate slope            */
     FCELL result;		/* result of interpolation       */
       FCELL val[4];		/* buffer for temporary values   */
-    FCELL cell[4][4];
+    FCELL c[4][4];
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx - 0.5);
-    col = (int)floor(*col_idx - 0.5);
+    row = (int)floor(row_idx - 0.5);
+    col = (int)floor(col_idx - 0.5);
 
     /* check for out of bounds of map - if out of bounds set NULL value     */
     if (row - 1 < 0 || row + 2 >= cellhd->rows ||
@@ -47,20 +47,20 @@ void p_cubic(struct cache *ibuffer,	/* input buffer                  */
 
     for (i = 0; i < 4; i++)
 	for (j = 0; j < 4; j++) {
-	    const FCELL *cellp = CPTR(ibuffer, row - 1 + i, col - 1 + j);
-	    if (Rast_is_f_null_value(cellp)) {
+	    const FCELL cell = CVAL(ibuffer, row - 1 + i, col - 1 + j);
+	    if (Rast_is_f_null_value(&cell)) {
 		Rast_set_null_value(obufptr, 1, cell_type);
 		return;
 	    }
-	    cell[i][j] = *cellp;
+	    c[i][j] = cell;
 	}
 
     /* do the interpolation  */
-    t = *col_idx - 0.5 - col;
-    u = *row_idx - 0.5 - row;
+    t = col_idx - 0.5 - col;
+    u = row_idx - 0.5 - row;
 
     for (i = 0; i < 4; i++) {
-	const FCELL *tmp = cell[i];
+	const FCELL *tmp = c[i];
 	val[i] = Rast_interp_cubic(t, tmp[0], tmp[1], tmp[2], tmp[3]);
     }
 

raster/r.proj/cubic_f.c

@@ -16,18 +16,18 @@
 void p_cubic_f(struct cache *ibuffer,	/* input buffer                  */
 		void *obufptr,	/* ptr in output buffer          */
 		int cell_type,	/* raster map type of obufptr    */
-		double *col_idx,	/* column index          */
-		double *row_idx,	/* row index                     */
+		double col_idx,	/* column index          */
+		double row_idx,	/* row index                     */
 	    struct Cell_head *cellhd	/* cell header of input layer    */
     )
 {
     /* start nearest neighbor to do some basic tests */
     int row, col;		/* row/col of nearest neighbor   */
-    FCELL *cellp, cell;
+    FCELL cell;
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx);
-    col = (int)floor(*col_idx);
+    row = (int)floor(row_idx);
+    col = (int)floor(col_idx);
 
     /* check for out of bounds - if out of bounds set NULL value     */
     if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
@@ -35,13 +35,12 @@ void p_cubic_f(struct cache *ibuffer,	/* input buffer                  */
         return;
     }
 
-    cellp = CPTR(ibuffer, row, col);
+    cell = CVAL(ibuffer, row, col);
     /* if nearest is null, all the other interps will be null */
-    if (Rast_is_f_null_value(cellp)) {
+    if (Rast_is_f_null_value(&cell)) {
         Rast_set_null_value(obufptr, 1, cell_type);
         return;
     }
-    cell = *cellp;
 
     p_cubic(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
     /* fallback to bilinear if cubic is null */

raster/r.proj/lanczos.c

@@ -20,8 +20,8 @@
 void p_lanczos(struct cache *ibuffer,	/* input buffer                  */
 	     void *obufptr,	/* ptr in output buffer          */
 	     int cell_type,	/* raster map type of obufptr    */
-	     double *col_idx,	/* column index (decimal)        */
-	     double *row_idx,	/* row index (decimal)           */
+	     double col_idx,	/* column index (decimal)        */
+	     double row_idx,	/* row index (decimal)           */
 	     struct Cell_head *cellhd	/* information of output map     */
     )
 {
@@ -30,11 +30,11 @@ void p_lanczos(struct cache *ibuffer,	/* input buffer                  */
     int i, j, k;
     double t, u;			/* intermediate slope            */
     FCELL result;		/* result of interpolation       */
-    DCELL cell[25];
+    DCELL c[25];
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx);
-    col = (int)floor(*col_idx);
+    row = (int)floor(row_idx);
+    col = (int)floor(col_idx);
 
     /* check for out of bounds of map - if out of bounds set NULL value     */
     if (row - 2 < 0 || row + 2 >= cellhd->rows ||
@@ -46,20 +46,20 @@ void p_lanczos(struct cache *ibuffer,	/* input buffer                  */
     k = 0;
     for (i = 0; i < 5; i++) {
 	for (j = 0; j < 5; j++) {
-	    const FCELL *cellp = CPTR(ibuffer, row - 2 + i, col - 2 + j);
-	    if (Rast_is_f_null_value(cellp)) {
+	    const FCELL cell = CVAL(ibuffer, row - 2 + i, col - 2 + j);
+	    if (Rast_is_f_null_value(&cell)) {
 		Rast_set_null_value(obufptr, 1, cell_type);
 		return;
 	    }
-	    cell[k++] = *cellp;
+	    c[k++] = cell;
 	}
     }
 
     /* do the interpolation  */
-    t = *col_idx - 0.5 - col;
-    u = *row_idx - 0.5 - row;
+    t = col_idx - 0.5 - col;
+    u = row_idx - 0.5 - row;
 
-    result = Rast_interp_lanczos(t, u, cell);
+    result = Rast_interp_lanczos(t, u, c);
 
     Rast_set_f_value(obufptr, result, cell_type);
 }
@@ -67,18 +67,18 @@ void p_lanczos(struct cache *ibuffer,	/* input buffer                  */
 void p_lanczos_f(struct cache *ibuffer,	/* input buffer                  */
 	     void *obufptr,	/* ptr in output buffer          */
 	     int cell_type,	/* raster map type of obufptr    */
-	     double *col_idx,	/* column index (decimal)        */
-	     double *row_idx,	/* row index (decimal)           */
+	     double col_idx,	/* column index (decimal)        */
+	     double row_idx,	/* row index (decimal)           */
 	     struct Cell_head *cellhd	/* information of output map     */
     )
 {
     int row;			/* row indices for interp        */
     int col;			/* column indices for interp     */
-    FCELL *cellp, cell;
+    FCELL cell;
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx);
-    col = (int)floor(*col_idx);
+    row = (int)floor(row_idx);
+    col = (int)floor(col_idx);
 
     /* check for out of bounds - if out of bounds set NULL value     */
     if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
@@ -86,13 +86,12 @@ void p_lanczos_f(struct cache *ibuffer,	/* input buffer                  */
         return;
     }
 
-    cellp = CPTR(ibuffer, row, col);
+    cell = CVAL(ibuffer, row, col);
     /* if nearest is null, all the other interps will be null */
-    if (Rast_is_f_null_value(cellp)) {
+    if (Rast_is_f_null_value(&cell)) {
         Rast_set_null_value(obufptr, 1, cell_type);
         return;
     }
-    cell = *cellp;
 
     p_lanczos(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
     /* fallback to bicubic if lanczos is null */

raster/r.proj/main.c

@@ -97,15 +97,13 @@ int main(int argc, char **argv)
       overwrite,		/* Overwrite                    */
       curr_proj;		/* output projection (see gis.h) */
 
-    void *obuffer,		/* buffer that holds one output row     */
-     *obufptr;			/* column ptr in output buffer  */
+    void *obuffer;		/* buffer that holds one output row     */
+
     struct cache *ibuffer;	/* buffer that holds the input map      */
     func interpolate;		/* interpolation routine        */
 
-    double xcoord1, xcoord2,	/* temporary x coordinates      */
-      ycoord1, ycoord2,		/* temporary y coordinates      */
-      col_idx,			/* column index in input matrix */
-      row_idx,			/* row index in input matrix    */
+    double xcoord2,		/* temporary x coordinates      */
+      ycoord2,			/* temporary y coordinates      */
       onorth, osouth,		/* save original border coords  */
       oeast, owest, inorth, isouth, ieast, iwest;
     char north_str[30], south_str[30], east_str[30], west_str[30];
@@ -487,44 +485,51 @@ int main(int argc, char **argv)
 
     cell_size = Rast_cell_size(cell_type);
 
-    xcoord1 = xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
-    /**/ ycoord1 = ycoord2 = outcellhd.north - (outcellhd.ns_res / 2);
-    /**/ G_important_message(_("Projecting..."));
+    xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
+    ycoord2 = outcellhd.north - (outcellhd.ns_res / 2);
+    G_important_message(_("Projecting..."));
     G_percent(0, outcellhd.rows, 2);
 
     for (row = 0; row < outcellhd.rows; row++) {
-	obufptr = obuffer;
+	/* obufptr = obuffer */;
+
+        #pragma omp parallel for schedule (static)
 
 	for (col = 0; col < outcellhd.cols; col++) {
+	    void *obufptr = (void *)((const unsigned char *)obuffer + col * cell_size);
+
+	    double xcoord1 = xcoord2 + (col) * outcellhd.ew_res;
+	    double ycoord1 = ycoord2;
+
 	    /* project coordinates in output matrix to       */
 	    /* coordinates in input matrix                   */
 	    if (pj_do_proj(&xcoord1, &ycoord1, &oproj, &iproj) < 0)
 		Rast_set_null_value(obufptr, 1, cell_type);
 	    else {
 		/* convert to row/column indices of input matrix */
-		col_idx = (xcoord1 - incellhd.west) / incellhd.ew_res;
-		row_idx = (incellhd.north - ycoord1) / incellhd.ns_res;
+
+		/* column index in input matrix */
+		double col_idx = (xcoord1 - incellhd.west) / incellhd.ew_res;
+		/* row index in input matrix    */
+		double row_idx = (incellhd.north - ycoord1) / incellhd.ns_res;
 
 		/* and resample data point               */
 		interpolate(ibuffer, obufptr, cell_type,
-			    &col_idx, &row_idx, &incellhd);
+			    col_idx, row_idx, &incellhd);
 	    }
 
-	    obufptr = G_incr_void_ptr(obufptr, cell_size);
-	    xcoord2 += outcellhd.ew_res;
-	    xcoord1 = xcoord2;
-	    ycoord1 = ycoord2;
+	    /* obufptr = G_incr_void_ptr(obufptr, cell_size); */
 	}
 
 	Rast_put_row(fdo, obuffer, cell_type);
 
-	xcoord1 = xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
+	xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
 	ycoord2 -= outcellhd.ns_res;
-	ycoord1 = ycoord2;
 	G_percent(row, outcellhd.rows - 1, 2);
     }
 
     Rast_close(fdo);
+    release_cache(ibuffer);
 
     if (have_colors > 0) {
 	Rast_write_colors(mapname, G_mapset(), &colr);

raster/r.proj/nearest.c

@@ -12,17 +12,17 @@
 void p_nearest(struct cache *ibuffer,	/* input buffer                  */
 	       void *obufptr,	/* ptr in output buffer          */
 	       int cell_type,	/* raster map type of obufptr    */
-	       double *col_idx,	/* column index in input matrix  */
-	       double *row_idx,	/* row index in input matrix     */
+	       double col_idx,	/* column index in input matrix  */
+	       double row_idx,	/* row index in input matrix     */
 	       struct Cell_head *cellhd	/* cell header of input layer    */
     )
 {
     int row, col;		/* row/col of nearest neighbor   */
-    FCELL *cellp;
+    FCELL cell;
 
     /* cut indices to integer */
-    row = (int)floor(*row_idx);
-    col = (int)floor(*col_idx);
+    row = (int)floor(row_idx);
+    col = (int)floor(col_idx);
 
     /* check for out of bounds - if out of bounds set NULL value     */
     if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
@@ -30,12 +30,12 @@ void p_nearest(struct cache *ibuffer,	/* input buffer                  */
 	return;
     }
 
-    cellp = CPTR(ibuffer, row, col);
+    cell = CVAL(ibuffer, row, col);
 
-    if (Rast_is_f_null_value(cellp)) {
+    if (Rast_is_f_null_value(&cell)) {
 	Rast_set_null_value(obufptr, 1, cell_type);
 	return;
     }
 
-    Rast_set_f_value(obufptr, *cellp, cell_type);
+    Rast_set_f_value(obufptr, cell, cell_type);
 }

raster/r.proj/r.proj.h

@@ -17,6 +17,7 @@ typedef FCELL block[BDIM][BDIM];
 struct cache
 {
     int fd;
+    char *fname;
     int stride;
     int nblocks;
     block **grid;
@@ -24,7 +25,7 @@ struct cache
     int *refs;
 };
 
-typedef void (*func) (struct cache *, void *, int, double *, double *,
+typedef void (*func) (struct cache *, void *, int, double, double,
 		      struct Cell_head *);
 
 struct menu
@@ -38,27 +39,28 @@ extern void bordwalk(struct Cell_head *, struct Cell_head *, struct pj_info *,
 		     struct pj_info *);
 extern struct cache *readcell(int, const char *);
 extern block *get_block(struct cache *, int);
+extern void release_cache(struct cache *);
 
 /* declare resampling methods */
 /* bilinear.c */
-extern void p_bilinear(struct cache *, void *, int, double *, double *,
+extern void p_bilinear(struct cache *, void *, int, double, double,
 		       struct Cell_head *);
 /* cubic.c */
-extern void p_cubic(struct cache *, void *, int, double *, double *,
+extern void p_cubic(struct cache *, void *, int, double, double,
 		    struct Cell_head *);
 /* nearest.c */
-extern void p_nearest(struct cache *, void *, int, double *, double *,
+extern void p_nearest(struct cache *, void *, int, double, double,
 		      struct Cell_head *);
 /* bilinear_f.c */
-extern void p_bilinear_f(struct cache *, void *, int, double *, double *,
+extern void p_bilinear_f(struct cache *, void *, int, double, double,
 		       struct Cell_head *);
 /* cubic_f.c */
-extern void p_cubic_f(struct cache *, void *, int, double *, double *,
+extern void p_cubic_f(struct cache *, void *, int, double, double,
 		    struct Cell_head *);
 /* lanczos.c */
-extern void p_lanczos(struct cache *, void *, int, double *, double *,
+extern void p_lanczos(struct cache *, void *, int, double, double,
 		    struct Cell_head *);
-extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
+extern void p_lanczos_f(struct cache *, void *, int, double, double,
 		    struct Cell_head *);
 
 #if 1
@@ -66,7 +68,7 @@ extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
 #define BKIDX(c,y,x) ((y) * (c)->stride + (x))
 #define BKPTR(c,y,x) ((c)->grid[BKIDX((c),(y),(x))])
 #define BLOCK(c,y,x) (BKPTR((c),(y),(x)) ? BKPTR((c),(y),(x)) : get_block((c),BKIDX((c),(y),(x))))
-#define CPTR(c,y,x) (&(*BLOCK((c),HI((y)),HI((x))))[LO((y))][LO((x))])
+#define CVAL(c,y,x) ((*BLOCK((c),HI((y)),HI((x))))[LO((y))][LO((x))])
 
 #else
 

raster/r.proj/readcell.c

@@ -22,7 +22,6 @@ struct cache *readcell(int fdi, const char *size)
     int nrows;
     int ncols;
     int row;
-    char *filename;
     int nx, ny;
     int nblocks;
     int i;
@@ -46,20 +45,21 @@ struct cache *readcell(int fdi, const char *size)
     c->nblocks = nblocks;
     c->grid = (block **) G_calloc(nx * ny, sizeof(block *));
     c->blocks = (block *) G_malloc(nblocks * sizeof(block));
-    c->refs = (int *)G_malloc(nblocks * sizeof(int));
+    c->refs = (int *)G_calloc(nblocks, sizeof(int));
 
     if (nblocks < nx * ny) {
 	/* Temporary file must be created in output location */
 	G__switch_env();
-	filename = G_tempfile();
+	c->fname = G_tempfile();
 	G__switch_env();
-	c->fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+	c->fd = open(c->fname, O_RDWR | O_CREAT | O_EXCL, 0600);
 	if (c->fd < 0)
 	    G_fatal_error(_("Unable to open temporary file"));
-	remove(filename);
     }
-    else
+    else {
 	c->fd = -1;
+	c->fname = NULL;
+    }
 
     G_important_message(_("Allocating memory and reading input map..."));
     G_percent(0, nrows, 5);
@@ -97,36 +97,57 @@ struct cache *readcell(int fdi, const char *size)
 
     G_free(tmpbuf);
 
-    if (c->fd < 0)
+    if (c->fd < 0) {
 	for (i = 0; i < c->nblocks; i++) {
 	    c->grid[i] = &c->blocks[i];
 	    c->refs[i] = i;
 	}
+    }
+    else {
+	close(c->fd);
+	c->fd = -1;
+    }
 
     return c;
 }
 
 block *get_block(struct cache * c, int idx)
 {
+    int fd;
     int replace = rand() % c->nblocks;
     block *p = &c->blocks[replace];
     int ref = c->refs[replace];
     off_t offset = (off_t) idx * sizeof(FCELL) << L2BSIZE;
 
-    if (c->fd < 0)
+    if (c->fname == NULL)
 	G_fatal_error(_("Internal error: cache miss on fully-cached map"));
 
+    fd = open(c->fname, O_RDONLY, 0600);
+    if (fd < 0)
+	G_fatal_error(_("Unable to open temporary file"));
+
     if (ref >= 0)
 	c->grid[ref] = NULL;
 
     c->grid[idx] = p;
     c->refs[replace] = idx;
 
-    if (lseek(c->fd, offset, SEEK_SET) < 0)
+    if (lseek(fd, offset, SEEK_SET) < 0)
 	G_fatal_error(_("Error seeking on segment file"));
 
-    if (read(c->fd, p, sizeof(block)) < 0)
-	G_fatal_error(_("Error writing segment file"));
+    if (read(fd, p, sizeof(block)) < 0)
+	G_fatal_error(_("Error reading segment file"));
+	
+    close(fd);
 
     return p;
 }
+
+void release_cache(struct cache *c)
+{
+    G_free(c->grid);
+    G_free(c->blocks);
+    G_free(c->refs);
+    remove(c->fname);
+    G_free(c);
+}