gfortran bug?

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.

As far as I can tell, your fortran code is in fixed format which is not immediately obvious.
I have put the code in free format and  have reformatted the code with more indentation (1 space of indent is really not enough).
See below after my signature.
It  now compiles with no errors.

I don't see anything immediately wrong.
It would be nice to have a small test program in fortran to check the code.

Berend

     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

BTW. in       do i = 1, nl(1), 1  and    do j = 1, ndim(3), 1
the steps 1 are not necessary: 1 is the default