[R-SIG-Mac] gfortran bug?

[Andreas Noack Jensen]

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
-- 

Andreas Noack Jensen        Øster Farimagsgade 5, bygning 26
Ph.d.-stipendiat            1353 København K
Økonomisk Institut          Danmark
Københavns Universitet      Tel. (+1)646-643-0985