1! Transform multivariate model using LDL decomposition of H
2subroutine ldlssm(yt, ydimt, yobs, timevar, zt, p, m, n, ichols,nh,hchol,dim,info,hobs,tol)
3
4    implicit none
5
6    integer, intent(in) ::  p, m, n,nh
7    integer, intent(in), dimension(nh) ::  dim
8    integer, intent(in), dimension(p,nh) :: hobs
9    integer, intent(in), dimension(n) :: ydimt
10    integer, intent(in), dimension(p,n) :: yobs
11    integer, intent(in), dimension(5) :: timevar
12    integer, intent(inout), dimension(n) :: hchol
13    integer, intent(inout) :: info
14    integer, dimension(nh) :: hdiagtest
15    integer ::  t, i, k
16    double precision, intent(inout), dimension(p,p,nh) :: ichols
17    double precision, intent(inout), dimension(p,n) :: yt
18    double precision, intent(inout), dimension(p,m,(n-1)*timevar(1)+1) :: zt
19    double precision, intent(in) :: tol
20    double precision, dimension(p,p) :: cholhelp
21    double precision, dimension(p) :: yhelp
22    double precision, dimension(p,m) :: zhelp
23    double precision, external :: ddot
24
25    external dtrmv, dtrmm, dtrtri, ldl
26
27    hdiagtest = 0
28    info=0
29
30    if(p.GT.1) then !testing diagonality
31        do t= 1, nh
32            test: do i = 1, dim(t)
33                do k = i+1, dim(t)
34                    if(abs(ichols(hobs(k,t),hobs(i,t),t)) .GT. tol) then
35                        hdiagtest(t)=1
36                        exit test
37                    end if
38                end do
39            end do test
40        end do
41        if(sum(hdiagtest).NE.0) then
42            do i = 1, nh
43                if(hdiagtest(i).EQ.1) then
44                    cholhelp(1:dim(i),1:dim(i)) = ichols(hobs(1:dim(i),i),hobs(1:dim(i),i),i)
45                    call ldl(cholhelp(1:dim(i),1:dim(i)),dim(i),tol,info)
46                    if(info .EQ. -1) then
47                        info=1
48                        return
49                    end if
50                    call dtrtri('l','u',dim(i),cholhelp(1:dim(i),1:dim(i)),dim(i),info)
51                    if(info .NE. 0) then
52                        info=2
53                        return
54                    end if
55                    ichols(hobs(1:dim(i),i),hobs(1:dim(i),i),i) = cholhelp(1:dim(i),1:dim(i))
56               end if
57
58            end do
59
60            do t = 1,(n-1)*max(timevar(1),timevar(2))+1
61                if(ydimt(t).GT.0 .AND. hdiagtest(hchol(t)).NE.0) then
62                    cholhelp(1:ydimt(t),1:ydimt(t)) = ichols(yobs(1:ydimt(t),t),yobs(1:ydimt(t),t),hchol(t))
63                    zhelp(1:ydimt(t),:) = zt(yobs(1:ydimt(t),t),:,(t-1)*timevar(1)+1)
64                    call dtrmm('l','l','n','u',ydimt(t),m,1.0d0,cholhelp(1:ydimt(t),1:ydimt(t)),&
65                    ydimt(t),zhelp(1:ydimt(t),:),ydimt(t))
66                    zt(yobs(1:ydimt(t),t),:,(t-1)*timevar(1)+1) = zhelp(1:ydimt(t),:)
67
68                end if
69            end do
70            do t= 1, n
71                if(ydimt(t).GT.0  .AND. hdiagtest(hchol(t)).NE.0) then
72                    cholhelp(1:ydimt(t),1:ydimt(t)) = ichols(yobs(1:ydimt(t),t),yobs(1:ydimt(t),t),hchol(t))
73                    yhelp(1:ydimt(t)) = yt(yobs(1:ydimt(t),t),t)
74                    call dtrmv('l','n','u',ydimt(t),cholhelp(1:ydimt(t),1:ydimt(t)),ydimt(t),yhelp(1:ydimt(t)),1)
75                    yt(yobs(1:ydimt(t),t),t) = yhelp(1:ydimt(t))
76                end if
77            end do
78        end if
79    end if
80
81
82end subroutine ldlssm
83