[R-SIG-Mac] gfortran bug?

[Andreas Noack Jensen]

I have been able to make an example that shows the weird behaviour. For the
Fortran code in last post run

rmfilter    <- function(x, coef)
{
    dyn.load('~/rmfilter.so')
    nl    <- dim(x)
    cdim    <- dim(coef)
    init    <- matrix(0, cdim[3], nl[2])
    ans    <- .Fortran('rmfilter', as.double(x), as.integer(nl),
as.double(coef), as.integer(cdim), as.double(init))[[1]]
    dim(ans)    <- dim(x)
    return(ans)
}

eps <- matrix(rnorm(300), ncol = 3)
tmp <- lapply(seq(5000), function(x) rmfilter(eps, array(diag(3)*0.5,
c(3,3,2))))
sum(!sapply(2:5000, function(x) all.equal(tmp[[1]], tmp[[x]])) == 'TRUE')

The three last runs I have done give

tmp <- lapply(seq(5000), function(x) rmfilter(eps, array(diag(3)*0.5,
c(3,3,2))))
> sum(!sapply(2:5000, function(x) all.equal(tmp[[1]], tmp[[x]])) == 'TRUE')
[1] 639
> tmp <- lapply(seq(5000), function(x) rmfilter(eps, array(diag(3)*0.5,
c(3,3,2))))
> sum(!sapply(2:5000, function(x) all.equal(tmp[[1]], tmp[[x]])) == 'TRUE')
[1] 0
> tmp <- lapply(seq(5000), function(x) rmfilter(eps, array(diag(3)*0.5,
c(3,3,2))))
> sum(!sapply(2:5000, function(x) all.equal(tmp[[1]], tmp[[x]])) == 'TRUE')
[1] 29

I hope you have some ideas.

Best
Andreas


Den 25/04/10 16.22 skrev "Andreas Noack Jensen"
<andreas.noack.jensen at econ.ku.dk> følgende:

> Dear list members
> 
> I have organised some code in a package including a Fortran subroutine for a
> multivariate recursive filter to simulate VAR-porcesses. I have worked with
> R for some time, but I am new to writing packages and coding Fortran so
> first I thought I had made some mistakes in the Fortran code but I have now
> tested the package on Ubuntu 10.4 and Windows XP and the problem is only
> present on Mac.
> 
> On Mac I quite randomly get very weird results from the subroutine. It is
> like the initial values explodes to fx 2e+290 even for identical calls
> without random number generation. Most results are identical to results I
> get from a pure R code but at some calls they explode in a non systematic
> fashion.
> 
> I would have tried to use the gfortran 4.2.4 from r.research.att.com, to
> isolate the error (right now I use their 4.2.3), but it seems like it
> installs 4.2.1 instead and I get the same error still. I have tried to build
> gfortran myself but was unsuccessful in doing so.
> 
> I cannot produce a small example you can run but as an example watch this:
> 
>> for(i in 1:500){
> +     cat(i)
> +     tmp <- I1(civecm:::simulate.I1(fit.postWW2, res =
> fit.postWW2$residuals), 2, 'drift')
> + }
> 1234567891011121314Error in svd(X) : infinite or missing values in 'x'
>> for(i in 1:500){
> +     cat(i)
> +     tmp <- I1(civecm:::simulate.I1(fit.postWW2, res =
> fit.postWW2$residuals), 2, 'drift')
> + }
> 1234Error in svd(X) : infinite or missing values in 'x'
> 
> which is weird because there is no random number generation in play.
> Completely identical calls every time. The error comes from the weird data
> generated in the Fortran subroutine.
> 
> Some information about my system:
> R version 2.11.0 (2010-04-22)
> x86_64-apple-darwin9.8.0
> 
> locale:
> [1] da_DK.UTF-8/da_DK.UTF-8/C/C/da_DK.UTF-8/da_DK.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] civecm_0.1-1   timsac_1.2.1   MASS_7.3-5     foreign_0.8-40
> 
> I use the latest prebuild version of R from r.research.att.com.
> 
> Finally the Fortran and R code for the subroute is:
> 
>       subroutine rmfilter(x, nl, coef, ndim, init)
> 
>       implicit none
>       integer, intent(in) :: nl(2), ndim(3)
>       double precision, intent(inout) :: x(nl(1),nl(2))
>       double precision, intent(in) :: coef(ndim(1),ndim(2),ndim(3)),
>      &init(ndim(3),ndim(nl(2)))
>       double precision :: y(nl(1),nl(2))
>       integer :: i, j
> 
>       y = 0.0d+0 
>       
>       do i = 1, nl(1), 1
>        do j = 1, ndim(3), 1
>         if ( i <= j ) then
>          call dgemm('N', 'T', 1, ndim(1), nl(2), 1.0d+0,
>      &init(ndim(3)-j+1,:), 1, coef(:,:,j), ndim(1), 1.0d+0, y(i,:), 1)
>         else
>          call dgemm('N', 'T', 1, ndim(1), nl(2), 1.0d+0, y(i-j,:), 1,
>      &coef(:,:,j), ndim(1), 1.0d+0, y(i,:), 1)
>         end if
>        end do
>        y(i,:) = y(i,:) + x(i,:)
>       end do
>       x = y
>       return
>       end subroutine rmfilter
> 
> and
> 
> 
> rmfilter    <- function(x, coef, init)
> {
>     if(any(is.na(x))) stop('Series includes NAs')
>     if(any(is.na(coef))) stop('Coefficients includes NAs')
>     if(is.null(dim(x))) stop('Input series has no dimensions')
>     nl        <- dim(x)
>     if(!is.array(coef)) stop('coef must be an array')
>     cdim    <- dim(coef)
>     if(nl[2] != cdim[2] | cdim[1] != cdim[2]) stop('coef has wrong
> dimensions')
>     if(missing(init)) init <- matrix(0, cdim[3], nl[2])
>     else if(any(is.na(init))) stop('Initial values includes NAs')
>     else if(any(dim(init) != c(cdim[3], nl[2]))) stop('init matrix has wrong
> dimensions')
>     ans    <- .Fortran('rmfilter', as.double(x), as.integer(nl),
> as.double(coef), as.integer(cdim), as.double(init))[[1]]
>     dim(ans)    <- dim(x)
>     dimnames(ans)<- dimnames(x)
>     tsp(ans)    <- tsp(hasTsp(x))
>     class(ans)    <- class(x)
>     return(ans)
> }
> 
> Best
> Andreas