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