[Rd] Ancient C /Fortran code linpack error
Berend Hasselman
bhh at xs4all.nl
Thu Feb 9 19:59:58 CET 2017
> On 9 Feb 2017, at 17:44, Martin Maechler <maechler at stat.math.ethz.ch> wrote:
>
[snip] ...
>> You should/could use macro R_FINITE to test each entry of the hessian.
>> In package nleqslv I test for a "correct" jacobian like this in file nleqslv.c in function fcnjac:
>
>> for (j = 0; j < *n; j++)
>> for (i = 0; i < *n; i++) {
>> if( !R_FINITE(REAL(sexp_fjac)[(*n)*j + i]) )
>> error("non-finite value(s) returned by jacobian (row=%d,col=%d)",i+1,j+1);
>> rjac[(*ldr)*j + i] = REAL(sexp_fjac)[(*n)*j + i];
>> }
>
> A minor hint on that: While REAL(.) (or INTEGER(.) ...) is really cheap in
> the R sources themselves, that is not the case in package code.
>
> Hence, not only nicer to read but even faster is
>
> double *fj = REAL(sexp_fjac);
> for (j = 0; j < *n; j++)
> for (i = 0; i < *n; i++) {
> if( !R_FINITE(fj[(*n)*j + i]) )
> error("non-finite value(s) returned by jacobian (row=%d,col=%d)",i+1,j+1);
> rjac[(*ldr)*j + i] = fj[(*n)*j + i];
> }
>
Tom, Martin,
Thanks for both of your suggestions. Very helpful.
Berend Hasselman
More information about the R-devel
mailing list