[Rd] Ancient C /Fortran code linpack error

[Göran Broström]

Thanks to all who answered my third question. I learned something, but:

On 2017-02-09 17:44, Martin Maechler wrote:
>
>>> 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];
>>         }
>
> 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];
>      }
>
[...]

isn't this even easier to read (and correct?):

     for (j = 0; j < n*; j++)
          for (i = 0; i < n*; i++){
               if ( !R_FINITE(hessian[i][j]) ) error("blah...")
          }

? In .C land, that is. (And sure, I'm afraid of ±Inf in this context.)

Thanks again, Göran