1!simulate gaussian state space model
2subroutine simgaussianuncond(timevar, zt, ht, tt, rtv, qt, a1, p1, &
3nnd,nsim, epsplus, etaplus, aplus1, p, n, m, r, info,&
4tol,sim,c,simwhat,simdim,antithetics)
5
6    implicit none
7
8    integer, intent(in) :: p, m, r, n, nsim,nnd,simdim,simwhat,antithetics
9    integer, intent(in), dimension(5) :: timevar
10    integer, intent(inout) :: info
11    integer ::  t, i,k,l
12    double precision, intent(in) :: tol
13    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
14    double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht
15    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
16    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
17    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
18    double precision, intent(in), dimension(m) :: a1
19    double precision, intent(in), dimension(m,m) ::  p1
20    double precision, intent(in),dimension(nsim) :: c
21    double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim
22    double precision, intent(inout), dimension(r,n,nsim) :: etaplus
23    double precision, intent(inout), dimension(p,n,nsim) :: epsplus
24    double precision, intent(inout), dimension(m,nsim) :: aplus1
25
26    double precision, dimension(m,n+1) :: aplus
27    double precision, dimension(r,r,(n-1)*timevar(5)+1) :: cholqt
28    double precision, dimension(m,m) :: cholp1
29    double precision, dimension(r,r) :: rcholtmp
30    double precision, dimension(r,n-1,4) :: etatmp
31    double precision, dimension(m,n,4) :: alphatmp
32    double precision, external :: ddot
33
34    external dsymm, dgemm, dsymv, ldl, dtrmv, dgemv
35
36
37
38    do t = 1, (n-1)*timevar(5)+1
39        if(r.EQ.1) then
40            cholqt(1,1,t)=sqrt(qt(1,1,t))
41        else
42            rcholtmp = qt(:,:,t)
43            call ldl(rcholtmp,r,tol,info)
44            if(info .NE. 0) then
45                info = -2
46                return
47            end if
48            do i=1,r
49                cholqt(i,i,t)=sqrt(rcholtmp(i,i))
50            end do
51            do i=1,r-1
52                cholqt((i+1):r,i,t) = rcholtmp((i+1):r,i)*cholqt(i,i,t)
53            end do
54        end if
55    end do
56
57
58    if(nnd.GT.0) then
59        if(m.EQ.1) then
60            cholp1(1,1)=sqrt(p1(1,1))
61        else
62            cholp1 = p1
63            call ldl(cholp1,m,tol,info)
64            if(info .NE. 0) then
65                info=-3
66                return
67            end if
68            do i=1,m
69                cholp1(i,i)=sqrt(cholp1(i,i))
70            end do
71            do i=1,m-1
72                cholp1((i+1):m,i) = cholp1((i+1):m,i)*cholp1(i,i)
73            end do
74        end if
75    end if
76
77
78    do i = 1, nsim
79        aplus = 0.0d0
80        aplus(:,1) = a1
81        if(nnd.GT.0) then
82            call dtrmv('l','n','n',m,cholp1,m,aplus1(:,i),1)
83            aplus(:,1) = aplus(:,1)+aplus1(:,i)
84        end if
85
86        do t = 1, n
87            do k = 1, p
88                epsplus(k,t,i) = epsplus(k,t,i)*sqrt(ht(k,k,(t-1)*timevar(2)+1))
89            end do
90            call dtrmv('l','n','n',r,cholqt(:,:,(t-1)*timevar(5)+1),r,etaplus(:,t,i),1)
91            call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,aplus(:,t),1,0.0d0,aplus(:,t+1),1)
92            call dgemv('n',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,etaplus(:,t,i),1,1.0d0,aplus(:,t+1),1)
93        end do
94
95
96        !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations
97        select case(simwhat)
98            case(1)
99                sim(:,:,i) = epsplus(:,:,i)
100                if(antithetics .EQ. 1) then
101                    sim(:,:,i+nsim) =  -epsplus(:,:,i)
102                    sim(:,:,i+2*nsim) = c(i)*sim(:,:,i)
103                    sim(:,:,i+3*nsim) = c(i)*sim(:,:,i+nsim)
104                end if
105            case(2)
106                sim(:,:,i) = etaplus(:,:,i)
107                if(antithetics .EQ. 1) then
108                    sim(:,:,i+nsim) = -etaplus(:,:,i)
109                    sim(:,:,i+2*nsim) = c(i)*sim(:,:,i)
110                    sim(:,:,i+3*nsim) = c(i)*sim(:,:,i+nsim)
111                end if
112            case(3)
113                sim(1:p,:,i) = epsplus(:,:,i)
114                sim((p+1):,:,i) = etaplus(:,:,i)
115                if(antithetics .EQ. 1) then
116                    sim(1:p,:,i+nsim) = -epsplus(:,:,i)
117                    sim(1:p,:,i+2*nsim) = c(i)*sim(1:p,:,i)
118                    sim(1:p,:,i+3*nsim) = c(i)*sim(1:p,:,i+nsim)
119                    sim((p+1):,:,i+nsim) = -etaplus(:,:,i)
120                    sim((p+1):,:,i+2*nsim) = c(i)*sim((p+1):,:,i)
121                    sim((p+1):,:,i+3*nsim) = c(i)*sim((p+1):,:,i+nsim)
122                end if
123            case(4)
124                sim(:,1,i) = aplus(:,1)
125                etatmp(:,:,1) = etaplus(:,1:(n-1),i)
126                if(antithetics .EQ. 1) then
127                    sim(:,1,i+nsim) = -aplus(:,1)
128                    sim(:,1,i+2*nsim) = c(i)*sim(:,1,i)
129                    sim(:,1,i+3*nsim) = c(i)*sim(:,1,i+nsim)
130
131                    etatmp(:,:,2) = -etaplus(:,1:(n-1),i)
132                    etatmp(:,:,3) = c(i)*etatmp(:,:,1)
133                    etatmp(:,:,4) = c(i)*etatmp(:,:,2)
134                end if
135                do k = 1, 3*antithetics+1
136                    do t = 2, n
137                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,sim(:,t-1,i+(k-1)*nsim),&
138                        1,0.0d0,sim(:,t,i+(k-1)*nsim),1)
139                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
140                        1.0d0,sim(:,t,i+(k-1)*nsim),1)
141                    end do
142                end do
143            case(5)
144
145                alphatmp(:,1,1) = aplus(:,1)
146                etatmp(:,:,1) =  etaplus(:,1:(n-1),i)
147                call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,1),1,0.0d0,sim(:,1,i),1)
148
149                if(antithetics .EQ. 1) then
150                    alphatmp(:,1,2) = -aplus(:,1)
151                    alphatmp(:,1,3) = c(i)*alphatmp(:,1,1)
152                    alphatmp(:,1,4) = c(i)*alphatmp(:,1,2)
153                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,2),1,0.0d0,sim(:,1,i+nsim),1)
154                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,3),1,0.0d0,sim(:,1,i+2*nsim),1)
155                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,4),1,0.0d0,sim(:,1,i+3*nsim),1)
156
157                    etatmp(:,:,2) = -etaplus(:,1:(n-1),i)
158                    etatmp(:,:,3) = c(i)*etatmp(:,:,1)
159                    etatmp(:,:,4) = c(i)*etatmp(:,:,2)
160
161                end if
162                do k = 1, 3*antithetics+1
163                    do t = 2, n
164                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),&
165                        1,0.0d0,alphatmp(:,t,k),1)
166                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
167                        1.0d0,alphatmp(:,t,k),1)
168                        call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,&
169                        alphatmp(:,t,k),1,0.0d0,sim(:,t,i+(k-1)*nsim),1)
170                    end do
171                end do
172            case(6)
173
174
175                alphatmp(:,1,1) = aplus(:,1)
176                etatmp(:,:,1) = etaplus(:,1:(n-1),i)
177
178                if(antithetics .EQ. 1) then
179                    alphatmp(:,1,2) = -aplus(:,1)
180                    alphatmp(:,1,3) = c(i)*alphatmp(:,1,1)
181                    alphatmp(:,1,4) = c(i)*alphatmp(:,1,2)
182
183                    etatmp(:,:,2) = -etaplus(:,1:(n-1),i)
184                    etatmp(:,:,3) = c(i)*etatmp(:,:,1)
185                    etatmp(:,:,4) = c(i)*etatmp(:,:,2)
186
187                end if
188                do k = 1, 3*antithetics+1
189                    do l = 1, p
190                          sim(l,1,i+(k-1)*nsim) = ddot(m,zt(l,:,1),1,alphatmp(:,1,k),1)
191
192                    end do
193                    do t = 2, n
194                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),&
195                        1,0.0d0,alphatmp(:,t,k),1)
196                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
197                        1.0d0,alphatmp(:,t,k),1)
198                        do l = 1, p
199
200                                sim(l,t,i+(k-1)*nsim) = ddot(m,zt(l,:,(t-1)*timevar(1)+1),1,alphatmp(:,t,k),1)
201
202                        end do
203                    end do
204                end do
205        end select
206    end do
207
208end subroutine simgaussianuncond
209
210