[Rd] DUP=T/F

[Paul Gilbert]

In some tests I am getting a difference in the 7th or 8th significant digit depending
on whether I set DUP=T or DUP=F in the .Fortran call below. While this is not a huge
error it is bigger than I would expect and may be a symptom of another problem. Any
comments or suggestions? (R1.2.2 on Solaris.)

Paul Gilbert
_______

genD.ARMA <- function(model, data, d=0.01, eps=1e-4, r=6, warn=F){
   n <-length(c(output.data(data))) # this has the same length as the residual
   sampleT <-periods(data)
   ps <-dim(model$A)[3]
      ms <-dim(model$C)[3]
      ns <-ps  # this could be 1 except z0 is used for TREND and F, Q and R for
scratch
     loc   <- match(model$location,       c("A","B","C","t"))
     cloc  <- match(model$const.location, c("A","B","C","t"))
   if(is.null(ms))
     {ms<-0
      C<-array(0,c(1,ns,1)) # can't call compiled with 0 length arrays
      u <- matrix(0,sampleT,1)
      G <- matrix(0,ns,1)
     }
   else
     {C <-model$C
      u <- input.data(data)
      G <- matrix(0,ns,ms)
     }
   x <-model$parms
   h0 <- abs(d*x)+eps*(x==0.0)
   D <- matrix(1e20, n,(length(x)*(length(x) + 3))/2)
   Daprox <-  matrix(0,n,r)
   Hdiag  <-  matrix(0,n,length(x))
   Haprox <-  matrix(0,n,r)
   storage.mode(D)     <-"double"
   storage.mode(Daprox)<-"double"
   storage.mode(Haprox)<-"double"
   storage.mode(Hdiag )<-"double"
   D <-.Fortran("gend",
            D=D,
            as.integer(is.ARMA(model)),
            p=as.integer(length(x)),
            x=as.double(x),
            as.double(h0),
            as.integer(n),    #6
            as.integer((length(x)*(length(x)+ 3))/2), #cols of D
            f0=double(n),
            r=as.integer(r),
            #                       work space for GEND
            Haprox=Haprox,     #10
            Hdiag=Hdiag,
            Daprox=Daprox,
            x=double(length(x)),
            delta=double(length(x)),
            f1=double(n),
            f2=double(n),
            #                   work space for ARMAp / KFp
         #     cov=matrix(1e20,ps,ps),
            # pred is f0,f1,f2 passed above
            as.integer(ms),     #  input dimension m  #17
            as.integer(ps),     # output dimension p
            as.integer(sampleT),
            as.integer(periods(data)),
            u=as.double(u),
            y=as.double(output.data(data)),
            as.integer(loc),   #as.character(model$location), #23
            as.integer(model$i),
            as.integer(model$j),
            as.integer(length(model$const)),
            const=as.double(model$const),
            as.integer(cloc),  #as.character(model$const.location),
            as.integer(model$const.i ),
            as.integer(model$const.j),
        #  for ARMA models:
            as.integer(model$l),       #31
            as.integer(model$const.l),
            as.integer( dim(model$A)[1]),#1+order of A
            as.integer( dim(model$B)[1]),#1+order of B
            as.integer( dim(C)[1]),#1+order of C
            A=model$A,
            B=model$B,
            C=C,
        #  for state space models:
            as.integer(ns),  # state dim.  #39
        #    state=as.double(matrix(0,sampleT,ns)),
        #    track=as.double(array(0,c(sampleT,ns,ns))),
            z0=as.double(rep(0,ps)), # note: this is TREND for ARMA
            P0=as.double(diag(0,ps)),
            F=as.double(matrix(0,ns,ns)),
            G=as.double(G),
            H=as.double(matrix(0,ps,ns)),
            K=as.double(matrix(0,ns,ps)),
            Q=as.double(matrix(0,ns,ns)),
            R=as.double(matrix(0,ps,ps)),
            gain=as.integer(F),   #48
            DUP=T)[c("D","p","f0", "x", "r")]
   D$d   <- d
   D$eps <- eps
   D$f0 <- l(model, data, result="pred") - c(output.data(data))
   invisible(classed(D,"Darray")) #constructor
}


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._