1!diffuse filtering for single time point
2subroutine dfilter1stepnv(ymiss, yt, zt, tt, at, vt, ft,kt,&
3finf,kinf,p,m,j,lik)
4
5    implicit none
6
7    integer, intent(in) ::  p, m, j
8    integer :: i
9    integer, intent(in), dimension(p) :: ymiss
10    double precision, intent(in), dimension(p) :: yt
11    double precision, intent(in), dimension(m,p) :: zt
12    double precision, intent(in), dimension(m,m) :: tt
13    double precision, intent(inout), dimension(m) :: at
14    double precision, intent(inout), dimension(p) :: vt
15    double precision, intent(in), dimension(p) :: ft,finf
16    double precision, intent(in), dimension(m,p) :: kt,kinf
17    double precision, intent(inout) :: lik
18    double precision, dimension(m) :: ahelp
19
20    double precision, external :: ddot
21    external dgemv
22
23    do i = 1, j
24        if(ymiss(i).EQ.0) then
25            vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
26            if (finf(i) .GT. 0.0d0) then
27                at = at + vt(i)/finf(i)*kinf(:,i)
28                lik = lik - 0.5d0*log(finf(i))
29            else
30                if (ft(i) .GT. 0.0d0) then
31                    at = at + vt(i)/ft(i)*kt(:,i)
32                    lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))
33                end if
34            end if
35        end if
36    end do
37    if(j .EQ. p) then
38        call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
39        at = ahelp
40    end if
41end subroutine dfilter1stepnv
42
43
44!non-diffuse filtering for single time point
45subroutine filter1stepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j,lik)
46
47    implicit none
48
49    integer, intent(in) ::  p, m,j
50    integer ::  i
51    integer, intent(in), dimension(p) :: ymiss
52    double precision, intent(in), dimension(p) :: yt
53    double precision, intent(in), dimension(m,p) :: zt
54    double precision, intent(in), dimension(m,m) :: tt
55    double precision, intent(inout), dimension(m) :: at
56    double precision, intent(in), dimension(p) :: ft
57    double precision, intent(inout), dimension(p) :: vt
58    double precision, intent(in), dimension(m,p) :: kt
59    double precision, intent(inout) :: lik
60    double precision, dimension(m) :: ahelp
61
62    double precision, external :: ddot
63
64    external dgemv
65
66    do i = j+1, p
67        if(ymiss(i).EQ.0) then
68            vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1)
69            if (ft(i) .GT. 0.0d0) then
70                at = at + vt(i)/ft(i)*kt(:,i)
71                lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i))
72            end if
73        end if
74    end do
75
76    call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1)
77    at = ahelp
78
79end subroutine filter1stepnv
80