[Rd] Ancient C /Fortran code linpack error

[Berend Hasselman]

> On 9 Feb 2017, at 16:00, Göran Broström <goran.brostrom at umu.se> wrote:
> 
> In my package 'glmmML' I'm using old C code and linpack in the optimizing procedure. Specifically, one part of the code looks like this:
> 
>    F77_CALL(dpoco)(*hessian, &bdim, &bdim, &rcond, work, info);
>    if (*info == 0){
>        F77_CALL(dpodi)(*hessian, &bdim, &bdim, det, &job);
>        ........
> 
> This usually works OK, but with an ill-conditioned data set (from a user of glmmML) it happened that the hessian was all nan. However, dpoco returned *info = 0 (no error!) and then the call to dpodi hanged R!
> 
> I googled for C and nan and found a work-around: Change 'if ...' to
> 
>   if (*info == 0 & (hessian[0][0] == hessian[0][0])){
> 
> which works as a test of hessian[0][0] (not) being NaN.
> 
> I'm using the .C interface for calling C code.
> 
> Any thoughts on how to best handle the situation? Is this a bug in dpoco? Is there a simple way to test for any NaNs in a vector?

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];
        }

There may be a more compact way with a macro in the R headers.
I feel that If other code can't handle non-finite values: then test.

Berend Hasselman