From bhh at xs4all.nl  Thu Feb  9 19:59:58 2017
From: bhh at xs4all.nl (Berend Hasselman)
Date: Thu, 9 Feb 2017 19:59:58 +0100
Subject: [Rd] Ancient C /Fortran code linpack error
In-Reply-To: <22684.40017.122891.70235@stat.math.ethz.ch>
References: <7240dae6-2dc9-88ec-c9c8-f552fc1c2dd5@umu.se>
	<7B0E2308-8715-4BAA-A9F3-4FFDADF990A2@xs4all.nl>
	<22684.40017.122891.70235@stat.math.ethz.ch>
Message-ID: <C6AEB448-97C1-4AC1-B19B-DC957B9B074D@xs4all.nl>


> 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



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