bug in gsw_ipv_vs_fnsquared_ratio() in C library

[richardsc]

Snipped from a developer email:

In the C library, file gsw_oceanographic_toolbox.c, the function gsw_ipv_vs_fnsquared_ratio() is using a value of p_ref that is neither defined in the function argument list nor in the code below it. I’ve put below a snippet. I think the behaviour in this case will be undefined. Some compilers might use 0 for the value of p_ref, but others might use some junk value that was in the system memory at the pointed-to address.
Note that p_ref is not defined before use.

void
gsw_ipv_vs_fnsquared_ratio(double *sa, double *ct, double *p, int nz,
double *ipv_vs_fnsquared_ratio, double *p_mid)
{
int k;

double dsa, sa_mid, dct, ct_mid, dp, p_ref;
double alpha_mid, beta_mid;
double alpha_pref, beta_pref, numerator, denominator;

for (k = 0; k < nz-1; k++) {
   dsa = (sa[k+1] - sa[k]);
   sa_mid = 0.5*(sa[k] + sa[k+1]);
   dct = (ct[k+1] - ct[k]);
   ct_mid = 0.5*(ct[k] + ct[k+1]);
   dp = (p[k+1] - p[k]);
   p_mid[k] = 0.5*(p[k] + p[k+1]);

   alpha_mid = gsw_alpha(sa_mid,ct_mid,p_mid[k]);
   beta_mid = gsw_beta(sa_mid,ct_mid,p_mid[k]);
   alpha_pref = gsw_alpha(sa_mid,ct_mid,p_ref);

[richardsc]

I can confirm that this segfaults even if run from the old teos() method.

> teos('gsw_ipv_vs_fnsquared_ratio', 34, 28, 10)

 *** caught segfault ***
address 0x40, cause 'memory not mapped'

Traceback:
 1: .C("gsw3a", as.character(lib), as.character(name), as.integer(ngood),     as.double(a1[good]), as.double(a2[good]), as.double(a3[good]),     rval = double(ngood), NAOK = TRUE, PACKAGE = "oce")
 2: teos("gsw_ipv_vs_fnsquared_ratio", 34, 28, 10)

[richardsc]

Let's wait to see what the higher ups say, but my vote is that if the fix is obvious and easy we do it for now, and document the change. It's not like we're dependent on the source from somewhere else, since we're directly including a particular snapshot in the package.

I know we don't want to just create another version of the library, but I think that outweighs the wait that would come from waiting for upstream changes to the C library (which in all likelihood is going to change anyway).

[fdelahoyde]

Hi Everyone,

I'm the developer/maintainer of the C library and am in the process of putting
it up on github (I've been using mercurial for my RCS so be patient). I've
indeed noted the bug and it's likely it is also in the F90 library. I believe
the better fix is to change the API to include pref as does the Matlab library.

Frank

            Frank Delahoyde | Phone: 858.534.9562
    Shipboard Technical Support | Computing Resources

Scripps Institution of Oceanography | Nimitz Marine Facility
297 Rosecrans Street |
San Diego, Ca. 92106 | fdelahoyde@ucsd.edu

On 12/30/14, 5:52 PM, Clark Richards wrote:

Let's wait to see what the higher ups say, but my vote is that if the fix is
obvious and easy we do it for now, and document the change. It's not like we're
dependent on the source from somewhere else, since we're directly including a
particular snapshot in the package.

I know we don't want to just create another version of the library, but I think
that outweighs the wait that would come from waiting for upstream changes to the
C library (which in all likelihood is going to change anyway).

—
Reply to this email directly or view it on GitHub
#9 (comment).

[dankelley]

Hi Frank. This is indeed the fix that will go into GSW-R (probably tomorrow morning). The fix looks pretty simple.

It's hard to know who is in charge of the C library. I had thought @PaulMBarker was doing that. Is it you?

[fdelahoyde]

Hi Dan,

Yes, it's me. We couldn't deal with the Fortran (sorry Paul).

Frank

On 12/30/14, 6:17 PM, Dan Kelley wrote:

Hi Frank. This is indeed the fix that will go into GSW-R (probably tomorrow
morning). The fix looks pretty simple.

It's hard to know who is in charge of the C library. I had thought @PaulMBarker
https://github.com/PaulMBarker was doing that. Is it you?

—
Reply to this email directly or view it on GitHub
#9 (comment).

[dankelley]

I've updated the library (C and h files) within GSW-R.

@fdelahoyde in case you want to see the differences in the C library, issue the following two commands in an OS shell:

git diff 7cf2cc e19735e src/gsw_oceanographic_toolbox.c
git diff 7cf2cc e19735e src/gswteos-10.h

you can also see all the differences (including in my R code) with

git diff 7cf2cc e19735e

in case that's of interest. Finally, if you don't like git commands, you can see colour-coded changes in the github webpage

[dankelley]

@richardsc Although this bug has been fixed, we might want to leave the issue open a while in case others don't realize they can read and comment upon closed issues. --Dan.

[dankelley]

Closing, in light of previous "2-day wait to close" argument.

[PaulMBarker]

There is another bug that needs fixing in this code that is how dsa and dct are calculated. Frank alerted me to it.

Paul.

Here is the latest matlab lines:

dSA = SA(Ishallow,:) - SA(Ideep,:);
dCT = CT(Ishallow,:) - CT(Ideep,:);

and here are the fortran lines from my laptop:
dsa = (sa(k+1) - sa(k))
dct = (ct(k+1) - ct(k))

and here are Glenn's lines
forall (i = 1: nz-1)
dsa(i) = sa(i+1) - sa(i)
sa_mid(i) = 0.5d0_(sa(i) + sa(i+1))
dct(i) = ct(i+1) - ct(i)
ct_mid(i) = 0.5d0_(ct(i) + ct(i+1))
dp(i) = p(i+1) - p(i)
p_mid(i) = 0.5d0*(p(i) + p(i+1))
end forall

they are different and this is wrong. The fortran needs changing to match the matlab, it should be
dsa = (sa(k) - sa(k+1))
dct = (ct(k) - ct(k+1))

and in Glenn's version to:
forall (i = 1: nz-1)
dsa(i) = sa(i) - sa(i+1)
sa_mid(i) = 0.5d0_(sa(i) + sa(i+1))
dct(i) = ct(i) - ct(i+1)
ct_mid(i) = 0.5d0_(ct(i) + ct(i+1))
p_mid(i) = 0.5d0*(p(i) + p(i+1))
end forall

Note that in the fortran code the line
dp = (p(k+1) - p(k))
or
dp(i) = p(i+1) - p(i)
is unnecessary and should be deleted.

---

From: Dan Kelley [notifications@github.com]
Sent: Saturday, January 03, 2015 3:48 AM
To: TEOS-10/GSW-R
Cc: Paul Barker
Subject: Re: [GSW-R] bug in gsw_ipv_vs_fnsquared_ratio() in C library (#9)

Closing, in light of previous "2-day wait to close" argument.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/9#issuecomment-68541106.

[dankelley]

This is fixed in the new version of the C library, incorporated into gsw today.