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