Skip to content

Commit

Permalink
further tidyed fortran codes
Browse files Browse the repository at this point in the history
  • Loading branch information
@helske
helske committed Apr 16, 2015
1 parent f752dff commit 585b4b0
Show file tree
Hide file tree
Showing 15 changed files with 781 additions and 686 deletions.
2 changes: 1 addition & 1 deletion R/rstandard.KFS.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#' @export
#' @details For object of class KFS with fully Gaussian observations, several
#' types of standardized residuals can be computed. Also the standardization
#' in for multivariate residuals can be done either by Cholesky decomposition
#' for multivariate residuals can be done either by Cholesky decomposition
#' \eqn{L^{-1}_t residual_t,}{L^(-1)[t]residual[t]} or component-wise
#' \eqn{residual_t/sd(residual_t),}{residual[t]/sd(residual[t])}.
#'
Expand Down
2 changes: 1 addition & 1 deletion man/rstandard.KFS.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ Extract Standardized Residuals from KFS output
\details{
For object of class KFS with fully Gaussian observations, several
types of standardized residuals can be computed. Also the standardization
in for multivariate residuals can be done either by Cholesky decomposition
for multivariate residuals can be done either by Cholesky decomposition
\eqn{L^{-1}_t residual_t,}{L^(-1)[t]residual[t]} or component-wise
\eqn{residual_t/sd(residual_t),}{residual[t]/sd(residual[t])}.

Expand Down
2 changes: 1 addition & 1 deletion src/approxloop.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ subroutine approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p
! compute new estimate of thetahat
rankp2 = rankp
call kfstheta(ytilde, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, tvrqr, a1, p1, p1inf, &
p, n, m, r,tol,rankp2,thetanew,lik)
p, n, m, r,tol,rankp2,thetanew,lik)



Expand Down
4 changes: 2 additions & 2 deletions src/covmeanw.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ subroutine covmeanw(x,w,m,n,k,meanx,covx)
meanx = meanx + x(:,:,i)*w(i)
end do
do i = 1, k
x(:,:,i) = sqrt(w(i))*(x(:,:,i) - meanx)
x(:,:,i) = sqrt(w(i))*(x(:,:,i) - meanx)
end do

do t = 1, n
Expand Down Expand Up @@ -46,7 +46,7 @@ subroutine covmeanwprotect(x,w,m,n,k,meanx,covx)
x2(:,:,i) = sqrt(w(i))*(x(:,:,i) - meanx)
end do

do t = 1, n
do t = 1, n
call dgemm('n','t',m,m,k,1.0d0,x2(:,t,:),m,x2(:,t,:),m,0.0d0,covx(:,:,t),m)
end do

Expand Down
45 changes: 19 additions & 26 deletions src/filtersimfast.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,55 +20,51 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &
double precision, dimension(m) :: arec
double precision, external :: ddot

external daxpy, dgemv
external dgemv

j=0
d=0
if(dt.GT.0) then
!diffuse filtering begins
arec = a1
diffuse: do while(d .LT. (dt-1))
d = d+1
do j=1, p
if(ymiss(d,j).EQ.0) then
vt(j,d) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1) !arec
if (finf(j,d) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j,d)/finf(j,d),kinf(:,j,d),1,arec,1) !a_rec = a_rec + kinf(:,i,t)*vt(:,t)/finf(j,d)
vt(j,d) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1)
if (finf(j,d) .GT. 0.0d0) then
arec = arec + vt(j,d)/finf(j,d)*kinf(:,j,d)
else
if(ft(j,d) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j,d)/ft(j,d),kt(:,j,d),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)/ft(i,t)
if(ft(j,d) .GT. 0.0d0) then
arec = arec + vt(j,d)/ft(j,d)*kt(:,j,d)
end if
end if
end if
end do

call dgemv('n',m,m,1.0d0,tt(:,:,(d-1)*timevar(3)+1),m,arec,1,0.0d0,at(:,d+1),1)
arec = at(:,d+1)

arec = at(:,d+1)
end do diffuse

