Skip to content

Commit

Permalink
further modularization
Browse files Browse the repository at this point in the history
  • Loading branch information
@helske
helske committed Apr 16, 2015
1 parent e78ea8d commit ddacec9
Show file tree
Hide file tree
Showing 3 changed files with 118 additions and 141 deletions.
75 changes: 75 additions & 0 deletions src/filteronestepnovar.f95
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
!non-diffuse filtering for single time point
subroutine filteronestepnv(ymiss, yt, zt, tt, at,vt, ft,kt,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(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, dimension(m) :: ahelp

double precision, external :: ddot

external dgemv

do i = j+1, p
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)
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


!diffuse filtering for single time point
subroutine diffusefilteronestepnv(ymiss, yt, zt, tt, at, vt, ft,kt,&
finf,kinf,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(m,m) :: tt
double precision, intent(inout), dimension(m) :: at
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, dimension(m) :: ahelp

double precision, external :: ddot

external dgemv

do i = 1, j
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
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
89 changes: 22 additions & 67 deletions src/filtersimfast.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &
implicit none

integer, intent(in) :: p, m,n,dt,jt
integer :: t, i,d,j
integer :: t, j
integer, intent(in), dimension(n,p) :: ymiss
integer, intent(in), dimension(5) :: timevar
double precision, intent(in), dimension(n,p) :: yt
Expand All @@ -22,76 +22,31 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, &

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)
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. 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)
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)

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. 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)
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

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)
end if
at(:,1) = a1
!diffuse filtering begins
do t = 1, (dt - 1)
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,p)
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
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)
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

call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,arec,1,0.0d0,at(:,t+1),1)
arec = at(:,t+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
95 changes: 21 additions & 74 deletions src/smoothsimfast.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,&

logical, intent(in) :: needeps
integer, intent(in) :: p, m, r,n,dt,jt
integer :: t, i,d,j
integer :: t, i,j
integer, intent(in), dimension(n,p) :: ymiss
integer, intent(in), dimension(5) :: timevar
double precision, intent(in), dimension(n,p) :: yt
Expand All @@ -21,88 +21,35 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,&
double precision, intent(inout), dimension(p,n) :: epshat
double precision, intent(inout), dimension(r,n) :: etahat
double precision, dimension(p,n) :: vt
double precision, dimension(m) :: arec,rrec,rrec1,rhelp,help
double precision, dimension(m,m) :: im,linf,l0
double precision, dimension(m) :: at
double precision, dimension(m,m) :: im
double precision, intent(inout), dimension(m) :: rt0,rt1
double precision :: finv
double precision, external :: ddot

external dgemv, dger, dsymv

j = 0
d = 0
if(dt.GT.0) then
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)
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. 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,help,1)
arec = help
end do diffuse
at = a1
!diffuse filtering begins
do t = 1, dt - 1
call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),&
finf(:,t),kinf(:,:,t),p,m,p)
end do

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)

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. 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)
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

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

end if
t = dt
call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at,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,vt(:,t),ft(:,t),kt(:,:,t),p,m,jt)

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)
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
call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,arec,1,0.0d0,help,1)
arec = help
do t = dt + 1, n
call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,0)
end do

end if


Expand All @@ -115,7 +62,7 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,&

rt0 = 0.0d0

do t = n, d+1, -1
do t = n, dt+1, -1
call smoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat(:,t),needeps)
Expand Down

0 comments on commit ddacec9

Please sign in to comment.