1! functions for computing p(theta)
2
3! Only differences with gloglik is what is stored/returned (ft,finf,kt,kinf)
4! Also no missing values as yt is actually the signals, and no ht
5subroutine pthetafirst(yt, timevar, zt, tt, rqr, a1, p1, p1inf,&
6p, m, n, lik, tol,rankp2,kt,kinf,ft,finf,d,j)
7
8
9    implicit none
10
11    integer, intent(in) ::  p, m, n
12    integer, intent(inout) :: rankp2,d,j
13    integer ::  t, tv,rankp
14    integer, intent(in), dimension(5) :: timevar
15    double precision, intent(in), dimension(n,p) :: yt
16    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
17    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
18    double precision, intent(in), dimension(m) :: a1
19    double precision, intent(in), dimension(m,m) ::  p1,p1inf
20    double precision, intent(in) :: tol
21    double precision, intent(inout) :: lik
22    double precision, dimension(m) :: at
23    double precision, dimension(p) :: vt
24    double precision, intent(inout), dimension(p,n) :: ft,finf
25    double precision, intent(inout), dimension(m,p,n) :: kt,kinf
26    double precision, dimension(m,m) :: pt,pinf
27    double precision, intent(inout), dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr
28    integer, dimension(p) :: ymiss
29    double precision, dimension(p,p) :: ht
30
31    external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2
32
33    tv= max(timevar(4),timevar(5))
34    ymiss = 0
35    ht = 0.0d0
36
37    rankp = rankp2
38    j=0
39    d=0
40    at = a1
41    pt = p1
42    pinf=p1inf
43
44    ! Diffuse initialization
45    if(rankp .GT. 0) then
46        diffuse: do while(d .LT. n .AND. rankp .GT. 0)
47            d = d+1
48            call dfilter1step(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,&
49            tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),&
50            at,pt,vt,ft(:,d),kt(:,:,d),pinf,finf(:,d),kinf(:,:,d),rankp,lik,tol,0.0d0,p,m,j)
51        end do diffuse
52        if(rankp .EQ. 0 .AND. j .LT. p) then
53            call filter1step(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,&
54            tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),&
55            at,pt,vt,ft(:,d),kt(:,:,d),lik,tol,0.0d0,p,m,j)
56        else
57            j = p
58        end if
59    end if
60
61    !Non-diffuse filtering continues from t=d+1, i=1
62
63    do t = d+1, n
64        call filter1step(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht,&
65        tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),&
66        at,pt,vt,ft(:,t),kt(:,:,t),lik,tol,0.0d0,p,m,0)
67    end do
68
69
70
71
72end subroutine pthetafirst
73
74
75! use output of pthetafirst, kt and ft do not change
76subroutine pthetarest(yt, timevar, zt, tt, a1,&
77p, m, n, lik, kt,kinf,ft,finf,dt,jt)
78
79
80    implicit none
81
82    integer, intent(in) ::  p, m, n,dt,jt
83    integer ::  t
84    integer, intent(in), dimension(5) :: timevar
85    double precision, intent(in), dimension(n,p) :: yt
86    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
87    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
88    double precision, intent(in), dimension(m) :: a1
89    double precision, intent(inout) :: lik
90    double precision, dimension(m) :: at
91    double precision, dimension(p) :: vt
92    double precision, intent(in), dimension(p,n) :: ft,finf
93    double precision, intent(in), dimension(m,p,n) :: kt,kinf
94    integer, dimension(p) :: ymiss
95    double precision, external :: ddot
96
97    external dgemv
98    ymiss = 0
99    at = a1
100    if(dt .GT. 0) then
101        !diffuse filtering begins
102        do t = 1, dt - 1
103            call dfilter1stepnv(ymiss,yt(t,:),&
104            transpose(zt(:,:,(t-1)*timevar(1)+1)),tt(:,:,(t-1)*timevar(3)+1),&
105            at,vt,ft(:,t),kt(:,:,t), finf(:,t),kinf(:,:,t),p,m,p,lik)
106        end do
107
108        t = dt
109        call dfilter1stepnv(ymiss,yt(t,:),&
110        transpose(zt(:,:,(t-1)*timevar(1)+1)),tt(:,:,(t-1)*timevar(3)+1),&
111        at,vt,ft(:,t),kt(:,:,t),finf(:,t),kinf(:,:,t),p,m,jt,lik)
112        !non-diffuse filtering begins
113        if(jt .LT. p) then
114            call filter1stepnv(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
115            tt(:,:,(t-1)*timevar(3)+1),at,vt,ft(:,t),kt(:,:,t),p,m,jt,lik)
116        end if
117    end if
118    !Non-diffuse filtering continues from t=d+1, i=1
119    do t = dt + 1, n
120        call filter1stepnv(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
121        tt(:,:,(t-1)*timevar(3)+1),at,vt,ft(:,t),kt(:,:,t),p,m,0,lik)
122    end do
123
124end subroutine pthetarest
125
126