1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2c Copyright (C) INRIA
3c
4c Copyright (C) 2012 - 2016 - Scilab Enterprises
5c
6c This file is hereby licensed under the terms of the GNU GPL v2.0,
7c pursuant to article 5.3.4 of the CeCILL v.2.1.
8c This file was originally licensed under the terms of the CeCILL v2.1,
9c and continues to be available under such terms.
10c For more information, see the COPYING file which you should have received
11c along with this program.
12c
13C/MEMBR ADD NAME=ICSE2,SSI=0
14c
15      subroutine icse2(ind,nu,u,co,g,icsef,icsec2,icsei,y0,tob,obs,
16     &ob,ytot,f,b,fy,fu,ipv2,itob,cof,ytob,c2y,y0u,dmy,smy,oldmu,
17     &y,oldp,p,p0,gt,yob,d,itu,dtu,
18     &t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
19     &itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
20c
21c     sous programme de icse.fortran:calcul du gradient par
22c     integration du systeme adjoint
23c
24      implicit double precision (a-h,o-z)
25      dimension iu(5), u(nu),g(nu),y0(ny),tob(ntob),obs(nob,ny),
26     &ob(nex,ntob,nob),ytot(ny,nti+ntf),f(ny),b(ny),fy(ny,ny),
27     &fu(ny,nuc+nuv),ipv2(ny),itob(ntob),cof(nob,ntob),ytob(ny,ntob),
28     &c2y(ny,ntob),y0u(ny,nu),dmy(ny,ny),smy(ny,ny),
29     &oldmu(ny,nuc+nuv),y(ny),oldp(ny),p(ny),p0(ny),
30     &gt(nu),yob(nob,ntob),d(nob),itu(nitu),dtu(ndtu)
31      external icsef,icsec2,icsei
32c
33      character*6 nomf,nomc,nomi
34c
35c     Initialisation
36c
37      call dset(nu,0.0d+0,g,1)
38      call dset(ny,0.0d+0,p,1)
39      kt=nti+ntf
40      ktob=ntob
41c     lui et nui servent quand l'etat initial depend du controle
42      if (iu(2).gt.0) lui=min(nu,1+nuc)
43      if (iu(1).gt.0) lui=1
44      nui=iu(1)*nuc+iu(2)*nuv*(nti+ntf+1)
45c
46c
47c     Calcul de itob,vecteur des indices des instants de mesure
48c     a partir de tob,vecteur des instants de mesure
49c     itob(j) est l'entier le plus proche de tob(j)/dt
50c
51      do 1 j=1,ntobi
521     itob(j)=int(((tob(j)-t0)/dti)+0.50d+0)
53      if (ntobi.lt.ntob) then
54        itob(ntobi+1)=nti+int(0.50d+0+(tob(ntobi+1)-t0-nti*dti)/dtf)
55      endif
56      if (ntobi+1.lt.ntob) then
57        do 2 j=ntobi+2,ntob
582       itob(j)=itob(ntobi+1)+int(0.50d+0+(tob(j)-tob(ntobi+1))/dtf)
59      endif
60c
61c     Ecriture de ytob tableau des valeurs de l'etat
62c     aux instants de mesure
63c
64      do 10 j=1,ntob
65      do 10 i=1,ny
6610    call dcopy(ny,ytot(1,itob(j)),1,ytob(1,j),1)
67c
68c     Si ind=2,on calcule seulement le cout
69c     Si ind=3,on calcule seulement le gradient
70c     Si ind=4,on calcule le cout et le gradient
71c
72      if (ind.ne.3) then
73        indc=1
74        call icsec2(indc,nu,tob,obs,cof,ytob,ob,u,
75     &  co,c2y,g,yob,d,itu,dtu,
76     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
77     & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
78        if (indc.le.0) then
79          ind=indc
80          return
81        endif
82      endif
83      if (ind.eq.2) return
84c
85c     Calcul du gradient du cout en utilisant l'etat adjoint
86c     discretise:
87c     calcul de la derivee partielle c2y du cout par rapport
88c     a l'etat
89c
90      indc=2
91      call icsec2(indc,nu,tob,obs,cof,ytob,ob,u,co,c2y,g,yob,d,itu,
92     &dtu,
93     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
94     & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
95      if (indc.le.0) then
96        ind=indc
97        return
98      endif
99c
100c     +Evaluations successives de la contribution au gradient
101c     a chaque pas de temps a l'aide de l'etat adjoint(non stocke)
102c
103      do 100 kt=nti+ntf,1,-1
104c
105c     *Calcul de l'etat adjoint p_kt au pas kt
106c       Calcul du second membre p
107c       Initialisation:
108c       nt=nti+ntf
109c       si kt=nt,p est nul
110c       si kt<nt,on a au depart p=p_kt+1,etat adjoint au pas kt+1;
111c       on prend p=p0,ou p0=(smy)t*I*p avec smy=Id+(dt/2).dfy(t,y,u)
112c       avec dt=t_(kt+1)-t_kt,y=y_kt,t=t_kt,I designant la matrice
113c       diagonale d'ordre ny dont les nea premiers coefficients
114c       valent 0 et les autres 1 et Id designant la matrice
115c       identite d'ordre ny;
116c       si le systeme est affine avec partie lineaire autonome
117c       (ilin=2) smy est calculee seulement aux premiers pas de temps
118c       pour dti et dtf,sinon (ilin=0 ou 1) smy est calculee a chaque
119c       pas de temps
120c
121c       stockage de l'etat adjoint et de dt/2 au pas kt+1
122c
123        call dcopy(ny,p,1,oldp,1)
124        luv=min(nu,1+nuc+kt*nuv)
125c
126c       calcul de y=y_kt,dt=t_(kt+1)-t_kt,dt2=dt/2,
127c                 dt2new=(t_kt-t_(kt-1))/2
128c
129        call dcopy(ny,ytot(1,kt),1,y,1)
130c
131        if (kt.lt.nti) then
132          t=kt*dti+t0
133          dt=dti
134        else
135          t=nti*dti+(kt-nti)*dtf+t0
136          dt=dtf
137        endif
138        dt2=dt/2.0d+0
139        if (kt.ne.nti) then
140          dt2new=dt2
141        else
142          dt2new=dti/2.0d+0
143        endif
144c
145c       Dans le cas ilin<=1,
146c         calcul de fy=dfy(t,y,u) puis de smy=Id+(dt/2).dmy
147c         lorsque kt<(nti+ntf)
148c       Sinon (ilin>1),fy=dfy a ete calcule dans icse1
149c
150c
151        if (ilin.le.1) then
152          indf=2
153          call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu,
154     &        t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
155     &        itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
156          if (indf.le.0) then
157            ind=indf
158            return
159          endif
160        endif
161c
162        if (kt.ne.nti+ntf) then
163          if (ilin.le.1.or.kt.eq.nti+ntf-1.or.kt.eq.nti-1) then
164            do 30 i=1,ny
165            do 30 j=1,ny
16630          smy(i,j)=dt2*fy(i,j)
167            do 35 i=1,ny
16835          smy(i,i)=1.0d+0+smy(i,i)
169          endif
170c
171c         calcul de p0=(smy)t*I*p puis p=p0
172c
173          if (nea.gt.0) then
174            do 40 i=1,nea
17540          p(i)=0.0d+0
176          endif
177          call dmmul(p,1,smy,ny,p0,1,1,ny,ny)
178c
179          call dcopy(ny,p0,1,p,1)
180        endif
181c
182c       Fin du calcul du second membre p
183c         si kt=itob(ktob),on ajoute c2y(.,ktob) au second membre p
184c
185        if (ktob.gt.0) then
186          if (kt.eq.itob(ktob)) then
187            do 50 i=1,ny
18850          p(i)=p(i)+c2y(i,ktob)
189            ktob=ktob-1
190          endif
191        endif
192c
193c       Calcul et factorisation de la matrice dmy du systeme
194c       de l'etat adjoint
195c       dmy=I-dt2new.dfy(t,y,u) avec dt2new=(t_kt-t_(kt-1))/2,
196c       y=y_kt,t=t_kt,Idesignant la matrice diagonale d'ordre
197c       ny dont les nea premiers coefficients valent 0 et les
198c       autres 1;
199c       si le systeme est affine avec partie lineaire autonome
200c       (ilin=2) dmy est calculee et factorisee aux premiers
201c       pas de temps pour dti et dtf,sinon (ilin=0 ou 1) dmy est
202c       calculee et factorisee a chaque pas de temps
203c
204        if (ilin.le.1.or.kt.eq.nti+ntf.or.kt.eq.nti) then
205          do 60 i=1,ny
206          do 60 j=1,ny
20760        dmy(i,j)=-dt2new*fy(i,j)
208          do 65 i=1+nea,ny
20965        dmy(i,i)=1.0d+0+dmy(i,i)
210          call dgefa(dmy,ny,ny,ipv2,info)
211        endif
212c
213c       Resolution de (dmy)t*X=p,la solution s'appelant p
214c       p est alors p_kt,etat adjoint au pas kt
215c
216        call dgesl(dmy,ny,ny,ipv2,p,1)
217c
218c     *Calcul du gradient g au pas kt+1
219c       calcul de la contribution gt au gradient au pas kt+1:
220c       gt=(dt/2).(I*dfu(t_kt,y_kt,u)+dfu(t_kt+1,y_kt+1,u))t*p_kt+1
221c       avec dt=t_(kt+1)-t_kt,I designant la matrice diagonale
222c       d'ordre ny dont les nea premiers coefficients valent 0
223c       et les autres 1
224c
225        if (nuv.gt.0.or.iu(3).eq.1) then
226          indf=3
227          call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu,
228     &        t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
229     &        itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
230          if (indf.le.0) then
231            ind=indf
232            return
233          endif
234          if (kt.lt.nti+ntf) then
235            call dmmul(oldp,1,oldmu,ny,gt,1,1,ny,nuc+nuv)
236            call dscal(nuc+nuv,dt2,gt,1)
237c           le gradient g est la somme des contributions
238            if(iu(3).gt.0) then
239              call dadd(nuc,gt,1,g,1)
240            endif
241            if(nuv.gt.0) then
242              luv=min(nu,1+nuc+(kt+1)*nuv)
243              call dadd(nuv,gt(nuc+1),1,g(luv),1)
244            endif
245            if (nea.gt.0) then
246              do 70 i=1,nea
24770            oldp(i)=0.0d+0
248            endif
249            call dmmul(oldp,1,fu,ny,gt,1,1,ny,nuc+nuv)
250            call dscal(nuc+nuv,dt2,gt,1)
251c           le gradient g est la somme des contributions
252            if(iu(3).gt.0) then
253              call dadd(nuc,gt,1,g,1)
254            endif
255            if(nuv.gt.0) then
256              luv=min(nu,1+nuc+kt*nuv)
257              call dadd(nuv,gt(nuc+1),1,g(luv),1)
258            endif
259          endif
260c
261c         stockage de dfu(t_kt,y_kt,u) dans oldmu
262c
263        call dcopy(ny*(nuc+nuv),fu,1,oldmu,1)
264c
265c     *On passe au pas de temps suivant:kt-1,sauf si kt=1,auquel cas
266c       on calcule la contribution gt au gradient au pas kt=1 et
267c       on l'ajoute a g;on a:
268c       gt=(dt/2).(I*dfu(t0,y0,u)+dfu(t_1,y_1,u))t*p_1,avec
269c       dt=t_1-t0=dti car nti n'est jamais nul par convention
270c
271          if (kt.eq.1) then
272            t=t0
273            dt2=dti/2.0d+0
274            call dcopy(ny,y0,1,y,1)
275            indf=3
276            call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu,
277     &        t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
278     &        itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
279            if (indf.le.0) then
280              ind=indf
281              return
282            endif
283            call dmmul(p,1,oldmu,ny,gt,1,1,ny,nuc+nuv)
284            call dscal(nuc+nuv,dt2,gt,1)
285c           le gradient g est la somme des contributions
286            if(iu(3).gt.0) then
287              call dadd(nuc,gt,1,g,1)
288            endif
289            if(nuv.gt.0) then
290              luv=min(nu,1+nuc+nuv)
291              call dadd(nuv,gt(nuc+1),1,g(luv),1)
292            endif
293            if (nea.gt.0) then
294              do 90 i=1,nea
29590            p(i)=0.0d+0
296            endif
297            call dmmul(p,1,fu,ny,gt,1,1,ny,nuc+nuv)
298            call dscal(nuc+nuv,dt2,gt,1)
299c           le gradient g est la somme des contributions
300            if(iu(3).gt.0) then
301              call dadd(nuc,gt,1,g,1)
302            endif
303            if(nuv.gt.0) then
304              luv=min(nu,1+nuc)
305              call dadd(nuv,gt(nuc+1),1,g(luv),1)
306            endif
307          endif
308        endif
309100   continue
310c
311c     gradient par rapport au controle initial
312c
313      if(max(iu(1),iu(2)) .gt.0) then
314c
315c     calcul de l'etat adjoint initial p0
316c
317        indf=2
318        call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu,
319     &        t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
320     &        itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
321        if (indf.eq.0) then
322          ind=indf
323          return
324        endif
325        do 120 i=1,ny
326        do 120 j=1,ny
327120     smy(i,j)=dt2*fy(i,j)
328        do 125 i=1,ny
329125     smy(i,i)=1.0d+0+smy(i,i)
330        if (nea.gt.0) then
331          do 130 i=1,nea
332130       p(i)=0.0d+0
333        endif
334        call dmmul(p,1,smy,ny,p0,1,1,ny,ny)
335c     incrementation du gradient
336        indi=2
337        call icsei(indi,nui,u(lui),y0,y0u,itu,dtu,
338     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
339     & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi)
340        if (indi.le.0) then
341           ind=indi
342           return
343        endif
344        call dmmul(p0,1,y0u,ny,gt,1,1,nui,nui)
345        do 150 i=1,nui
346150     g(lui+i-1)=g(lui+i-1)+gt(i)
347c
348      endif
349      end
350