1! compute the additional term for marginal likelihood 2subroutine marginalxx(p1inf,zt,tt,m,p,n,k,timevar,lik,info) 3 4 integer, intent(inout) :: info 5 integer, intent(in) :: m,p,n,k 6 integer, intent(in), dimension(5) :: timevar 7 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 8 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 9 double precision, intent(in), dimension(m,m) :: p1inf 10 double precision, intent(inout) :: lik 11 integer ::j,i 12 double precision, dimension(m,k) :: a,a2 13 double precision, dimension(k,k) :: s 14 double precision, dimension(p,k) :: v 15 16 external dgemm, dpotrf, dsyrk 17 18 a=0.0d0 19 j=1 20 do i=1, m 21 if(sum(p1inf(:,i)).GT.0.0d0) then 22 a(i,j) = 1.0d0 23 j = j+1 24 end if 25 end do 26 s=0.0d0 27 do i = 1, n 28 call dgemm('n','n',p,k,m,1.0d0,zt(:,:,(i-1)*timevar(1)+1),p,a,m,0.0d0,v,p) 29 call dgemm('n','n',m,k,m,1.0d0,tt(:,:,(i-1)*timevar(3)+1),m,a,m,0.0d0,a2,m) 30 a = a2 31 call dsyrk('u','t',k,p,1.0d0,v,p,1.0d0,s,k) 32 end do 33 call dpotrf('u', k, s, k, info) 34 if(info.EQ.0) then 35 do i=1, k 36 lik = lik + log(s(i,i)) 37 end do 38 else 39 info = -1 40 end if 41 42end subroutine marginalxx 43