[Rd] Ancient C /Fortran code linpack error

[Tomas Kalibera]

On 02/09/2017 05:05 PM, Berend Hasselman 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];
>          }
>
> 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
And if performance was of importance, you could use the trick from 
mayHaveNaNOrInf in array.c (originally from pqR), but be careful to test 
the individual operands of the sum.
mayHaveNaNOrInf does not do the test for performance reasons, but as a 
result it can have false positives.

Rboolean hasNaNOrInf(double *x, R_xlen_t n)
{
     if ((n&1) != 0 && !R_FINITE(x[0]))
         return TRUE;
     for (R_xlen_t i = n&1; i < n; i += 2)
         if (!R_FINITE(x[i]+x[i+1])&& (!R_FINITE(x[i]) || !R_FINITE(x[i+1]))
             return TRUE;
     return FALSE;
}

Tomas
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel