Andreas,
I can't even get your code to compile - can you send the actual files, please?
Also you may want to ask on R-devel since I'm not sure there are many Fortran
specialists reading this mailing list... (Although there is a possibility of
miscompilations in majority the cases the cause is user code. The fact that it
works on another platform is not necessarily a proof that it is not a bug in
the code - it is often a coincidence).
Cheers,
Simon
PS: GNU Fortran 4.2.4 is really 4.2.4 - even if it is used in gcc 4.2.1 Xcode -
that's what you see on the command line.
On Apr 27, 2010, at 8:42 AM, Andreas Noack Jensen wrote:
> 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.jen...@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
>
> _______________________________________________
> R-SIG-Mac mailing list
> R-SIG-Mac@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-mac
>
>
_______________________________________________
R-SIG-Mac mailing list
R-SIG-Mac@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
