3ad847c4e2638ac4cd0f352deb8bdea332c01c3d
[scilab.git] / scilab / modules / optimization / src / fortran / zgcbd.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) INRIA
3 c
4 c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 c
6 c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 c This file was originally licensed under the terms of the CeCILL v2.1,
9 c and continues to be available under such terms.
10 c For more information, see the COPYING file which you should have received
11 c along with this program.
12 c
13       subroutine zgcbd(simul,n,binf,bsup,x,f,g,zero,napmax,itmax,indgc
14      &  ,ibloc,nfac,iprint,io,epsx,epsf,epsg,dir,df0,diag,x2,
15      &izs,rzs,dzs,y,s,z,ys,zs,nt,index,wk1,wk2,alg,ialg,nomf)
16 c
17       implicit double precision (a-h,o-z)
18       real rzs(*)
19       double precision dzs(*)
20       dimension x2(n),dir(n),epsx(n)
21       dimension binf(n),bsup(n),x(n),g(n),diag(n),ibloc(n),izs(*)
22       dimension y(nt,n),s(nt,n),z(nt,n),ys(nt),zs(nt)
23       dimension wk1(n),wk2(n),alg(15)
24       character*6 nomf
25       integer index(nt),ialg(15)
26       external simul
27       character bufstr*(4096)
28 c
29       if(iprint.ge.4) then
30         write(bufstr,10000)
31         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
32 10000 format (' dans gcbd. algorithme utilise: ')
33       if(ialg(1).eq.1) then
34         write(bufstr,10001)
35         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
36         endif
37 10001 format ('        emploi correction de powell ')
38       if(ialg(2).eq.1) then
39         write(bufstr,10002)
40         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
41         endif
42 10002 format ('  mise a jour de diag par la methode bfgs')
43       if(ialg(3).eq.1) then
44         write(bufstr,10003)
45         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
46         endif
47 10003 format ('  mise a echelle de diag par methode de shanno-phua')
48       if(ialg(3).eq.2) then
49         write(bufstr,10004)
50         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
51         endif
52 10004 format ('  mise a echelle de diag seulement a la 2e iter')
53       if(ialg(4).eq.1) then
54         write(bufstr,10005)
55         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
56         endif
57 10005 format ('      memorisation pour choix iteration ')
58       if(ialg(5).eq.1) then
59         write(bufstr,10006)
60         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
61         endif
62 10006 format ('      memorisation par variable')
63       if(ialg(6).eq.1) then
64         write(bufstr,10007)
65         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
66         endif
67 10007 format ('      relachememt de variables a toutes les iteration')
68       if(ialg(6).eq.2) then
69         write(bufstr,10008)
70         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
71         endif
72 10008 format ('      relachement de vars si decroissance g_norme')
73       if(ialg(6).eq.10) then
74         write(bufstr,10009)
75         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
76         endif
77 10009 format ('      relachement de vars si dec f % iter_init du cycle')
78       if(ialg(6).eq.11) then
79         write(bufstr,10010)
80         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
81         endif
82 10010 format ('      relachement de vars si dec f % dec du cycle')
83       if(ialg(7).eq.1) then
84         write(bufstr,10011)
85         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
86         endif
87 10011 format ('      choix de vars a relacher par bertsekas modifiee')
88       if(ialg(8).eq.1) then
89         write(bufstr,10012)
90         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
91         endif
92 10012 format ('      choix de dir descente par methode de gradient')
93       if(ialg(8).eq.2) then
94         write(bufstr,10013)
95         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
96         endif
97 10013 format ('      choix de dir descente par methode qn')
98       if(ialg(8).eq.3) then
99         write(bufstr,10014)
100         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
101         endif
102 10014 format ('      choix de dir descente par qn sans memoire.nt depl')
103       if(ialg(8).eq.4) then
104         write(bufstr,10015)
105         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
106         endif
107 10015 format ('      choix de dir descente par qn -mem,redem,sans acc.')
108       if(ialg(8).eq.5) then
109         write(bufstr,10016)
110         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
111         endif
112 10016 format ('     choix de dir descente par qn -mem,redem,avec acc.')
113       if(ialg(9).eq.2) then
114         write(bufstr,10017)
115         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
116         endif
117 10017 format ('      redem si relachement de vars')
118       if(ialg(9).eq.10) then
119         write(bufstr,10018)
120         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
121         endif
122 10018 format ('      redem si dec f % dec iter_init du cycle')
123       if(ialg(9).eq.11) then
124         write(bufstr,10019)
125         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
126         endif
127 10019 format ('      redem si dec f % dec totale du cycle.')
128       if(ialg(9).eq.12) then
129         write(bufstr,10020)alg(9)
130         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
131         endif
132 10020 format ('    redem si diminution du gradient des var libres d un',
133      & 'facteur',d11.4)
134       endif
135 c
136 c     section 1  initialisations
137 c     irl nombre de rech lin 'lentes'
138 c     nred nombre de redemarrage de la direction de descente
139 c     icycl nombre de cycles de minimisation
140 c
141       epsgcp=1.0d-5
142       indsim=4
143       indrl=1
144       irl=0
145       irl=0
146       nred=1
147       icycl=1
148       nap=0
149 c
150       iresul=1
151       call proj(n,binf,bsup,x)
152       indsim=4
153       call simul(indsim,n,x,f,g,izs,rzs,dzs)
154       nap=nap+1
155       if(indsim.gt.0)go to 99
156       indgc=-1
157       if(indsim.eq.0)indgc=0
158       if(iprint.gt.0) then
159         write(bufstr,123)indgc
160         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
161         endif
162       go to 900
163 99    continue
164       ceps0=20.0d+0
165       eps0=0.0d+0
166       do 100 i=1,n
167 100   eps0=eps0+epsx(i)
168       eps0=ceps0*eps0/n
169 c
170 c     calcul de zng
171       znog0=rednor(n,binf,bsup,x,epsx,g)
172       zng=znog0
173       zngrit=znog0
174       zngred=znog0
175 c
176       do 130 i=1,n
177 130   ibloc(i)=0
178       izag=3
179       izag1=izag
180       nap=0
181       iter=0
182       scal=1.0d+0
183       nfac=n
184       np=0
185       lb=1
186       nb=2
187       if(ialg(8).eq.3) nb=1
188       do 140 i=1,nt
189 140   index(i)=i
190       tetaq=alg(9)
191       condm=alg(2)
192       param=alg(1)
193       indgc1=indgc
194 c     si indgc=0 on init diag a k*ident puis scal a it=2
195 c
196       if(indgc.eq.1.or.indgc.ge.100)go to 150
197       if(indgc.eq.2)go to 180
198       indgc=-13
199       if(iprint.gt.0) then
200         write(bufstr,123) indgc
201         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
202         endif
203       go to 900
204 c
205 150   continue
206 c     on initialise diag par approximation quadratique
207 c     df0 decroissance prevue . si mod quad df0=((dh)-1g,g)/2
208 c     et on cherche dh diag de la forme cst/(dx)**2
209 c     donc cst=som((g(i)*(dx))**2))/(2*df0)
210       sy=0.0d+0
211       do 160 i=1,n
212 160    sy=sy+(g(i)*epsx(i))**2
213       sy=sy/(2.0d+0*df0)
214       do 170 i=1,n
215 170   diag(i)=(sy + zero)/(epsx(i)**2 + zero)
216 180   continue
217 c
218 c
219 c     bouclage
220 200   iter=iter +1
221       indgc=1
222       if(iter.gt.itmax)then
223          indgc=5
224          go to 900
225       endif
226 201   continue
227       if(iprint.ge.2) then
228         write(bufstr,1210)iter,f
229         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
230         endif
231 1210  format(' dans gcbd  iter=',i3,'  f=',d15.7)
232       if(iter.eq.1)then
233          irit=1
234          goto 301
235       endif
236 c
237       call majysa(n,nt,np,y,s,ys,lb,g,x,wk2,wk1,index,ialg,nb)
238       inp=index(np)
239 c
240 c
241 c     correction powell sur y si (y,s) trop petit
242       if(ialg(1).ne.1) go to 290
243       param1=1.-param
244       bss=0.0d+0
245       do 260 i=1,n
246 260   bss=bss + diag(i)*s(inp,i)**2
247       bss2=param*bss
248       if(ys(inp).gt.bss2)go to 290
249       if(iprint.gt.2) then
250         write(bufstr,1270) ys(inp)
251         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
252         endif
253 1270  format(' gcbd. emploi correction powell (y,s)=',d11.4)
254       teta=param1*bss/(bss-ys(inp))
255       teta1=1.0d+0-teta
256       do 274 i=1,n
257 274   y(inp,i)=teta*y(inp,i)+teta1*diag(i)*s(inp,i)
258       ys(inp)=bss2
259 c     verif correction powell (facultatif; faire go to 300)
260       ys1=ddot(n,s(inp,1),1,y(inp,1),1)
261       ys1=abs(bss2-ys1)/bss2
262       if(iprint.gt.2) then
263         write(bufstr,1280) ys1
264         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
265         endif
266 1280  format(' erreur relative correction powell =',d11.4)
267 c
268 c mise a jour de diag
269 290   continue
270       if(ialg(2).eq.1)
271      &  call bfgsd(diag,n,nt,np,y,s,ys,condm,param,zero,index)
272 c
273       if(ialg(3).eq.1.or.(ialg(3).eq.2.and.iter.eq.2))
274      &  call shanph(diag,n,nt,np,y,s,ys,scal,index,io,iprint)
275 c
276       call majz(n,np,nt,y,s,z,ys,zs,diag,index)
277 c
278 c     section 3 determination des variables libres et bloquees
279 300   continue
280 c     -----decision de relachement a l'iteration courante
281 c          relachement si irit=1 (sinon irit=0)
282       irit=0
283       if(ialg(6).eq.1) irit=1
284       if(ialg(6).eq.2.and.znglib.le.alg(6)*zngrit)irit=1
285       if(ialg(6).eq.10.and.diff.le.dfrit1*alg(6))irit=1
286       if(ialg(6).eq.11.and.diff.le.difrit*alg(6))irit=1
287       if(irit.eq.1) nred=nred+1
288 c    ----choix des variables a relacher
289       iprint1=iprint
290 301   if(ialg(7).eq.1)call relvar(ind,n,x,binf,bsup,x2,g,diag,
291      &  iprint,io,ibloc,izag,iter,nfac,irit)
292 c
293 c
294 c     section 4 expression de dir
295       if (np.eq.0) then
296          do 400 i=1,n
297            dir(i)=-g(i)/diag(i)
298 400      continue
299       else
300          do 410 i=1,n
301             dir(i)=-scal*g(i)
302 410      continue
303          call gcp(n,index,ibloc,np,nt,y,s,z,ys,zs,diag,g,dir,wk1,
304      &            wk2,epsgcp)
305       endif
306 c
307 c     section 5  redemarrage
308 c
309       if(ialg(8).eq.4.or.ialg(8).eq.5) then
310          ired=0
311          if(ialg(9).eq.2.and.ind.eq.1) ired=1
312          if(ialg(9).eq.10.and.diff.lt.dfred1*tetaq) ired=1
313          if(ialg(9).eq.11.and.diff.lt.difred*tetaq)  ired=1
314          if(ialg(9).eq.12.and.znglib.le.tetaq*zngred) ired=1
315          if(ired.eq.1) then
316             icycl=icycl+1
317             np=0
318             lb=1
319             if(iprint.gt.2) then
320               write(bufstr,1000) icycl
321               call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
322               endif
323 1000        format ('   redemarrage. icycl=',i5)
324          endif
325       endif
326 c
327 c     section 6 annulation de d(i) , i dans ib
328       if(ialg(6).eq.1)go to 640
329       do 630 i=1,n
330 630      if(ibloc(i).gt.0) dir(i)=0.0d+0
331 640   continue
332 c
333 c     recherche lineaire
334 c     conservation de x et g dans wk1 et wk2
335       call dcopy(n,x,1,wk1,1)
336       call dcopy(n,g,1,wk2,1)
337 c     calcul de la derivee dans la direction dir
338       ifp=0
339       fn=f
340       znog0=zng
341 702   dfp=0.0d+0
342       do 710 i=1,n
343       epsxi=epsx(i)
344       xi=x(i)
345       diri=dir(i)
346       if(xi-binf(i).le.epsxi.and.diri.lt.0.0d+0)dir(i)=0.0d+0
347 710   if(bsup(i)-xi.le.epsxi.and.diri.gt.0.0d+0)dir(i)=0.0d+0
348       dfp=ddot(n,g,1,dir,1)
349       if(-dfp.gt.0)go to 715
350       if(ifp.eq.1) then
351          indgc=6
352          go to 900
353       endif
354 c     restauration dir
355       if(iprint.ge.3) then
356          write(bufstr,1712)dfp,zero
357          call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
358          endif
359 1712  format(' gcbd : restauration dir ; fp,zero',2d11.4)
360       do 712 i=1,n
361 712   dir(i)=-scal*g(i)
362       ifp=1
363       go to 702
364 715   continue
365 c     pas initial suivant idee fletcher
366       t=-2.0d+0*diff/dfp
367       if(iter.eq.1)t=-2.0d+0*df0/dfp
368       tmax=1.0d+10
369       t=min(t,tmax)
370       t=max(t,1.0d+10*zero)
371       napm=15
372       napm1=nap + napm
373       if(napm1.gt.napmax) napm1=napmax
374       napav=nap
375       amd=0.70d+0
376       amf=0.10d+0
377 c
378       call rlbd(indrl,n,simul,x,binf,bsup,f,dfp,t,tmax,dir,g,tproj,
379      &  amd,amf,iprint,io,zero,nap,napm1,x2,izs,rzs,dzs)
380       if(iprint.gt.2) then
381         write(bufstr,750)indrl,t,f
382         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
383         endif
384 750   format(' retour mlibd indrl=',i6,' pas= ',d11.4,' f= ',d11.4)
385       if(nap-napav.ge.5) irl=irl+1
386       if(indrl.ge.10)then
387          indsim=4
388          nap=nap + 1
389          call simul(indsim,n,x,f,g,izs,rzs,dzs)
390          if(indsim.le.0)then
391             indgc=-3
392             if(indsim.eq.0)indgc=0
393             if(iprint.gt.0) then
394               write(bufstr,123)indgc
395               call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
396               endif
397 123         format(' gcbd : retour avec indgc=',i8)
398             go to 900
399          endif
400       endif
401       if(indrl.le.0)then
402          indgc=10
403          if(indrl.eq.0)indgc=0
404          if(indrl.eq.-3)indgc=13
405          if(indrl.eq.-4)indgc=12
406          if(indrl.le.-1000)indgc=11
407          if(iprint.gt.0) then
408            write(bufstr,123)indgc
409            call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
410            endif
411          go to 900
412       endif
413       if(iprint.ge.5) then
414          do 760 i=1,n
415 760      if(iprint.gt.2) write(io,777)i,x(i),g(i),dir(i)
416 777      format(' i=',i2,' xgd ',3f11.4)
417
418       endif
419 c
420       if(nap.lt.napmax)go to 758
421       if(iprint.gt.0)  then
422         write(bufstr,755)
423         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
424         endif
425 755   format(' gcbd max appels simul')
426       indgc=4
427       go to 900
428 758   continue
429 c
430 c     section 8 test de convergence
431       do 805 i=1,n
432       if(abs(x(i)-wk1(i)).gt.epsx(i))go to 806
433 805   continue
434       if(iprint.gt.0) then
435         write(bufstr,1805)
436         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
437         endif
438 1805  format(' gcbd. retour apres convergence sur x')
439       indgc=3
440       go to 900
441 c     calcul grad residuel,norme l2
442 806   continue
443       difg=rednor(n,binf,bsup,x,epsx,g)
444       diff=fn-f
445       if(iprint.ge.2) then
446         write(bufstr,860)epsg,difg,epsf,diff,nap
447         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
448         endif
449 860   format(' gcbd. epsg,difg=',2d11.4,'  epsf,diff=',2d11.4
450      &,'  nap=',i3)
451 c
452       if(diff.le.epsf) then
453          indgc=2
454          go to 900
455       endif
456       if(difg.le.epsg) then
457          indgc=1
458          go to 900
459       endif
460 c
461 c     -----mise a jour de difrit,dfrit1,difred,dfred1
462       if(irit.eq.1) then
463          difrit=diff
464          dfrit1=diff
465        else
466          difrit=difrit+diff
467       endif
468       if(ired.eq.1) then
469          difred=diff
470          dfred1=diff
471        else
472          difred=difred + diff
473       endif
474 c
475       znglib=0.0d+0
476       do 884 i=1,n
477          if(ibloc(i).gt.0)go to 884
478          aa=g(i)
479          if(x(i)-binf(i).le.epsx(i)) aa=min(0.0d+0,aa)
480          if(bsup(i)-x(i).le.epsx(i)) aa=max(0.0d+0,aa)
481          znglib=znglib+aa**2
482 884   continue
483       znglib=sqrt(znglib)
484       if(ired.eq.1)zngred=znglib
485       if(irit.eq.1)zngrit=znglib
486       go to 200
487 c
488 c      fin des calculs
489 900   if(indrl.eq.0)indgc=0
490       if(indgc.eq.1.and.indrl.le.0)  indgc=indrl
491       if(iprint.gt.0) then
492          write(bufstr,123)indgc
493          call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
494          endif
495       if(iprint.ge.1.and.indrl.le.zero) then
496         write(bufstr,1910)
497         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
498         write(bufstr,1911) indrl
499         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
500         endif
501 1910  format(' arret impose par la recherche lineaire. cf notice rlbd')
502 1911  format(' indicateur de rlbd=',i6)
503       if(iprint.ge.1) then
504         write(bufstr,950)f,difg,nap,iter,indgc
505         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
506         endif
507 950   format(' f,norme grad,nap,iter,indgc=',2e11.4,3i6)
508 c
509 c     autres impressions finales
510       if(indgc1.lt.100) return
511       zrl=0.
512       if(iter.gt.0) zrl=dble(nap)/dble(iter)
513 2000  format('     nom    n       f        norm2g   nf   iter  rl/it ',
514      &          ' irl   cpu  cycl   red')
515
516       write(bufstr,2001)nomf,f,difg,nap,iter,zrl,irl
517       call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
518 2001  format(1x,a6,2e11.4,2i5,f6.2,i5)
519       end