From c391dea58e97e82b3d39c3ca6e3c674866c12358 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 29 Aug 2021 14:57:47 -1000 Subject: [PATCH 1/7] Update greenspline -C for internediate grids Currently, the modifiers +m and +M are mutually exclusive and thus it is only possible to output intermediary grids for the cumulative or incremental solution. However, since the main calculation time is in determining the solution, it is helpful to be able to write both types of grids at the same time. This PR relaces the +m and +M with two new modifiers that both can be set: +c for cumulative grids and +i for incremental grids. We now use the regular -Ggrdfile setting and automaticallly insert _inc_### or _cum_### before the file extension, and we start with eigenvalue 0 rather than calling it 1 (averything else in GMT starts at 0 so starting at 1 causes grief in batch and movie). The old +m|M modifier is honored in a backwards compatiple way with the original file template and mumbering scheme. --- doc/rst/source/greenspline.rst | 21 +++++---- src/greenspline.c | 85 ++++++++++++++++++++++++---------- 2 files changed, 73 insertions(+), 33 deletions(-) diff --git a/doc/rst/source/greenspline.rst b/doc/rst/source/greenspline.rst index 23374693966..a47ea95d0a5 100644 --- a/doc/rst/source/greenspline.rst +++ b/doc/rst/source/greenspline.rst @@ -15,7 +15,7 @@ Synopsis **gmt greenspline** [ *table* ] |-G|\ *grdfile* [ |-A|\ *gradfile*\ **+f**\ **1**\|\ **2**\|\ **3**\|\ **4**\|\ **5** ] -[ |-C|\ [**n**]\ *value*\ [%][**+f**\ *file*][**+m**\|\ **M**] ] +[ |-C|\ [**n**]\ *value*\ [%][**+c**][**+f**\ *file*][**+i**] ] [ |SYN_OPT-D3| ] [ |-E|\ [*misfitfile*] ] [ |-I|\ *xinc*\ [/*yinc*\ [/*zinc*]] ] @@ -142,7 +142,7 @@ Optional Arguments .. _-C: -**-C**\ [**n**]\ *value*\ [%][**+f**\ *file*][**+m**\|\ **M**] +**-C**\ [**n**]\ *value*\ [%][**+c**][**+f**\ *file*][**+i**] Find an approximate surface fit: Solve the linear system for the spline coefficients by SVD and eliminate the contribution from all eigenvalues whose ratio to the largest eigenvalue is less than *value* @@ -153,12 +153,15 @@ Optional Arguments execution will stop after saving the eigenvalues, i.e., no surface output is produced. Specify **-Cn** to retain only the *value* largest eigenvalues; append % if *value* is the percentage of eigenvalues - to use instead. The two last modifiers (**+m**\|\ **M**) are only + to use instead. The two last modifiers (**+c** and **i**) are only available for 2-D gridding and can be used to write intermediate grids, - one per eigenvalue, and thus require a file name template with a C-format - integer specification to be given via **-G**. The **+m** modifier will - write the contributions to the grid for each eigenvalue, while **+M** - will instead produce the cumulative sum of these contributions. + one per eigenvalue, and thus require a file name with a suitable extension + to be given via **-G** (we automatically insert "_cum_###" or "_inc_###" + before the extension, using a fixed integer format for the eigenvalue + number starting at 0). The **+i** modifier will write the **i**\ ncremental + contributions to the grid for each eigenvalue, while **+c** will instead + produce the **c**\ umulative sum of these contributions. Use both modifiers + to write both types of intermediate grids. .. _-D: @@ -369,9 +372,9 @@ of the surface slope in the NW direction, try:: gmt greenspline @Table_5_11.txt -R0/6.5/-0.2/6.5 -I0.1 -Sr0.95 -V -Z1 -Q-45 -Gslopes.nc To use Cartesian cubic splines and evaluate the cumulative solution as a function of eigenvalue, -using the output template with three digits for the eigenvalue, try:: +using output file based on the main grid name (such as contribution_cum_033.nc), try:: - gmt greenspline @Table_5_11.txt -R0/6.5/-0.2/6.5 -I0.1 -Gcontribution_%3.3d.nc -Sc -Z1 -C+M + gmt greenspline @Table_5_11.txt -R0/6.5/-0.2/6.5 -I0.1 -Gcontribution.nc -Sc -Z1 -C+c Finally, to use Cartesian minimum curvature splines in recovering a surface where the input data is a single surface value (pt.txt) and the diff --git a/src/greenspline.c b/src/greenspline.c index 78e7759041c..0ee104dc4cc 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -227,7 +227,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *C) { /* De static int usage (struct GMTAPI_CTRL *API, int level) { const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); - GMT_Usage (API, 0, "usage: %s [] -G [-A+f] [-C[n][%%][+f][+m|M]] " + GMT_Usage (API, 0, "usage: %s [
] -G [-A+f] [-C[n][%%][+c][+f][+i]] " "[-D] [-E[]] [-I[/[/]]] [-L] [-N] [-Q] " "[-R//[//]]] [-Sc|l|t|r|p|q[]] [-T] " "[%s] [-W[w]] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s [%s]%s[%s] [%s]\n", @@ -257,17 +257,17 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "4: X, Vcomponents."); GMT_Usage (API, 3, "5: X, Vunit-vector, Vmagnitude."); GMT_Usage (API, -2, "Here, X = (x, y[, z]) is the position vector, V = (Vx, Vy[, Vz]) is the gradient vector."); - GMT_Usage (API, 1, "\n-C[n][%%][+f][+m|M]"); + GMT_Usage (API, 1, "\n-C[n][%%][+c][+f][+i]"); GMT_Usage (API, -2, "Solve by SVD and eliminate eigenvalues whose ratio to largest eigenvalue is less than [0]. " "A negative cutoff will stop execution after reporting the eigenvalues. " "Use -Cn to select only the largest eigenvalues [all]. " "Use %% to select a percentage of the eigenvalues instead " "[Default uses Gauss-Jordan elimination to solve the linear system]."); + GMT_Usage (API, 3, "+c Valid for 2-D gridding only and will create a series of intermediate " + "grids for each eigenvalue holding the cumulative result. Requires -G with a valid filename " + "and extension and we will insert _cum_### before the extension."); GMT_Usage (API, 3, "+f Save the eigenvalues to ."); - GMT_Usage (API, 3, "+m Valid for 2-D gridding only and will create a series of intermediate " - "grids for each eigenvalue holding the incremental result. Requires -G with a filename " - "template containing an integer C formatting specifier (e.g., %%d)."); - GMT_Usage (API, 3, "+M Same, but for cumulative results."); + GMT_Usage (API, 3, "+i As +c but save incremental results, inserting _inc_### before the extension."); gmt_grdcube_info_syntax (API->GMT, 'D'); GMT_Usage (API, 1, "\n-E[]"); GMT_Usage (API, -2, "Evaluate solution at input locations and report misfit statistics. " @@ -458,7 +458,11 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM k = (Ctrl->C.mode) ? 1 : 0; if (gmt_get_modifier (opt->arg, 'f', p)) Ctrl->C.file = strdup (p); - if (gmt_get_modifier (opt->arg, 'm', p)) + if (gmt_get_modifier (opt->arg, 'i', p)) + Ctrl->C.movie |= GREENSPLINE_INC_MOVIE; + if (gmt_get_modifier (opt->arg, 'c', p)) + Ctrl->C.movie |= GREENSPLINE_CUM_MOVIE; + if (gmt_get_modifier (opt->arg, 'm', p)) /* Deprecated for +i only and no +c allowed */ Ctrl->C.movie = GREENSPLINE_INC_MOVIE; else if (gmt_get_modifier (opt->arg, 'M', p)) Ctrl->C.movie = GREENSPLINE_CUM_MOVIE; @@ -468,7 +472,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM Ctrl->C.file = strdup (p); } else { - GMT_Report (API, GMT_MSG_ERROR, "Option -C: Expected -C[n][+]\n"); + GMT_Report (API, GMT_MSG_ERROR, "Option -C: Expected -C[n][+c][+f][+i]\n"); n_errors++; } } @@ -702,7 +706,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && Ctrl->A.mode > 5, "Option -A: format must be in 0-5 range\n"); n_errors += gmt_M_check_condition (GMT, !(GMT->common.R.active[RSET] || Ctrl->N.active || Ctrl->T.active), "No output locations specified (use either [-R -I], -N, or -T)\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->R3.mode && dimension != 2, "The -R or -T option only applies to 2-D gridding\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C +m|M modifier only applies to 2-D gridding\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C +m|M modifiers only apply to 2-D gridding\n"); #ifdef DEBUG n_errors += gmt_M_check_condition (GMT, !TEST && !Ctrl->N.active && Ctrl->R3.dimension != dimension, "The -R and -Z options disagree on the dimension\n"); #else @@ -720,8 +724,6 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->N.file, "Option -N: Must specify node file name\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->N.file && gmt_access (GMT, Ctrl->N.file, R_OK), "Option -N: Cannot read file %s!\n", Ctrl->N.file); n_errors += gmt_M_check_condition (GMT, (Ctrl->I.active + GMT->common.R.active[RSET]) == 1 && dimension == 2, "Must specify -R, -I, [-r], -G for gridding\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && strchr (Ctrl->G.file, '%') == NULL, "Must give a filename template via -G when -C..+m|M is selected\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C..+m|M modifier is only available for 2-D gridding\n"); n_errors += gmt_M_check_condition (GMT, dimension == 3 && Ctrl->G.active && API->external && strchr (Ctrl->G.file, '%'), "Option -G: Cannot contain format-specifiers when not used on the command line\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); @@ -1449,6 +1451,17 @@ GMT_LOCAL void greenspline_dump_system (double *A, double *b, uint64_t nm, char fprintf (stderr, "\t||\t%12.6f\n", b[row]); } } + +GMT_LOCAL void greenspline_set_filename (char *name, unsigned int k, unsigned int width, unsigned int mode, char *file) { + /* Turn name, k, mode into a filename, e.g. for "solution.grd", 33, inc will give solution_inc_033.grd */ + unsigned int s = strlen (name) - 1; + static char *type[3] = {"", "inc", "cum"}; + while (name[s] != '.') s--; /* Wind backwards to start of extension */ + name[s] = '\0'; /* Temporarily chop off extension */ + sprintf (file, "%s_%s_%*.*d.%s", name, type[mode], width, width, k, &name[s+1]); + name[s] = '.'; /* Restore original name */ +} + #define bailout(code) {gmt_M_free_options (mode); return (code);} #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} @@ -2448,17 +2461,18 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } /* Else we are writing a grid or cube */ if (Ctrl->C.movie) { /* 2-D only: Write out grid after adding contribution for each eigenvalue separately */ - gmt_grdfloat *tmp = NULL; + unsigned int width = urint (floor (log10 ((double)nm))) + 1; /* Width of maximum integer needed */ + gmt_grdfloat *current = NULL, *previous = NULL; static char *mkind[3] = {"", "Incremental", "Cumulative"}; char file[PATH_MAX] = {""}; - if (Ctrl->C.movie == GREENSPLINE_INC_MOVIE) tmp = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); + if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) { + current = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); + previous = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); + } gmt_grd_init (GMT, Out->header, options, true); for (k = 0; k < nm; k++) { - GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate %s spline for eigenvalue # %d\n", mkind[Ctrl->C.movie], (int)k+1); - snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s). %s contribution for eigenvalue # %d", method[Ctrl->S.mode], Ctrl->S.arg, mkind[Ctrl->C.movie], (int)k+1); - if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) - Return (API->error); /* Update solution for k eigenvalues only */ + GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)k+1); gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */ (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)k, 1); for (row = 0; row < Grid->header->n_rows; row++) { @@ -2482,16 +2496,39 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { Out->data[ij] = (gmt_grdfloat)V[GMT_Z]; } } - sprintf (file, Ctrl->G.file, (int)k+1); - if (Ctrl->C.movie == GREENSPLINE_INC_MOVIE) { - gmt_M_grd_loop (GMT, Out, row, col, ij) Out->data[ij] -= tmp[ij]; /* Incremental improvement since last time */ - gmt_M_grd_loop (GMT, Out, row, col, ij) tmp[ij] += Out->data[ij]; /* Current solution */ + if (Ctrl->C.movie) gmt_M_memcpy (current, Out->data, Out->header->size, gmt_grdfloat); /* Save current solution */ + + if (Ctrl->C.movie & GREENSPLINE_CUM_MOVIE) { /* Write out the cumulative solution first */ + if (strchr (Ctrl->G.file, '%')) /* Gave a template, use it to write one of the two types of grids */ + sprintf (file, Ctrl->G.file, (int)k+1); + else /* Create the appropriate cumulative gridfile name from static file name */ + greenspline_set_filename (Ctrl->G.file, k, width, GREENSPLINE_CUM_MOVIE, file); + snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s). %s contribution for eigenvalue # %d", method[Ctrl->S.mode], Ctrl->S.arg, mkind[GREENSPLINE_CUM_MOVIE], (int)k+1); + if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) + Return (API->error); /* Update solution for k eigenvalues only */ + if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Out) != GMT_NOERROR) { + Return (API->error); + } } - if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Out) != GMT_NOERROR) { - Return (API->error); + if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) { /* Want to write out incremental solution due to this eigenvalue */ + gmt_M_grd_loop (GMT, Out, row, col, ij) Out->data[ij] = current[ij] - previous[ij]; /* Incremental improvement since last time */ + gmt_M_memcpy (previous, current, Out->header->size, gmt_grdfloat); /* Save current solution which will be previous for next eigenvalue */ + if (strchr (Ctrl->G.file, '%')) /* Gave a template, use it to write one of the two types of grids */ + sprintf (file, Ctrl->G.file, (int)k+1); + else /* Create the appropriate cumulative gridfile name from static file name */ + greenspline_set_filename (Ctrl->G.file, k, width, GREENSPLINE_INC_MOVIE, file); + snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s). %s contribution for eigenvalue # %d", method[Ctrl->S.mode], Ctrl->S.arg, mkind[GREENSPLINE_INC_MOVIE], (int)k+1); + if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) + Return (API->error); /* Update solution for k eigenvalues only */ + if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Out) != GMT_NOERROR) { + Return (API->error); + } } } - if (Ctrl->C.movie == GREENSPLINE_INC_MOVIE) gmt_M_free_aligned (GMT, tmp); /* Free original column-oriented grid */ + if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) { /* Free temporary arrays */ + gmt_M_free_aligned (GMT, current); + gmt_M_free_aligned (GMT, previous); + } gmt_M_free (GMT, A); gmt_M_free (GMT, s); gmt_M_free (GMT, v); From 3567580f8ece6570a414b400130f8602e169379e Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 29 Aug 2021 15:15:40 -1000 Subject: [PATCH 2/7] typos --- doc/rst/source/greenspline.rst | 2 +- src/greenspline.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/rst/source/greenspline.rst b/doc/rst/source/greenspline.rst index a47ea95d0a5..84d1c906ce7 100644 --- a/doc/rst/source/greenspline.rst +++ b/doc/rst/source/greenspline.rst @@ -153,7 +153,7 @@ Optional Arguments execution will stop after saving the eigenvalues, i.e., no surface output is produced. Specify **-Cn** to retain only the *value* largest eigenvalues; append % if *value* is the percentage of eigenvalues - to use instead. The two last modifiers (**+c** and **i**) are only + to use instead. The two other modifiers (**+c** and **i**) are only available for 2-D gridding and can be used to write intermediate grids, one per eigenvalue, and thus require a file name with a suitable extension to be given via **-G** (we automatically insert "_cum_###" or "_inc_###" diff --git a/src/greenspline.c b/src/greenspline.c index 0ee104dc4cc..7ce99d78bbb 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -706,7 +706,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && Ctrl->A.mode > 5, "Option -A: format must be in 0-5 range\n"); n_errors += gmt_M_check_condition (GMT, !(GMT->common.R.active[RSET] || Ctrl->N.active || Ctrl->T.active), "No output locations specified (use either [-R -I], -N, or -T)\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->R3.mode && dimension != 2, "The -R or -T option only applies to 2-D gridding\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C +m|M modifiers only apply to 2-D gridding\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C +c+i modifiers only apply to 2-D gridding\n"); #ifdef DEBUG n_errors += gmt_M_check_condition (GMT, !TEST && !Ctrl->N.active && Ctrl->R3.dimension != dimension, "The -R and -Z options disagree on the dimension\n"); #else From 64a1a813e5a44f699d7bcd0e95ce82b3ac1a0d7c Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 29 Aug 2021 15:20:53 -1000 Subject: [PATCH 3/7] Update greenspline.c --- src/greenspline.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/greenspline.c b/src/greenspline.c index 7ce99d78bbb..b079a0c6b0b 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -707,6 +707,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition (GMT, !(GMT->common.R.active[RSET] || Ctrl->N.active || Ctrl->T.active), "No output locations specified (use either [-R -I], -N, or -T)\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->R3.mode && dimension != 2, "The -R or -T option only applies to 2-D gridding\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && dimension != 2, "The -C +c+i modifiers only apply to 2-D gridding\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.movie && strchr (Ctrl->G.file, '%') == NULL && strchr (Ctrl->G.file, '.') == NULL, "Option -G: When -C +i|c is used your grid file must have an extension\n"); #ifdef DEBUG n_errors += gmt_M_check_condition (GMT, !TEST && !Ctrl->N.active && Ctrl->R3.dimension != dimension, "The -R and -Z options disagree on the dimension\n"); #else @@ -1453,7 +1454,8 @@ GMT_LOCAL void greenspline_dump_system (double *A, double *b, uint64_t nm, char } GMT_LOCAL void greenspline_set_filename (char *name, unsigned int k, unsigned int width, unsigned int mode, char *file) { - /* Turn name, k, mode into a filename, e.g. for "solution.grd", 33, inc will give solution_inc_033.grd */ + /* Turn name, eigenvalue number k, precision width and mode into a filename, e.g., + * ("solution.grd", 33, 3, GREENSPLINE_INC_MOVIE, file) will give solution_inc_033.grd */ unsigned int s = strlen (name) - 1; static char *type[3] = {"", "inc", "cum"}; while (name[s] != '.') s--; /* Wind backwards to start of extension */ From 683c8e4b98b496629ea2a2b3f6f724b26735dcb8 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 30 Aug 2021 15:45:56 -1000 Subject: [PATCH 4/7] Fix alloc and freeing of temp grids --- src/greenspline.c | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/greenspline.c b/src/greenspline.c index b079a0c6b0b..9f47e7059f4 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -2467,9 +2467,10 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { gmt_grdfloat *current = NULL, *previous = NULL; static char *mkind[3] = {"", "Incremental", "Cumulative"}; char file[PATH_MAX] = {""}; - if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) { + if (Ctrl->C.movie) { current = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); - previous = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); + if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) + previous = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); } gmt_grd_init (GMT, Out->header, options, true); @@ -2527,9 +2528,10 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } } } - if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) { /* Free temporary arrays */ + if (Ctrl->C.movie) { /* Free temporary arrays */ gmt_M_free_aligned (GMT, current); - gmt_M_free_aligned (GMT, previous); + if (Ctrl->C.movie & GREENSPLINE_INC_MOVIE) + gmt_M_free_aligned (GMT, previous); } gmt_M_free (GMT, A); gmt_M_free (GMT, s); From 0c3956481e98a68e3c1d2560d3d82db5f97afe9c Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 30 Aug 2021 17:16:01 -1000 Subject: [PATCH 5/7] Address final issues --- doc/rst/source/greenspline.rst | 8 ++--- src/gmt_vector.c | 3 +- src/greenspline.c | 65 ++++++++++++++++++++++------------ 3 files changed, 49 insertions(+), 27 deletions(-) diff --git a/doc/rst/source/greenspline.rst b/doc/rst/source/greenspline.rst index 84d1c906ce7..0084f2351cf 100644 --- a/doc/rst/source/greenspline.rst +++ b/doc/rst/source/greenspline.rst @@ -15,7 +15,7 @@ Synopsis **gmt greenspline** [ *table* ] |-G|\ *grdfile* [ |-A|\ *gradfile*\ **+f**\ **1**\|\ **2**\|\ **3**\|\ **4**\|\ **5** ] -[ |-C|\ [**n**]\ *value*\ [%][**+c**][**+f**\ *file*][**+i**] ] +[ |-C|\ [[**n**]\ *value*\ [%]][**+c**][**+f**\ *file*][**+i**] ] [ |SYN_OPT-D3| ] [ |-E|\ [*misfitfile*] ] [ |-I|\ *xinc*\ [/*yinc*\ [/*zinc*]] ] @@ -142,7 +142,7 @@ Optional Arguments .. _-C: -**-C**\ [**n**]\ *value*\ [%][**+c**][**+f**\ *file*][**+i**] +**-C**\ [[**n**]\ *value*\ [%]][**+c**][**+f**\ *file*][**+i**] Find an approximate surface fit: Solve the linear system for the spline coefficients by SVD and eliminate the contribution from all eigenvalues whose ratio to the largest eigenvalue is less than *value* @@ -151,8 +151,8 @@ Optional Arguments eigenvalues to the specified file for further analysis. If a negative *value* is given then **+f**\ *file* is required and execution will stop after saving the eigenvalues, i.e., no surface - output is produced. Specify **-Cn** to retain only the *value* largest - eigenvalues; append % if *value* is the percentage of eigenvalues + output is produced. Specify **-Cn**\ *value* to retain only the *value* largest + eigenvalues; append % if *value* is the *percentage* of eigenvalues to use instead. The two other modifiers (**+c** and **i**) are only available for 2-D gridding and can be used to write intermediate grids, one per eigenvalue, and thus require a file name with a suitable extension diff --git a/src/gmt_vector.c b/src/gmt_vector.c index eac50e54cf7..96d035a7024 100644 --- a/src/gmt_vector.c +++ b/src/gmt_vector.c @@ -1099,6 +1099,7 @@ int gmt_svdcmp (struct GMT_CTRL *GMT, double *a, unsigned int m_in, unsigned int int gmt_solve_svd (struct GMT_CTRL *GMT, double *u, unsigned int m, unsigned int nu, double *v, double *w, double *b, unsigned int k, double *x, double cutoff, unsigned int mode) { /* Mode = 0: Use all singular values s_j for which s_j/s_0 > cutoff [0 = all] * mode = 1: Use the first cutoff singular values only. If cutoff is < 1 we assume this is the fraction of eigenvalues we want. + * We return the number of eigenvalues used. */ double w_abs, sing_max; int i, j, n_use = 0, n = (int)nu; /* Because OpenMP cannot handle unsigned loop variables */ @@ -1164,7 +1165,7 @@ int gmt_solve_svd (struct GMT_CTRL *GMT, double *u, unsigned int m, unsigned int } if (mode == 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, - "gmt_solve_svd: Ratio limit %g ratained %d singular values\n", cutoff, n_use); + "gmt_solve_svd: Ratio limit %g retained %d singular values\n", cutoff, n_use); if (mode == 2) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "gmt_solve_svd: Selected first %d singular values\n", n_use); diff --git a/src/greenspline.c b/src/greenspline.c index 9f47e7059f4..8e7afc24ac4 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -227,7 +227,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *C) { /* De static int usage (struct GMTAPI_CTRL *API, int level) { const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); - GMT_Usage (API, 0, "usage: %s [
] -G [-A+f] [-C[n][%%][+c][+f][+i]] " + GMT_Usage (API, 0, "usage: %s [
] -G [-A+f] [-C[[n][%%]][+c][+f][+i]] " "[-D] [-E[]] [-I[/[/]]] [-L] [-N] [-Q] " "[-R//[//]]] [-Sc|l|t|r|p|q[]] [-T] " "[%s] [-W[w]] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s [%s]%s[%s] [%s]\n", @@ -257,11 +257,11 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "4: X, Vcomponents."); GMT_Usage (API, 3, "5: X, Vunit-vector, Vmagnitude."); GMT_Usage (API, -2, "Here, X = (x, y[, z]) is the position vector, V = (Vx, Vy[, Vz]) is the gradient vector."); - GMT_Usage (API, 1, "\n-C[n][%%][+c][+f][+i]"); + GMT_Usage (API, 1, "\n-C[[n][%%]][+c][+f][+i]"); GMT_Usage (API, -2, "Solve by SVD and eliminate eigenvalues whose ratio to largest eigenvalue is less than [0]. " "A negative cutoff will stop execution after reporting the eigenvalues. " "Use -Cn to select only the largest eigenvalues [all]. " - "Use %% to select a percentage of the eigenvalues instead " + "Use %% to select a percentage of the eigenvalues instead [100] " "[Default uses Gauss-Jordan elimination to solve the linear system]."); GMT_Usage (API, 3, "+c Valid for 2-D gridding only and will create a series of intermediate " "grids for each eigenvalue holding the cumulative result. Requires -G with a valid filename " @@ -2227,7 +2227,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { /* Sort singular values into ascending order */ gmt_sort_array (GMT, eig, nm, GMT_DOUBLE); for (i = 0, j = nm-1; i < nm; i++, j--) { - E->table[0]->segment[0]->data[GMT_X][i] = i + 1.0; /* Let 1 be x-value of the first eigenvalue */ + E->table[0]->segment[0]->data[GMT_X][i] = i; /* Let 0 be x-value of the first eigenvalue */ E->table[0]->segment[0]->data[GMT_Y][i] = eig[j]; } E->table[0]->segment[0]->n_rows = nm; @@ -2463,7 +2463,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } /* Else we are writing a grid or cube */ if (Ctrl->C.movie) { /* 2-D only: Write out grid after adding contribution for each eigenvalue separately */ - unsigned int width = urint (floor (log10 ((double)nm))) + 1; /* Width of maximum integer needed */ + unsigned int width = urint (floor (log10 ((double)n_use))) + 1; /* Width of maximum integer needed */ gmt_grdfloat *current = NULL, *previous = NULL; static char *mkind[3] = {"", "Incremental", "Cumulative"}; char file[PATH_MAX] = {""}; @@ -2474,31 +2474,52 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } gmt_grd_init (GMT, Out->header, options, true); - for (k = 0; k < nm; k++) { - GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)k+1); + for (k = 0; k < n_use; k++) { + GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)k); gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */ (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)k, 1); - for (row = 0; row < Grid->header->n_rows; row++) { - V[GMT_Y] = yp[row]; - for (col = 0; col < Grid->header->n_columns; col++) { - ij = gmt_M_ijp (Grid->header, row, col); - if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */ - V[GMT_X] = xp[col]; - /* Here, V holds the current output coordinates */ - for (p = 0, wp = 0.0; p < nm; p++) { - r = greenspline_get_radius (GMT, V, X[p], 2U); - if (Ctrl->Q.active) { + if (Ctrl->Q.active) { /* Derivatives of solution */ +#ifdef _OPENMP +#pragma omp parallel for private(row,V,col,ij,p,wp,r,C,part) shared(Grid,yp,xp,n_use,GMT,Ctrl,X,G,par,Lz,alpha,Out,normalize,norm) +#endif + for (row = 0; row < Grid->header->n_rows; row++) { + V[GMT_Y] = yp[row]; + for (col = 0; col < Grid->header->n_columns; col++) { + ij = gmt_M_ijp (Grid->header, row, col); + if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */ + V[GMT_X] = xp[col]; + /* Here, V holds the current output coordinates */ + for (p = 0, wp = 0.0; p < n_use; p++) { + r = greenspline_get_radius (GMT, V, X[p], 2U); C = greenspline_get_dircosine (GMT, Ctrl->Q.dir, V, X[p], 2U, false); part = dGdr (GMT, r, par, Lz) * C; + wp += alpha[p] * part; } - else - part = G (GMT, r, par, Lz); - wp += alpha[p] * part; + V[GMT_Z] = greenspline_undo_normalization (V, wp, normalize, norm, 2U); + Out->data[ij] = (gmt_grdfloat)V[GMT_Z]; } - V[GMT_Z] = greenspline_undo_normalization (V, wp, normalize, norm, 2U); - Out->data[ij] = (gmt_grdfloat)V[GMT_Z]; } } + else { +#ifdef _OPENMP +#pragma omp parallel for private(row,V,col,ij,p,wp,r,part) shared(Grid,yp,xp,n_use,GMT,X,G,par,Lz,alpha,Out,normalize,norm) +#endif + for (row = 0; row < Grid->header->n_rows; row++) { + V[GMT_Y] = yp[row]; + for (col = 0; col < Grid->header->n_columns; col++) { + ij = gmt_M_ijp (Grid->header, row, col); + if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */ + V[GMT_X] = xp[col]; + /* Here, V holds the current output coordinates */ + for (p = 0, wp = 0.0; p < n_use; p++) { + r = greenspline_get_radius (GMT, V, X[p], 2U); + part = G (GMT, r, par, Lz); + wp += alpha[p] * part; + } + Out->data[ij] = (gmt_grdfloat)greenspline_undo_normalization (V, wp, normalize, norm, 2U); + } + } /* End of row-loop [OpenMP] */ + } if (Ctrl->C.movie) gmt_M_memcpy (current, Out->data, Out->header->size, gmt_grdfloat); /* Save current solution */ if (Ctrl->C.movie & GREENSPLINE_CUM_MOVIE) { /* Write out the cumulative solution first */ From 529b1e656ac10fa741c949fc8747f7924c396135 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 31 Aug 2021 11:39:10 -1000 Subject: [PATCH 6/7] Add comments regarding unsorted eigenvalues vector --- src/gmt_vector.c | 2 +- src/greenspline.c | 15 ++++++++++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/gmt_vector.c b/src/gmt_vector.c index 96d035a7024..583620be5aa 100644 --- a/src/gmt_vector.c +++ b/src/gmt_vector.c @@ -1128,7 +1128,7 @@ int gmt_solve_svd (struct GMT_CTRL *GMT, double *u, unsigned int m, unsigned int /* mode = 1: Find the m largest singular values, with m = cutoff (if <1 it is the fraction of values). * Either case requires sorted singular values so we need to do some work first. * It also assumes that the matrix passed is a squared normal equation kind of matrix - * so that the singular values are the individual variace contributions. */ + * so that the singular values are the individual variance contributions. */ struct GMT_SINGULAR_VALUE { double value; unsigned int order; diff --git a/src/greenspline.c b/src/greenspline.c index 18193294e0a..1d689312a28 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -2505,6 +2505,9 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } /* Else we are writing a grid or cube */ if (Ctrl->C.movie) { /* 2-D only: Write out grid after adding contribution for each eigenvalue separately */ + /* Note: Because the SVD decomposition is not sorting the vectors from largest to smallest eigenvalue the + * gmt_solve_svd sets to zero those we don't want but we must still loop over its full length to ensure we + * include the eigenvalues we want. */ unsigned int width = urint (floor (log10 ((double)n_use))) + 1; /* Width of maximum integer needed */ gmt_grdfloat *current = NULL, *previous = NULL; static char *mkind[3] = {"", "Incremental", "Cumulative"}; @@ -2516,13 +2519,13 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } gmt_grd_init (GMT, Out->header, options, true); - for (k = 0; k < n_use; k++) { + for (k = 0; k < n_use; k++) { /* Only loop over the first n_use eigenvalues (if restricted) */ GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)k); gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */ (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)k, 1); if (Ctrl->Q.active) { /* Derivatives of solution */ #ifdef _OPENMP -#pragma omp parallel for private(row,V,col,ij,p,wp,r,C,part) shared(Grid,yp,xp,n_use,GMT,Ctrl,X,G,par,Lz,alpha,Out,normalize,norm) +#pragma omp parallel for private(row,V,col,ij,p,wp,r,C,part) shared(Grid,yp,xp,nm,GMT,Ctrl,X,G,par,Lz,alpha,Out,normalize,norm) #endif for (row = 0; row < Grid->header->n_rows; row++) { V[GMT_Y] = yp[row]; @@ -2531,7 +2534,8 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */ V[GMT_X] = xp[col]; /* Here, V holds the current output coordinates */ - for (p = 0, wp = 0.0; p < n_use; p++) { + for (p = 0, wp = 0.0; p < nm; p++) { + if (alpha[p] == 0.0) continue; /* Note: The alpha's are not sorted so must loop over all then skip the zeros */ r = greenspline_get_radius (GMT, V, X[p], 2U); C = greenspline_get_dircosine (GMT, Ctrl->Q.dir, V, X[p], 2U, false); part = dGdr (GMT, r, par, Lz) * C; @@ -2544,7 +2548,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } else { #ifdef _OPENMP -#pragma omp parallel for private(row,V,col,ij,p,wp,r,part) shared(Grid,yp,xp,n_use,GMT,X,G,par,Lz,alpha,Out,normalize,norm) +#pragma omp parallel for private(row,V,col,ij,p,wp,r,part) shared(Grid,yp,xp,nm,GMT,X,G,par,Lz,alpha,Out,normalize,norm) #endif for (row = 0; row < Grid->header->n_rows; row++) { V[GMT_Y] = yp[row]; @@ -2553,7 +2557,8 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */ V[GMT_X] = xp[col]; /* Here, V holds the current output coordinates */ - for (p = 0, wp = 0.0; p < n_use; p++) { + for (p = 0, wp = 0.0; p < nm; p++) { /* Note: The alpha's are not sorted so must loop over all then skip the zeros */ + if (alpha[p] == 0.0) continue; /* Note: The alpha's are not sorted so must loop over all then skip the zeros */ r = greenspline_get_radius (GMT, V, X[p], 2U); part = G (GMT, r, par, Lz); wp += alpha[p] * part; From f0f851f61d12f145b7b72297fad6e7d3638fb0fc Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 31 Aug 2021 12:02:50 -1000 Subject: [PATCH 7/7] Fix sign loop vars for Windows OpenMP --- src/greenspline.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/greenspline.c b/src/greenspline.c index 1d689312a28..0b72c7bd3a9 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -2509,6 +2509,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { * gmt_solve_svd sets to zero those we don't want but we must still loop over its full length to ensure we * include the eigenvalues we want. */ unsigned int width = urint (floor (log10 ((double)n_use))) + 1; /* Width of maximum integer needed */ + int64_t col, row, p; /* On Windows, the 'for' index variables must be signed, so redefine these 3 inside this block only */ gmt_grdfloat *current = NULL, *previous = NULL; static char *mkind[3] = {"", "Incremental", "Cumulative"}; char file[PATH_MAX] = {""};