From bhh at xs4all.nl  Fri Feb 10 16:13:53 2017
From: bhh at xs4all.nl (Berend Hasselman)
Date: Fri, 10 Feb 2017 16:13:53 +0100
Subject: [Rd] Ancient C /Fortran code linpack error
In-Reply-To: <8b7e4792-f306-87e6-3cb8-7d1611d1a640@umu.se>
References: <7240dae6-2dc9-88ec-c9c8-f552fc1c2dd5@umu.se>
	<7B0E2308-8715-4BAA-A9F3-4FFDADF990A2@xs4all.nl>
	<22684.40017.122891.70235@stat.math.ethz.ch>
	<8b7e4792-f306-87e6-3cb8-7d1611d1a640@umu.se>
Message-ID: <810D55DC-9BA8-4F86-9DFD-C32C7E8697FA@xs4all.nl>


> On 10 Feb 2017, at 14:53, G?ran Brostr?m <goran.brostrom at umu.se> wrote:
> 
> 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.)
> 

Only if you have lda and n equal (which you indeed have; but still worth mentioning) when calling dpoco.

Berend



[Rd] Ancient C /Fortran code linpack error
Berend Hasselman