Skip to content

Commit

Permalink
Modularization complete
Browse files Browse the repository at this point in the history
  • Loading branch information
@helske
helske committed Apr 17, 2015
1 parent ddacec9 commit 35fb413
Show file tree
Hide file tree
Showing 20 changed files with 2,054 additions and 439 deletions.
109 changes: 55 additions & 54 deletions src/filteronestep.f95 → src/filter1step.f95
View file
Original file line number Diff line number Diff line change
@@ -1,58 +1,6 @@
!non-diffuse filtering for single time point
subroutine filteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,lik,tol,meps,c,p,m,j)

implicit none

integer, intent(in) :: p, m,j
integer :: i
integer, intent(in), dimension(p) :: ymiss
double precision, intent(in), dimension(p) :: yt
double precision, intent(in), dimension(m,p) :: zt
double precision, intent(in), dimension(p,p) :: ht
double precision, intent(in), dimension(m,m) :: tt
double precision, dimension(m,m) :: rqr
double precision, intent(in) :: tol,c,meps
double precision, intent(inout) :: lik
double precision, intent(inout), dimension(m) :: at
double precision, intent(inout), dimension(p) :: vt,ft
double precision, intent(inout), dimension(m,p) :: kt
double precision, intent(inout), dimension(m,m) :: pt
double precision, dimension(m,m) :: mm
double precision, dimension(m) :: ahelp
double precision :: finv

double precision, external :: ddot

external dgemm, dsymm, dgemv, dsymv, dsyr

do i = j+1, p
call dsymv('u',m,1.0d0,pt,m,zt(:,i),1,0.0d0,kt(:,i),1)
ft(i) = ddot(m,zt(:,i),1,kt(:,i),1) + ht(i,i)
if(ymiss(i).EQ.0) then
vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
if (ft(i) .GT. tol*maxval(zt(:,i)**2)) then
finv = 1.0d0/ft(i)
at = at + vt(i)*finv*kt(:,i)
call dsyr('u',m,-finv,kt(:,i),1,pt,m)
lik = lik - c - 0.5d0*(log(ft(i)) + vt(i)**2*finv)
else
ft(i)=0.0d0
end if
end if
end do

call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp

call dsymm('r','u',m,m,1.0d0,pt,m,tt,m,0.0d0,mm,m)
call dgemm('n','t',m,m,m,1.0d0,mm,m,tt,m,0.0d0,pt,m)
pt = pt + rqr

end subroutine filteronestep


!diffuse filtering for single time point
subroutine diffusefilteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,&
subroutine dfilter1step(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,&
pinf,finf,kinf,rankp,lik,tol,meps,c,p,m,i)

implicit none
Expand Down Expand Up @@ -127,4 +75,57 @@ subroutine diffusefilteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,&
end if
end do

end subroutine diffusefilteronestep
end subroutine dfilter1step

!non-diffuse filtering for single time point
subroutine filter1step(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,lik,tol,c,p,m,j)

implicit none

integer, intent(in) :: p, m,j
integer :: i
integer, intent(in), dimension(p) :: ymiss
double precision, intent(in), dimension(p) :: yt
double precision, intent(in), dimension(m,p) :: zt
double precision, intent(in), dimension(p,p) :: ht
double precision, intent(in), dimension(m,m) :: tt
double precision, dimension(m,m) :: rqr
double precision, intent(in) :: tol,c
double precision, intent(inout) :: lik
double precision, intent(inout), dimension(m) :: at
double precision, intent(inout), dimension(p) :: vt,ft
double precision, intent(inout), dimension(m,p) :: kt
double precision, intent(inout), dimension(m,m) :: pt
double precision, dimension(m,m) :: mm
double precision, dimension(m) :: ahelp
double precision :: finv

double precision, external :: ddot

external dgemm, dsymm, dgemv, dsymv, dsyr

do i = j+1, p
call dsymv('u',m,1.0d0,pt,m,zt(:,i),1,0.0d0,kt(:,i),1)
ft(i) = ddot(m,zt(:,i),1,kt(:,i),1) + ht(i,i)
if(ymiss(i).EQ.0) then
vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
if (ft(i) .GT. tol*maxval(zt(:,i)**2)) then
finv = 1.0d0/ft(i)
at = at + vt(i)*finv*kt(:,i)
call dsyr('u',m,-finv,kt(:,i),1,pt,m)
lik = lik - c - 0.5d0*(log(ft(i)) + vt(i)**2*finv)
else
ft(i)=0.0d0
end if
end if
end do