d = dt
do j=1, jt
if(ymiss(d,j).EQ.0) then
vt(j,d) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1) !arec
vt(j,d) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1)

if (finf(j,d) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j,d)/finf(j,d),kinf(:,j,d),1,arec,1) !a_rec = a_rec + kinf(:,i,t)*vt(:,t)/finf(j,d)
if (finf(j,d) .GT. 0.0d0) then
arec = arec + vt(j,d)/finf(j,d)*kinf(:,j,d)
else
if(ft(j,d) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j,d)/ft(j,d),kt(:,j,d),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)/ft(i,t)
if(ft(j,d) .GT. 0.0d0) then
arec = arec + vt(j,d)/ft(j,d)*kt(:,j,d)
end if
end if
end if
end do


!non-diffuse filtering begins

do i = jt+1, p
if(ymiss(d,i).EQ.0) then
vt(i,d) = yt(d,i) - ddot(m,zt(i,:,(d-1)*timevar(1)+1),1,arec,1) !vt
if (ft(i,d) .GT. tol*maxval(zt(i,:,(d-1)*timevar(1)+1)**2)) then !ft.NE.0
call daxpy(m,vt(i,d)/ft(i,d),kt(:,i,d),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)
vt(i,d) = yt(d,i) - ddot(m,zt(i,:,(d-1)*timevar(1)+1),1,arec,1)
if (ft(i,d) .GT. 0.0d0) then
arec = arec + vt(i,d)/ft(i,d)*kt(:,i,d)
end if
end if
end do
Expand All @@ -78,19 +74,16 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &
end if

if(dt.LT.n) then

!Non-diffuse filtering continues from t=d+1, i=1


if(dt.EQ.0) then
arec = a1
end if
do t = dt+1, n
do i = 1, p
if(ymiss(t,i).EQ.0) then
vt(i,t) = yt(t,i) - ddot(m,zt(i,:,(t-1)*timevar(1)+1),1,arec,1) !variate vt
if (ft(i,t) .GT. tol*maxval(zt(i,:,(t-1)*timevar(1)+1)**2)) then !ft.NE.0
call daxpy(m,vt(i,t)/ft(i,t),kt(:,i,t),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)
vt(i,t) = yt(t,i) - ddot(m,zt(i,:,(t-1)*timevar(1)+1),1,arec,1)
if (ft(i,t) .GT. 0.0d0) then
arec = arec + vt(i,t)/ft(i,t)*kt(:,i,t)
end if
end if
end do
Expand Down
42 changes: 23 additions & 19 deletions src/gloglik.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
double precision, dimension(m,p) :: kt,kinf
double precision, dimension(m,m) :: pt,pinf,mm
double precision, dimension(m,r) :: mr
double precision :: c, meps
double precision :: c, meps, finv
double precision, external :: ddot
double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr

external dgemm, dsymm, dgemv, dsymv, daxpy, dsyr, dsyr2, marginalxx
external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2, marginalxx

meps = epsilon(meps)

Expand All @@ -55,26 +55,28 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
pinf=p1inf
diffuse: do while(d .LT. n .AND. rankp .GT. 0)
d = d+1
do j=1, p !sequential processing
do j=1, p
if(ymiss(d,j).EQ.0) then
vt(j) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1)
call dsymv('u',m,1.0d0,pt,m,zt(j,:,(d-1)*timevar(1)+1),1,0.0d0,kt(:,j),1) ! kt_t,i = pt_t,i*t(z_t,i)
call dsymv('u',m,1.0d0,pt,m,zt(j,:,(d-1)*timevar(1)+1),1,0.0d0,kt(:,j),1)
ft(j) = ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,kt(:,j),1)+ ht(j,j,(d-1)*timevar(2)+1)
call dsymv('u',m,1.0d0,pinf,m,zt(j,:,(d-1)*timevar(1)+1),1,0.0d0,kinf(:,j),1) ! kinf_t,i = pinf_t,i*t(z_t,i)
call dsymv('u',m,1.0d0,pinf,m,zt(j,:,(d-1)*timevar(1)+1),1,0.0d0,kinf(:,j),1)
finf(j) = ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,kinf(:,j),1)
if (finf(j) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j)/finf(j),kinf(:,j),1,arec,1) !a_rec = a_rec + kinf(:,i,t)*vt(:,t)/finf(j,d)
call dsyr('u',m,ft(j)/finf(j)**2,kinf(:,j),1,pt,m) !pt = pt + kinf*kinf'*ft/finf^2
call dsyr2('u',m,-1.0d0/finf(j),kt(:,j),1,kinf(:,j),1,pt,m) !pt = pt -(kt*kinf'+kinf*kt')/finf
call dsyr('u',m,-1.0d0/finf(j),kinf(:,j),1,pinf,m) !pirec = pirec -kinf*kinf'/finf
finv = 1.0d0/finf(j)
arec = arec +vt(j)*finv*kinf(:,j)
call dsyr('u',m,ft(j)*finv**2,kinf(:,j),1,pt,m)
call dsyr2('u',m,-finv,kt(:,j),1,kinf(:,j),1,pt,m)
call dsyr('u',m,-finv,kinf(:,j),1,pinf,m)

