[optimization] windows link fixed && valgrind error fixed
[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       znglib=0.0d+0
150 c
151       iresul=1
152       call proj(n,binf,bsup,x)
153       indsim=4
154       call simul(indsim,n,x,f,g,izs,rzs,dzs)
155       nap=nap+1
156       if(indsim.gt.0)go to 99
157       indgc=-1
158       if(indsim.eq.0)indgc=0
159       if(iprint.gt.0) then
160         write(bufstr,123)indgc
161         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
162         endif
163       go to 900
164 99    continue
165       ceps0=20.0d+0
166       eps0=0.0d+0
167       do 100 i=1,n
168 100   eps0=eps0+epsx(i)
169       eps0=ceps0*eps0/n
170 c
171 c     calcul de zng
172       znog0=rednor(n,binf,bsup,x,epsx,g)
173       zng=znog0
174       zngrit=znog0
175       zngred=znog0
176 c
177       do 130 i=1,n
178 130   ibloc(i)=0
179       izag=3
180       izag1=izag
181       nap=0
182       iter=0
183       scal=1.0d+0
184       nfac=n
185       np=0
186       lb=1
187       nb=2
188       if(ialg(8).eq.3) nb=1
189       do 140 i=1,nt
190 140   index(i)=i
191       tetaq=alg(9)
192       condm=alg(2)
193       param=alg(1)
194       indgc1=indgc
195 c     si indgc=0 on init diag a k*ident puis scal a it=2
196 c
197       if(indgc.eq.1.or.indgc.ge.100)go to 150
198       if(indgc.eq.2)go to 180
199       indgc=-13
200       if(iprint.gt.0) then
201         write(bufstr,123) indgc
202         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
203         endif
204       go to 900
205 c
206 150   continue
207 c     on initialise diag par approximation quadratique
208 c     df0 decroissance prevue . si mod quad df0=((dh)-1g,g)/2
209 c     et on cherche dh diag de la forme cst/(dx)**2
210 c     donc cst=som((g(i)*(dx))**2))/(2*df0)
211       sy=0.0d+0
212       do 160 i=1,n
213 160    sy=sy+(g(i)*epsx(i))**2
214       sy=sy/(2.0d+0*df0)
215       do 170 i=1,n
216 170   diag(i)=(sy + zero)/(epsx(i)**2 + zero)
217 180   continue
218 c
219 c
220 c     bouclage
221 200   iter=iter +1
222       indgc=1
223       if(iter.gt.itmax)then
224          indgc=5
225          go to 900
226       endif
227 201   continue
228       if(iprint.ge.2) then
229         write(bufstr,1210)iter,f
230         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
231         endif
232 1210  format(' dans gcbd  iter=',i3,'  f=',d15.7)
233       if(iter.eq.1)then
234          irit=1
235          goto 301
236       endif
237 c
238       call majysa(n,nt,np,y,s,ys,lb,g,x,wk2,wk1,index,ialg,nb)
239       inp=index(np)
240 c
241 c
242 c     correction powell sur y si (y,s) trop petit
243       if(ialg(1).ne.1) go to 290
244       param1=1.-param
245       bss=0.0d+0
246       do 260 i=1,n
247 260   bss=bss + diag(i)*s(inp,i)**2
248       bss2=param*bss
249       if(ys(inp).gt.bss2)go to 290
250       if(iprint.gt.2) then
251         write(bufstr,1270) ys(inp)
252         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
253         endif
254 1270  format(' gcbd. emploi correction powell (y,s)=',d11.4)
255       teta=param1*bss/(bss-ys(inp))
256       teta1=1.0d+0-teta
257       do 274 i=1,n
258 274   y(inp,i)=teta*y(inp,i)+teta1*diag(i)*s(inp,i)
259       ys(inp)=bss2
260 c     verif correction powell (facultatif; faire go to 300)
261       ys1=ddot(n,s(inp,1),1,y(inp,1),1)
262       ys1=abs(bss2-ys1)/bss2
263       if(iprint.gt.2) then
264         write(bufstr,1280) ys1
265         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
266         endif
267 1280  format(' erreur relative correction powell =',d11.4)
268 c
269 c mise a jour de diag
270 290   continue
271       if(ialg(2).eq.1)
272      &  call bfgsd(diag,n,nt,np,y,s,ys,condm,param,zero,index)
273 c
274       if(ialg(3).eq.1.or.(ialg(3).eq.2.and.iter.eq.2))
275      &  call shanph(diag,n,nt,np,y,s,ys,scal,index,io,iprint)
276 c
277       call majz(n,np,nt,y,s,z,ys,zs,diag,index)
278 c
279 c     section 3 determination des variables libres et bloquees
280 300   continue
281 c     -----decision de relachement a l'iteration courante
282 c          relachement si irit=1 (sinon irit=0)
283       irit=0
284       if(ialg(6).eq.1) irit=1
285       if(ialg(6).eq.2.and.znglib.le.alg(6)*zngrit)irit=1
286       if(ialg(6).eq.10.and.diff.le.dfrit1*alg(6))irit=1
287       if(ialg(6).eq.11.and.diff.le.difrit*alg(6))irit=1
288       if(irit.eq.1) nred=nred+1
289 c    ----choix des variables a relacher
290       iprint1=iprint
291 301   if(ialg(7).eq.1)call relvar(ind,n,x,binf,bsup,x2,g,diag,
292      &  iprint,io,ibloc,izag,iter,nfac,irit)
293 c
294 c
295 c     section 4 expression de dir
296       if (np.eq.0) then
297          do 400 i=1,n
298            dir(i)=-g(i)/diag(i)
299 400      continue
300       else
301          do 410 i=1,n
302             dir(i)=-scal*g(i)
303 410      continue
304          call gcp(n,index,ibloc,np,nt,y,s,z,ys,zs,diag,g,dir,wk1,
305      &            wk2,epsgcp)
306       endif
307 c
308 c     section 5  redemarrage
309 c
310       if(ialg(8).eq.4.or.ialg(8).eq.5) then
311          ired=0
312          if(ialg(9).eq.2.and.ind.eq.1) ired=1
313          if(ialg(9).eq.10.and.diff.lt.dfred1*tetaq) ired=1
314          if(ialg(9).eq.11.and.diff.lt.difred*tetaq)  ired=1
315          if(ialg(9).eq.12.and.znglib.le.tetaq*zngred) ired=1
316          if(ired.eq.1) then
317             icycl=icycl+1
318             np=0
319             lb=1
320             if(iprint.gt.2) then
321               write(bufstr,1000) icycl
322               call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
323               endif
324 1000        format ('   redemarrage. icycl=',i5)
325          endif
326       endif
327 c
328 c     section 6 annulation de d(i) , i dans ib
329       if(ialg(6).eq.1)go to 640
330       do 630 i=1,n
331 630      if(ibloc(i).gt.0) dir(i)=0.0d+0
332 640   continue
333 c
334 c     recherche lineaire
335 c     conservation de x et g dans wk1 et wk2
336       call dcopy(n,x,1,wk1,1)
337       call dcopy(n,g,1,wk2,1)
338 c     calcul de la derivee dans la direction dir
339       ifp=0
340       fn=f
341       znog0=zng
342 702   dfp=0.0d+0
343       do 710 i=1,n
344       epsxi=epsx(i)
345       xi=x(i)
346       diri=dir(i)
347       if(xi-binf(i).le.epsxi.and.diri.lt.0.0d+0)dir(i)=0.0d+0
348 710   if(bsup(i)-xi.le.epsxi.and.diri.gt.0.0d+0)dir(i)=0.0d+0
349       dfp=ddot(n,g,1,dir,1)
350       if(-dfp.gt.0)go to 715
351       if(ifp.eq.1) then
352          indgc=6
353          go to 900
354       endif
355 c     restauration dir
356       if(iprint.ge.3) then
357          write(bufstr,1712)dfp,zero
358          call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
359          endif
360 1712  format(' gcbd : restauration dir ; fp,zero',2d11.4)
361       do 712 i=1,n
362 712   dir(i)=-scal*g(i)
363       ifp=1
364       go to 702
365 715   continue
366 c     pas initial suivant idee fletcher
367       t=-2.0d+0*diff/dfp
368       if(iter.eq.1)t=-2.0d+0*df0/dfp
369       tmax=1.0d+10
370       t=min(t,tmax)
371       t=max(t,1.0d+10*zero)
372       napm=15
373       napm1=nap + napm
374       if(napm1.gt.napmax) napm1=napmax
375       napav=nap
376       amd=0.70d+0
377       amf=0.10d+0
378 c
379       call rlbd(indrl,n,simul,x,binf,bsup,f,dfp,t,tmax,dir,g,tproj,
380      &  amd,amf,iprint,io,zero,nap,napm1,x2,izs,rzs,dzs)
381       if(iprint.gt.2) then
382         write(bufstr,750)indrl,t,f
383         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
384         endif
385 750   format(' retour mlibd indrl=',i6,' pas= ',d11.4,' f= ',d11.4)
386       if(nap-napav.ge.5) irl=irl+1
387       if(indrl.ge.10)then
388          indsim=4
389          nap=nap + 1
390          call simul(indsim,n,x,f,g,izs,rzs,dzs)
391          if(indsim.le.0)then
392             indgc=-3
393             if(indsim.eq.0)indgc=0
394             if(iprint.gt.0) then
395               write(bufstr,123)indgc
396               call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
397               endif
398 123         format(' gcbd : retour avec indgc=',i8)
399             go to 900
400          endif
401       endif
402       if(indrl.le.0)then
403          indgc=10
404          if(indrl.eq.0)indgc=0
405          if(indrl.eq.-3)indgc=13
406          if(indrl.eq.-4)indgc=12
407          if(indrl.le.-1000)indgc=11
408          if(iprint.gt.0) then
409            write(bufstr,123)indgc
410            call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
411            endif
412          go to 900
413       endif
414       if(iprint.ge.5) then
415          do 760 i=1,n
416 760      if(iprint.gt.2) write(io,777)i,x(i),g(i),dir(i)
417 777      format(' i=',i2,' xgd ',3f11.4)
418
419       endif
420 c
421       if(nap.lt.napmax)go to 758
422       if(iprint.gt.0)  then
423         write(bufstr,755)
424         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
425         endif
426 755   format(' gcbd max appels simul')
427       indgc=4
428       go to 900
429 758   continue
430 c
431 c     section 8 test de convergence
432       do 805 i=1,n
433       if(abs(x(i)-wk1(i)).gt.epsx(i))go to 806
434 805   continue
435       if(iprint.gt.0) then
436         write(bufstr,1805)
437         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
438         endif
439 1805  format(' gcbd. retour apres convergence sur x')
440       indgc=3
441       go to 900
442 c     calcul grad residuel,norme l2
443 806   continue
444       difg=rednor(n,binf,bsup,x,epsx,g)
445       diff=fn-f
446       if(iprint.ge.2) then
447         write(bufstr,860)epsg,difg,epsf,diff,nap
448         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
449         endif
450 860   format(' gcbd. epsg,difg=',2d11.4,'  epsf,diff=',2d11.4
451      &,'  nap=',i3)
452 c
453       if(diff.le.epsf) then
454          indgc=2
455          go to 900
456       endif
457       if(difg.le.epsg) then
458          indgc=1
459          go to 900
460       endif
461 c
462 c     -----mise a jour de difrit,dfrit1,difred,dfred1
463       if(irit.eq.1) then
464          difrit=diff
465          dfrit1=diff
466        else
467          difrit=difrit+diff
468       endif
469       if(ired.eq.1) then
470          difred=diff
471          dfred1=diff
472        else
473          difred=difred + diff
474       endif
475 c
476       znglib=0.0d+0
477       do 884 i=1,n
478          if(ibloc(i).gt.0)go to 884
479          aa=g(i)
480          if(x(i)-binf(i).le.epsx(i)) aa=min(0.0d+0,aa)
481          if(bsup(i)-x(i).le.epsx(i)) aa=max(0.0d+0,aa)
482          znglib=znglib+aa**2
483 884   continue
484       znglib=sqrt(znglib)
485       if(ired.eq.1)zngred=znglib
486       if(irit.eq.1)zngrit=znglib
487       go to 200
488 c
489 c      fin des calculs
490 900   if(indrl.eq.0)indgc=0
491       if(indgc.eq.1.and.indrl.le.0)  indgc=indrl
492       if(iprint.gt.0) then
493          write(bufstr,123)indgc
494          call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
495          endif
496       if(iprint.ge.1.and.indrl.le.zero) then
497         write(bufstr,1910)
498         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
499         write(bufstr,1911) indrl
500         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
501         endif
502 1910  format(' arret impose par la recherche lineaire. cf notice rlbd')
503 1911  format(' indicateur de rlbd=',i6)
504       if(iprint.ge.1) then
505         write(bufstr,950)f,difg,nap,iter,indgc
506         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
507         endif
508 950   format(' f,norme grad,nap,iter,indgc=',2e11.4,3i6)
509 c
510 c     autres impressions finales
511       if(indgc1.lt.100) return
512       zrl=0.
513       if(iter.gt.0) zrl=dble(nap)/dble(iter)
514 2000  format('     nom    n       f        norm2g   nf   iter  rl/it ',
515      &          ' irl   cpu  cycl   red')
516
517       write(bufstr,2001)nomf,f,difg,nap,iter,zrl,irl
518       call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
519 2001  format(1x,a6,2e11.4,2i5,f6.2,i5)
520       end