From ddacec94cc994aafa2b7fda97493ff84e60640cb Mon Sep 17 00:00:00 2001 From: "jouni.helske@jyu.fi" Date: Thu, 16 Apr 2015 23:03:14 +0300 Subject: [PATCH] further modularization --- src/filteronestepnovar.f95 | 75 ++++++++++++++++++++++++++++++ src/filtersimfast.f95 | 89 +++++++++-------------------------- src/smoothsimfast.f95 | 95 +++++++++----------------------------- 3 files changed, 118 insertions(+), 141 deletions(-) create mode 100644 src/filteronestepnovar.f95 diff --git a/src/filteronestepnovar.f95 b/src/filteronestepnovar.f95 new file mode 100644 index 0000000..7bf86a4 --- /dev/null +++ b/src/filteronestepnovar.f95 @@ -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 diff --git a/src/filtersimfast.f95 b/src/filtersimfast.f95 index cf4aa1c..9f27766 100644 --- a/src/filtersimfast.f95 +++ b/src/filtersimfast.f95 @@ -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 @@ -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 diff --git a/src/smoothsimfast.f95 b/src/smoothsimfast.f95 index 4a54e42..c87bc67 100644 --- a/src/smoothsimfast.f95 +++ b/src/smoothsimfast.f95 @@ -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 @@ -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 @@ -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)