lik = lik - 0.5d0*log(finf(j))
rankp = rankp -1
else
if (ft(j) .GT. tol*maxval(zt(j,:,(d-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(j)/ft(j),kt(:,j),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)/ft(i,t)
call dsyr('u',m,-1.0d0/ft(j),kt(:,j),1,pt,m) !pt = pt -kt*kt'/ft
lik = lik - c - 0.5d0*(log(ft(j)) + vt(j)**2/ft(j))
finv = 1.0d0/ft(j)
arec = arec + vt(j)*finv*kt(:,j)
call dsyr('u',m,-finv,kt(:,j),1,pt,m)
lik = lik - c - 0.5d0*(log(ft(j)) + vt(j)**2*finv)
end if
end if
if(rankp .EQ. 0) then
Expand Down Expand Up @@ -105,10 +107,11 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
vt(i) = yt(d,i) - ddot(m,zt(i,:,(d-1)*timevar(1)+1),1,arec,1)
call dsymv('u',m,1.0d0,pt,m,zt(i,:,(d-1)*timevar(1)+1),1,0.0d0,kt(:,i),1)
ft(i) = ddot(m,zt(i,:,(d-1)*timevar(1)+1),1,kt(:,i),1) + ht(i,i,(d-1)*timevar(2)+1)
if (ft(i).GT. tol*maxval(zt(i,:,(d-1)*timevar(1)+1)**2)) then !ft.NE.0
call daxpy(m,vt(i)/ft(i),kt(:,i),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)
call dsyr('u',m,-1.0d0/ft(i),kt(:,i),1,pt,m) !p_rec = p_rec - kt*kt'*ft(i,t)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))-c
if (ft(i).GT. tol*maxval(zt(i,:,(d-1)*timevar(1)+1)**2)) then
finv = 1.0d0/ft(i)
arec = arec + vt(i)*finv*kt(:,i)
call dsyr('u',m,-finv,kt(:,i),1,pt,m)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2*finv)-c
end if
end if
end do
Expand All @@ -133,9 +136,10 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
call dsymv('u',m,1.0d0,pt,m,zt(i,:,(t-1)*timevar(1)+1),1,0.0d0,kt(:,i),1)
ft(i) = ddot(m,zt(i,:,(t-1)*timevar(1)+1),1,kt(:,i),1) + ht(i,i,(t-1)*timevar(2)+1)
if (ft(i).GT. tol*maxval(zt(i,:,(t-1)*timevar(1)+1)**2)) then
call daxpy(m,vt(i)/ft(i),kt(:,i),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)
call dsyr('u',m,-1.0d0/ft(i),kt(:,i),1,pt,m) !pt = pt - kt*kt'*ft(i,i,t)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))-c
finv = 1.0d0/ft(i)
arec = arec + vt(i)*finv*kt(:,i)
call dsyr('u',m,-finv,kt(:,i),1,pt,m)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2*finv)-c
end if
end if
end do
Expand Down
Loading

0 comments on commit 585b4b0

Please sign in to comment.