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