call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp

call dsymm('r','u',m,m,1.0d0,pt,m,tt,m,0.0d0,mm,m)
call dgemm('n','t',m,m,m,1.0d0,mm,m,tt,m,0.0d0,pt,m)
pt = pt + rqr

end subroutine filter1step

74 changes: 39 additions & 35 deletions src/filteronestepnovar.f95 → src/filter1stepnovar.f95
View file
Original file line number Diff line number Diff line change
@@ -1,75 +1,79 @@
!non-diffuse filtering for single time point
subroutine filteronestepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j)
!diffuse filtering for single time point
subroutine dfilter1stepnv(ymiss, yt, zt, tt, at, vt, ft,kt,&
finf,kinf,p,m,j,lik)

implicit none

integer, intent(in) :: p, m,j
integer :: i
integer, intent(in) :: p, m, j
integer :: i
integer, intent(in), dimension(p) :: ymiss
double precision, intent(in), dimension(p) :: yt
double precision, intent(in), dimension(m,p) :: zt
double precision, intent(in), dimension(m,m) :: tt
double precision, intent(inout), dimension(m) :: at
double precision, intent(in), dimension(p) :: ft
double precision, intent(inout), dimension(p) :: vt
double precision, intent(in), dimension(m,p) :: kt
double precision, intent(in), dimension(p) :: ft,finf
double precision, intent(in), dimension(m,p) :: kt,kinf
double precision, intent(inout) :: lik
double precision, dimension(m) :: ahelp

double precision, external :: ddot

double precision, external :: ddot
external dgemv

do i = j+1, p
do i = 1, j
if(ymiss(i).EQ.0) then
vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
if (ft(i) .GT. 0.0d0) then
at = at + vt(i)/ft(i)*kt(:,i)
if (finf(i) .GT. 0.0d0) then
at = at + vt(i)/finf(i)*kinf(:,i)
lik = lik - 0.5d0*log(finf(i))
else
if (ft(i) .GT. 0.0d0) then
at = at + vt(i)/ft(i)*kt(:,i)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))
end if
end if
end if
end do

call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp

end subroutine filteronestepnv
if(j .EQ. p) then
call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp
end if
end subroutine dfilter1stepnv


!diffuse filtering for single time point
subroutine diffusefilteronestepnv(ymiss, yt, zt, tt, at, vt, ft,kt,&
finf,kinf,p,m,j)
!non-diffuse filtering for single time point
subroutine filter1stepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j,lik)

implicit none

integer, intent(in) :: p, m, j
integer :: i
integer, intent(in) :: p, m,j
integer :: i
integer, intent(in), dimension(p) :: ymiss
double precision, intent(in), dimension(p) :: yt
double precision, intent(in), dimension(m,p) :: zt
double precision, intent(in), dimension(m,m) :: tt
double precision, intent(inout), dimension(m) :: at
double precision, intent(in), dimension(p) :: ft
double precision, intent(inout), dimension(p) :: vt
double precision, intent(in), dimension(p) :: ft,finf
double precision, intent(in), dimension(m,p) :: kt,kinf
double precision, intent(in), dimension(m,p) :: kt
double precision, intent(inout) :: lik
double precision, dimension(m) :: ahelp

double precision, external :: ddot

external dgemv

do i = 1, j
do i = j+1, p
if(ymiss(i).EQ.0) then
vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
if (finf(i) .GT. 0.0d0) then
at = at + vt(i)/finf(i)*kinf(:,i)
else
if (ft(i) .GT. 0.0d0) then
at = at + vt(i)/ft(i)*kt(:,i)
end if
if (ft(i) .GT. 0.0d0) then
at = at + vt(i)/ft(i)*kt(:,i)
lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))
end if
end if
end do
if(j .EQ. p) then
call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp
end if
end subroutine diffusefilteronestepnv

call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
at = ahelp

