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
_______________________________________________
R-SIG-Mac mailing list
R-SIG-Mac@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
