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