end subroutine filter1stepnv
53 changes: 28 additions & 25 deletions src/filtersimfast.f95
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
!fast filtering algorithm used in simulation filter
subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &
a1, ft,kt,finf, kinf, dt, jt, p, m, n,tol,at)
a1, ft,kt,finf, kinf, dt, jt, p, m, n,at)

implicit none

integer, intent(in) :: p, m,n,dt,jt
integer :: t, j
integer :: t
integer, intent(in), dimension(n,p) :: ymiss
integer, intent(in), dimension(5) :: timevar
double precision, intent(in), dimension(n,p) :: yt
Expand All @@ -14,39 +14,42 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &
double precision, intent(in), dimension(m) :: a1
double precision, intent(in), dimension(p,n) :: ft,finf
double precision, intent(in), dimension(m,p,n) :: kt,kinf
double precision, intent(in) :: tol
double precision, intent(inout), dimension(m,n+1) :: at
double precision, dimension(p,n) :: vt
double precision, dimension(m) :: arec
double precision :: lik
double precision, external :: ddot

external dgemv

lik = 0.0d0

at(:,1) = a1
!diffuse filtering begins
do t = 1, (dt - 1)
if(dt .GT. 0) then
!diffuse filtering begins
do t = 1, (dt - 1)
at(:,t+1) = at(:,t)
call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),&
finf(:,t),kinf(:,:,t),p,m,p,lik)
end do

t = dt
at(:,t+1) = at(:,t)
call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),&
finf(:,t),kinf(:,:,t),p,m,p)
finf(:,t),kinf(:,:,t),p,m,jt,lik)
!non-diffuse filtering begins
if(jt .LT. p) then
call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,jt,lik)
end if
end if
!Non-diffuse filtering continues from t=d+1, i=1
do t = dt + 1, n
at(:,t+1) = at(:,t)
call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,0,lik)
end do

t = dt
at(:,t+1) = at(:,t)
call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),&
finf(:,t),kinf(:,:,t),p,m,jt)
!non-diffuse filtering begins
call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,jt)

if(dt.LT.n) then
!Non-diffuse filtering continues from t=d+1, i=1
do t = dt + 1, n
at(:,t+1) = at(:,t)
call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,0)
end do
end if

end subroutine filtersimfast
18 changes: 9 additions & 9 deletions src/gloglik.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&

integer, intent(in) :: p, m, r, n
integer, intent(inout) :: rankp,marginal
integer :: t, i,d,j,tv
integer :: t,d,j,tv
integer, intent(in), dimension(n,p) :: ymiss
integer, intent(in), dimension(5) :: timevar
double precision, intent(in), dimension(n,p) :: yt
Expand All @@ -21,12 +21,12 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
double precision, intent(in), dimension(m,m) :: p1,p1inf
double precision, intent(in) :: tol
double precision, intent(inout) :: lik
double precision, dimension(m) :: at,arec
double precision, dimension(m) :: at
double precision, dimension(p) :: vt,ft,finf
double precision, dimension(m,p) :: kt,kinf
double precision, dimension(m,m) :: pt,pinf,mm
double precision, dimension(m,m) :: pt,pinf
double precision, dimension(m,r) :: mr
double precision :: c, meps, finv
double precision :: c, meps
double precision, external :: ddot
double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr

Expand Down Expand Up @@ -56,15 +56,15 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
diffuse: do while(d .LT. n .AND. rankp .GT. 0)
d = d+1

call diffusefilteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1), at,pt,vt,ft,kt,pinf,finf,kinf,rankp,lik,tol,meps,c,p,m,j)

end do diffuse

if(rankp .EQ. 0 .AND. j .LT. p) then
!non-diffuse filtering begins
call filteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),at,pt,vt,ft,kt,lik,tol,meps,c,p,m,j)
call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),at,pt,vt,ft,kt,lik,tol,c,p,m,j)

else
j = p
Expand All @@ -77,8 +77,8 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,&
!Non-diffuse filtering continues from t=d+1, i=1

do t = d+1, n
call filteronestep(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),&
tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),at,pt,vt,ft,kt,lik,tol,meps,c,p,m,0)
call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),&
tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),at,pt,vt,ft,kt,lik,tol,c,p,m,0)

end do

Expand Down
Loading

0 comments on commit 35fb413

Please sign in to comment.