442672c8cc7ff96735f973558e0a834d2a1a1897
[scilab.git] / scilab / modules / optimization / sci_gateway / fortran / sci_f_optim.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) INRIA
3
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at    
8 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9 c
10 c     --------------------------
11 c     SCILAB function : optim
12 c     --------------------------
13       subroutine scioptim(fname)
14       
15       character*(*) fname
16         INCLUDE 'stack.h'
17       double precision zero,df0,zng,dxmin
18       double precision epsg,epsg1,epsf,dzs
19       integer top2,topin,topind,topx,top3
20       character*80   nomsub
21       common /optim/ nomsub
22       integer       nizs,nrzs,ndzs
23       common /nird/ nizs,nrzs,ndzs
24       external foptim,boptim,fuclid,ctonb,ctcab
25       integer coin,coar,coti,cotd,cosi,cosd,nfac
26 c     
27       character*(nlgh+1) namef,namej
28       common/csolve/namef,namej
29 c
30       integer impn(nsiz)
31       logical eqid, getscalar
32       integer iadr, sadr
33 c     
34       parameter (nsiz1=nsiz-1)
35       data impn/672732690,nsiz1*673720360/
36 c
37       data coin,coar,coti,cotd,cosi,cosd,nfac
38      &     /   5906,6922,4637,3357,4636,3356, 0/
39 c     
40       
41       iadr(l)=l+l-1
42       sadr(l)=(l/2)+1
43
44       napm=100
45       itmax=100
46       epsf=0.0d+0
47       iepsx=1
48       indtes=0
49       imp=0
50       io=wte
51       zero=stk(leps)
52       df0=1.0d+0
53       immx=0
54       epsg=zero
55       indtv=0
56       icontr=1
57       ialg=1
58       irecd=0
59       ireci=0
60
61 c
62       if(infstk(top).eq.1) then
63          infstk(top)=0
64          if(eqid(idstk(1,top),impn)) then
65             if (.not.getscalar('optim',top,top,lr)) return
66             imp=stk(lr)
67             top=top-1
68             rhs=rhs-1
69          endif
70       endif
71
72       lf=lstk(top+1)
73       ldisp=lf+1
74       lizs=iadr(ldisp)
75       lrzs=ldisp
76       ldzs=ldisp
77       ldisp=ldzs+1
78 c     
79       top2=top+1-rhs
80       top3=top2
81       topin=top2-1
82       il=iadr(lstk(top2))
83       if (rhs.lt.2) then
84          call error(39)
85          return
86       endif
87 c     
88 c     
89 c     traitement de simul
90       l1=istk(il)
91       if((l1-10)*(l1-11)*(l1-13)*(l1-15).ne.0) then
92          err=top2-topin
93          call error(80)
94          return
95       endif
96 c
97 c Compute isim :
98 c - isim=1 if the cost function is given as a string (code 10), that is, the cost function is 
99 c   a C or Fortran routine.
100 c   In that case, the cost function is computed with a call to "foptim".
101 c - isim=2 if the cost function is given as a script (code 11), a compiled script (code 13),
102 c   or a list (code 15)
103 c   In that case, the cost function is computed with a call to "boptim".
104 c
105 c     cas simul=liste
106       if(istk(il).eq.10) then
107          if (istk(il+1)*istk(il+2).ne.1) then
108             err=top2-topin
109             call error(89)
110             return
111          endif
112          nc=min(istk(il+5)-1,80)
113          nomsub(1:80)= ' '
114          call cvstr(nc,istk(il+6),nomsub,1)
115          if(err.gt.0) return
116          nomsub(nc+1:nc+1)=char(0)
117          call setfoptim(nomsub,irep)
118          if ( irep.eq.1) then 
119             buf = nomsub
120             call error(50)
121             return
122          endif
123          isim=1
124       endif
125 c     cas simul=macro
126       if(istk(il).eq.11.or.istk(il).eq.13.or.istk(il).eq.15) then
127          kopt=top2
128          isim=2
129       endif
130       top2=top2+1
131       il=iadr(lstk(top2))
132 c     
133 c Compute icontr
134 c - icontr=1 if without constraints
135 c - icontr=2 if with bound constraints
136 c
137 c     contraintes de borne :chaine "b" (code 11), xinf , xsup
138       if(istk(il).eq.10.and.istk(il+5).eq.2.and.
139      +     abs(istk(il+6)).eq.11) then
140          if (rhs.lt.5) then
141             call error(39)
142             return
143          endif
144          top2=top2+1
145          il=iadr(lstk(top2))
146          if(istk(il).gt.2.or.istk(il).eq.0)  then
147             err=top2-topin
148             call error(54)
149             return
150          endif
151          nbi=istk(il+1)*istk(il+2)
152          if(istk(il).eq.1) then
153             nbi=nbi*(istk(il+3)+1)
154             lbi=sadr(il+4)
155          else
156             lbi=sadr(il+9+nbi)
157             nbi=(istk(il+8+nbi)-1)*(istk(il+3)+1)
158          endif
159          top2=top2+1
160          il=iadr(lstk(top2))
161          if(istk(il).gt.2.or.istk(il).eq.0)  then
162             err=top2-topin
163             call error(54)
164             return
165          endif
166          nbs=istk(il+1)*istk(il+2)
167          if(istk(il).eq.1) then
168             nbs=nbs*(istk(il+3)+1)
169             lbs=sadr(il+4)
170          else
171             lbs=sadr(il+9+nbs)
172             nbs=(istk(il+8+nbs)-1)*(istk(il+3)+1)
173          endif
174          if((nbs.ne.nbi)) then
175             call error(139)
176             return
177          endif
178          icontr=2
179          top2=top2+1
180          il=iadr(lstk(top2))
181       end if
182 c     
183 c     point initial
184       if(istk(il).gt.2.or.istk(il).eq.0)  then
185          err=top2-topin
186          call error(54)
187          return
188       endif
189       topx=top2
190       nx1=istk(il+1)
191       nx2=istk(il+2)
192       itvx=istk(il)
193       ilx=il
194       if(istk(il).eq.1) then
195          nx=nx1*nx2*(istk(il+3)+1)
196          lx=sadr(il+4)
197       else
198          nx=(istk(il+8+nx1*nx2)-1)*(istk(il+3)+1)
199          lx=sadr(il+9+nx1*nx2)
200       endif
201       if (icontr.ne.1.and.nx.ne.nbi) then
202          call error(135)
203          return
204       endif
205 c     a quoi servent les 2 lignes suivantes. elle pose pb pour le nom de la macro
206 c     simulateur dans les messages d'erreur
207 c     idstk(1,top-1)=nx1
208 c     idstk(2,top-1)=nx2
209 c     
210 c     stockage de g
211       lg=ldisp
212       ldisp=lg + nx
213       err=ldisp - lstk(bot)
214       if (err.gt.0) then
215          call error(17)
216          return
217       endif
218       if (top2.eq.top) go to 200
219       top2=top2+1
220       il=iadr(lstk(top2))
221 c     
222 c     choix d algorithme
223 c
224 c Compute ialg1
225 c     (26,23) > "qn" > Quasi-Newton > ialg1=1
226 c     (16,12) > "gc" > Gradient Conjugate > ialg1=2
227 c     (23,13) > "nd" > Non Differentiable > ialg1=10
228 c
229       if(istk(il).eq.10) then
230          if (istk(il+5)-1.ne.2) then
231             err=top2-topin
232             call error(36)
233             return
234          endif
235          ic1=abs(istk(il+6))
236          ic2=abs(istk(il+7))
237          ialg1=0
238          if (ic1.eq.26.and.ic2.eq.23) ialg1=1
239          if (ic1.eq.16.and.ic2.eq.12) ialg1=2
240          if (ic1.eq.23.and.ic2.eq.13) ialg1=10
241          if (ialg1.ne.0) then
242             ialg=ialg1
243             if (top2.eq.top) go to 200
244             top2=top2+1
245             il=iadr(lstk(top2))
246          end if
247       endif
248 c     
249 c     df0
250       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
251          df0=stk(sadr(il+4))
252          if (df0.le.0) then
253             call error(143)
254             return
255          endif
256          if (top.eq.top2) go to 200
257          top2=top2 + 1
258          il=iadr(lstk(top2))
259       endif
260 c     
261 c     mmx
262       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
263          l=sadr(il+4)
264          mmx=int(stk(l))
265          immx=1
266          if (top2.eq.top) go to 200
267          top2=top2+1
268          il=iadr(lstk(top2))
269       end if
270 c     
271 c     hot start (optimiseurs n1qn1 et qnbd)
272       if(istk(il).eq.1.and.istk(il+3).eq.0) then
273          if (ialg.ne.1) then
274             call error(137)
275             return
276          endif
277          ntv=istk(il+1)*istk(il+2)
278          if(icontr.eq.1.and.ntv.ne.nx*(nx+13)/2) then
279             err=top2-topin
280             call error(142)
281             return
282          endif
283          ltv=sadr(il+4)
284          indtv=1
285          if (top2.eq.top) go to 200
286          top2=top2+1
287          il=iadr(lstk(top2))
288       end if
289 c     
290 c     chaine 'ar'
291       if(istk(il).ne.10)  then
292          err=top2-topin
293          call error(55)
294          return
295       endif
296       if (istk(il+1)*istk(il+2).ne.1) then
297          err=top2-topin
298          call error(89)
299          return
300       endif
301       if(istk(il+5)-1.ne.2 ) then
302          err=top2-topin
303          call error(36)
304          return
305       endif
306       if(abs(istk(il+6))+256*abs(istk(il+7)).ne.coar) goto 150
307       if (top2.eq.top) go to 200
308       top2=top2+1
309       il=iadr(lstk(top2))
310 c     
311 c     napm et itmax
312       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
313          l=sadr(il+4)
314          napm=int(stk(l))
315          if (top2.eq.top) go to 200
316          top2=top2+1
317          il=iadr(lstk(top2))
318       end if
319       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
320          l=sadr(il+4)
321          itmax=int(stk(l))
322          if (top2.eq.top) go to 200
323          top2=top2+1
324          il=iadr(lstk(top2))
325       end if
326 c     
327 c     epsg,epsf,epsx (note : epsx est un vecteur)
328       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
329          epsg=stk(sadr(il+4))
330          if (top2.eq.top) go to 200
331          top2=top2+1
332          il=iadr(lstk(top2))
333       endif
334       if(istk(il).eq.1.and.istk(il+1)*istk(il+2).eq.1) then
335          epsf=stk(sadr(il+4))
336          if (top2.eq.top) go to 200
337          top2=top2+1
338          il=iadr(lstk(top2))
339       endif
340       if(istk(il).eq.1.and.istk(il+3).eq.0) then
341          if(istk(il+1)*istk(il+2).ne.nx) then
342             call error(138)
343             return
344          endif
345          iepsx=0
346          lepsx=sadr(il+4)
347          if (top2.eq.top) go to 200
348          top2=top2+1
349          il=iadr(lstk(top2))
350       endif
351 c     
352 c     chaine 'in'
353 c
354  150  if(istk(il).ne.10) then
355          err=top2-topin
356          call error(55)
357          return
358       endif
359       if (istk(il+1)*istk(il+2).ne.1) then
360          err=top2-topin
361          call error(89)
362          return
363       endif
364       if (istk(il+5)-1.ne.2) then
365          err=top2-topin
366          call error(36)
367          return
368       endif
369       if(abs(istk(il+6))+256*abs(istk(il+7)).eq. coin) then
370          if(isim.ne.1) then
371             buf='''in'' not allowed with simulator defined '//
372      $           'by a function'
373             call error(9999)
374             return
375          endif
376
377 c     on initialise nizs,nrzs,ndzs
378          indsim=10
379          if(isim.eq.1) then
380             call foptim(indsim,nx,stk(lx),stk(lf),stk(lg),
381      &           istk(lizs),sstk(lrzs),stk(ldzs))
382          else
383             call boptim(indsim,nx,stk(lx),stk(lf),stk(lg),
384      &           izs,rzs,dzs)
385          endif
386          if(err.gt.0.or.err1.gt.0) return
387 c     
388          if(indsim.le.0) then
389             indopt=-7
390             if(indsim.eq.0) indopt=0
391             go to 350
392          endif
393 c     stockage de izs,rzs,dzs dans la pile
394          l1=ldisp
395          lizs=iadr(l1)
396          lrzs=lizs+nizs
397          ldzs=sadr(lrzs+nrzs)
398          ldisp=ldzs + ndzs
399          err=ldisp - lstk(bot)
400          if (err.gt.0) then
401             call error(17)
402             return
403          endif
404          indsim=11
405          lstk(top+1)=ldisp
406          if(isim.eq.1) then
407             call  foptim(indsim,nx,stk(lx),stk(lf),stk(lg),
408      &           istk(lizs),sstk(lrzs),stk(ldzs))
409          endif
410          if(indsim.le.0) then
411             if(indsim.eq.0) indopt=0
412             go to 350
413          endif
414          if (top2.eq.top) go to 200
415          top2=top2 + 1
416          il=iadr(lstk(top2))
417       endif
418 c     
419 c     izs et dzs en entree (chaine 'ti' et/ou 'td' suivie du tableau)
420       if(istk(il).ne.10) then
421          err=top2-topin
422          call error(55)
423          return
424       endif
425       if (istk(il+1)*istk(il+2).ne.1) then
426          err=top2-topin
427          call error(89)
428          return
429       endif
430       if(istk(il+5)-1.ne.2) then
431          err=top2-topin
432          call error(36)
433          return
434       endif
435       if(abs(istk(il+6))+256*abs(istk(il+7)).eq.coti) then
436          if (top2.eq.top) go to 200
437          top2=top2+1
438          il=iadr(lstk(top2))
439          if(istk(il).eq.1.and.istk(il+3).eq.0) then
440             nizs=istk(il+1)*istk(il+2)
441             lizs1=sadr(il+4)
442             lizs=iadr(lizs1)
443             do 185 i=0,nizs-1
444                istk(lizs+i)=int(stk(lizs1+i))
445  185        continue
446             if (top2.eq.top) go to 200
447             top2=top2+1
448             il=iadr(lstk(top2))
449          endif
450       endif
451       if(istk(il).ne.10) then
452          err=top2-topin
453          call error(55)
454          return
455       endif
456       if (istk(il+1)*istk(il+2).ne.1) then
457          err=top2-topin
458          call error(89)
459          return
460       endif
461       if(istk(il+5)-1.ne.2) then
462          err=top2-topin
463          call error(36)
464          return
465       endif
466       if(abs(istk(il+6))+256*abs(istk(il+7)).eq.cotd)  then
467          if (top2.eq.top) go to 200
468          top2=top2+1
469          il=iadr(lstk(top2))
470          if(istk(il).eq.1.and.istk(il+3).eq.0) then
471             ndzs=istk(il+1)*istk(il+2)
472             ldzs=sadr(il+4)
473             if (top2.eq.top) go to 200
474             top2=top2+1
475             il=iadr(lstk(top2))
476          endif
477       endif
478 c     
479 c     mettre ti = izs et/ou td = dzs en sortie (chaine 'si' ou 'sd')
480       if(istk(il).ne.10) then
481          err=top2-topin
482          call error(55)
483          return
484       endif
485       if (istk(il+1)*istk(il+2).ne.1) then
486          err=top2-topin
487          call error(89)
488          return
489       endif
490       if(istk(il+5)-1.ne.2) then
491          err=top2-topin
492          call error(36)
493          return
494       endif
495       if(abs(istk(il+6))+256*abs(istk(il+7)).eq.cosi) then
496          ireci=1
497          if (top2.eq.top) go to 200
498          top2=top2+1
499          il=iadr(lstk(top2))
500       endif
501 c     
502       if(istk(il).ne.10) then
503          err=top2-topin
504          call error(55)
505          return
506       endif
507       if (istk(il+1)*istk(il+2).ne.1) then
508          err=top2-topin
509          call error(89)
510          return
511       endif
512       if(istk(il+5)-1.ne.2) then
513          err=top2-topin
514          call error(36)
515          return
516       endif
517       if(abs(istk(il+6))+256*abs(istk(il+7)).ne.cosd) then
518          err=top2-topin
519          call error(36)
520          return
521       endif
522       irecd=1
523 c     
524 c     fin epluchage liste appel
525  200  if (top.ne.top2) then
526          call error(39)
527          return
528       endif
529 c     
530 c     creation des variables contenant le simulateur et ind
531       top=top+1
532       topind=top
533       lstk(top)=ldisp
534       il=iadr(ldisp)
535       istk(il)=1
536       istk(il+1)=1
537       istk(il+2)=1
538       istk(il+3)=0
539       ldisp=sadr(il+4)
540       ldisp=ldisp+1
541 c     
542       top=top+1
543       lstk(top)=ldisp
544       il=iadr(ldisp)
545       istk(il)=1
546       istk(il+1)=il+2
547       istk(il+2)=kopt
548       istk(il+3)=topx
549       istk(il+4)=topind
550       ldisp=sadr(il+5)
551 c     
552 c     initialisation eventuelle de f et g
553       iarret=0
554       if (napm.lt.2.or.itmax.lt.1) iarret=1
555 c     
556       if((icontr.eq.1.and.(ialg.eq.2.or.ialg.eq.10)).or.
557      &     (icontr.eq.2.and.ialg.eq.1.and.indtv.eq.1) .or.
558      &     (iarret.eq.1)    )   then
559          indsim=4
560          lstk(top+1)=ldisp
561          if(isim.eq.1) then
562             call foptim(indsim,nx,stk(lx),stk(lf),stk(lg),
563      &           istk(lizs),sstk(lrzs),stk(ldzs))
564          else
565             call boptim(indsim,nx,stk(lx),stk(lf),stk(lg),
566      &           istk(lizs),sstk(lrzs),stk(ldzs))
567          endif
568          if(err.gt.0.or.err1.gt.0)return
569          if(indsim.le.0) then
570             indopt=-7
571             if(indsim.eq.0) indopt=0
572             go to 350
573          endif
574          if (napm.lt.2.or.itmax.lt.1) go to 300
575       endif
576 c     
577 c     appel de l optimiseur
578 c     
579 c     optimiseur n1qn1 : Quasi-Newton without constraints
580       if(icontr.eq.1.and.ialg.eq.1) then
581          lvar=ldisp
582          mode=3
583          ntv=nx*(nx+13)/2
584          ldisp=lvar + nx
585          if(indtv.eq.0) then
586             mode=1
587             ltv=lvar  + nx
588             ldisp=ltv + ntv
589          endif
590          err=ldisp - lstk(bot)
591          if (err.gt.0) then
592             call error(17)
593             return
594          endif
595          do 50 i=0,nx-1
596             stk(lvar+i)=0.10d+0
597  50      continue
598          nitv=0
599          lstk(top + 1)=ldisp
600 c     
601 c     mise en memoire de parametres d entree pour l affectation de indop
602          itmax1=itmax
603          napm1=napm
604          epsg1=epsg
605 c     
606          if(isim.eq.1) then
607             call n1qn1(foptim,nx,stk(lx),stk(lf),stk(lg),
608      &           stk(lvar),epsg,mode,itmax,napm,imp,io,stk(ltv),
609      &           istk(lizs),sstk(lrzs),stk(ldzs))
610          else
611             call n1qn1(boptim,nx,stk(lx),stk(lf),stk(lg),
612      &           stk(lvar),epsg,mode,itmax,napm,imp,io,stk(ltv),
613      &           istk(lizs),sstk(lrzs),stk(ldzs))
614          endif
615          if(err.gt.0.or.err1.gt.0) return
616 c     affectation de indopt
617          epsg=sqrt(epsg)
618          indopt=9
619          if(itmax.ge.itmax1) indopt=5
620          if(napm.ge.napm1) indopt=4
621          if(epsg1.ge.epsg) indopt=1
622          go to 300
623       endif
624
625 c     algorithme n1qn3 : Gradient Conjugate without constraints
626       if(icontr.eq.1.and.ialg.eq.2) then
627 c     calcul de epsrel
628          zng=0.0d+0
629          do 230 i=0,nx-1
630             zng=zng + stk(lg+i)**2
631  230     continue
632          zng=sqrt(zng)
633          if (zng.gt.0.0d+0) epsg=epsg/zng
634 c     calcul du scalaire dxmin
635          dxmin=stk(leps)
636          if (iepsx.eq.0) then
637             dxmin=stk(lepsx)
638             if (nx.gt.1) then
639                do 235 i=1,nx-1
640                   dxmin=min(dxmin,stk(lepsx+i))
641  235           continue
642             endif
643          endif
644 c     tableaux de travail (mmx=nombre de mises a jour)
645          if (immx.eq.0) mmx=10
646          ntv=4*nx + mmx*(2*nx + 1)
647          ltv=ldisp
648          ldisp=ltv + ntv
649          err=ldisp - lstk(bot)
650          if (err.gt.0) then
651             call error(17)
652             return
653          endif
654          lstk(top+1)=ldisp
655 c     
656          if(isim.eq.1) then
657             indsim=4
658             call foptim(indsim,nx,stk(lx),stk(lf),stk(lg),
659      &           istk(lizs),sstk(lrzs),stk(ldzs))
660             call n1qn3(foptim,fuclid,ctonb,ctcab,nx,stk(lx),stk(lf),
661      $           stk(lg),dxmin,df0,epsg,imp,io,mode,itmax,napm,
662      &           stk(Ltv),Ntv,istk(lizs),sstk(lrzs),stk(ldzs) )
663          else
664             indsim=4
665             call boptim(indsim,nx,stk(lx),stk(lf),stk(lg),
666      &           istk(lizs),sstk(lrzs),stk(ldzs))
667             call n1qn3(boptim,fuclid,ctonb,ctcab,nx,stk(lx),stk(lf),
668      &           stk(lg),dxmin,df0,epsg,imp,io,mode,itmax,napm,
669      &           stk(ltv),ntv,istk(lizs),sstk(lrzs),stk(ldzs) )
670          endif
671          if (err.gt.0.or.err1.gt.0) return
672          indopt=9
673          if (mode.eq.0) indopt=0
674          if (mode.eq.1) indopt=1
675          if (mode.eq.2) indopt=-10
676          if (mode.eq.3) indopt=7
677
678          if (mode.eq.4) indopt=5
679          if (mode.eq.5) indopt=4
680          if (mode.eq.6) indopt=3
681          if (mode.eq.7) indopt=7
682          go to 300
683       endif
684 c     
685 c     optimiseur n1fc1 : non smooth without constraints
686       if(icontr.eq.1.and.ialg.eq.10) then
687          if (immx.eq.0) mmx=10
688          nitv=2*mmx + 2
689          itv1=5*nx + (nx+4)*mmx
690          itv2=(mmx+9)*mmx + 8
691          err=ldisp + iepsx*nx + nitv/2 +1 +itv1 +itv2 -lstk(bot)
692          if (err.gt.0) then
693             call error(17)
694             return
695          endif
696          if (iepsx.eq.1) then
697             lepsx=ldisp
698             do 115 i=1,nx
699                stk(lepsx+i - 1)=zero
700  115        continue
701             ldisp=lepsx+nx
702          endif
703          litv=iadr(ldisp)
704          ltv1=sadr(litv+nitv)
705          ltv2=ltv1  + itv1
706          ldisp=ltv2 + itv2
707          lstk(top+1)=ldisp
708          if (isim.eq.1) then
709             call n1fc1(foptim,fuclid,nx,stk(lx),stk(lf),stk(lg),
710      &           stk(lepsx),df0,epsf,zero,imp,io,mode,itmax,napm,mmx,
711      &           istk(litv),stk(ltv1),stk(ltv2),istk(lizs),sstk(lrzs),
712      $           stk(ldzs))
713          else
714             call n1fc1(boptim,fuclid,nx,stk(lx),stk(lf),stk(lg),
715      &           stk(lepsx),df0,epsf,zero,imp,io,mode,itmax,napm,mmx,
716      &           istk(litv),stk(ltv1),stk(ltv2),istk(lizs),sstk(lrzs),
717      $           stk(ldzs))
718          endif
719          if (err.gt.0.or.err1.gt.0) return
720 c     interpretation de la cause de retour
721          indopt=9
722          if (mode.eq.0)indopt=0
723          if (mode.eq.1)indopt=2
724          if (mode.eq.2)indopt=-10
725          if (mode.eq.4)indopt=5
726          if (mode.eq.5)indopt=4
727          if (mode.eq.6)indopt=3
728          go to 300
729       endif
730 c     
731 c     optimiseur qnbd : quasi-newton with bound constraints
732       if(icontr.eq.2.and.ialg.eq.1) then
733          if (iepsx.eq.1) then
734             err=ldisp +nx - lstk(bot)
735             if (err.gt.0) then
736                call error(17)
737                return
738             endif
739             lepsx=ldisp
740             ldisp=lepsx + nx
741             do 118 i=0,nx-1
742                stk(lepsx+i)=zero
743  118        continue
744          endif
745          ntv1=nx*(nx+1)/2 + 4*nx + 1
746          nitv= 2*nx
747          if (indtv.eq.0) then
748             ntv=ntv1
749             err= ldisp + ntv + nitv/2 +1 - lstk(bot)
750             if (err.gt.0) then
751                call error(17)
752                return
753             endif
754             ltv=ldisp
755             litv=iadr(ltv+ntv)
756             lstk(top+1)=sadr(litv+nitv)
757          else
758             if (ntv.ne.ntv1+nitv) then
759                err=top2-topin
760                call error(142)
761                return
762             endif
763             ntv=ntv1
764             litv1=ltv+ntv
765             litv=iadr(litv1)
766             do 117 i=0,nitv-1
767                istk(litv+i)=int(stk(litv1+i))
768  117        continue
769          endif
770          indopt=1 +indtv
771          if(isim.eq.1) then
772             call qnbd(indopt,foptim,nx,stk(lx),stk(lf),
773      &           stk(lg),imp,io,zero,napm,itmax,epsf,epsg,
774      &           stk(lepsx),df0,stk(lbi),stk(lbs),nfac,
775      &           stk(ltv),ntv,istk(litv),nitv,
776      &           istk(lizs),sstk(lrzs),stk(ldzs))
777          else
778             call qnbd(indopt,boptim,nx,stk(lx),stk(lf),
779      &           stk(lg),imp,io,zero,napm,itmax,epsf,epsg,
780      &           stk(lepsx),df0,stk(lbi),stk(lbs),nfac,
781      &           stk(ltv),ntv,istk(litv),nitv,
782      &           istk(lizs),sstk(lrzs),stk(ldzs))
783          endif
784          if(err.gt.0.or.err1.gt.0) return
785          go to 300
786       endif
787 c     
788 c     optimiseur gcbd : Gradient Conjugate with bound constraints
789       if(icontr.eq.2.and.ialg.eq.2) then
790          nt=2
791          if (immx.eq.1) nt= max(1,mmx/3)
792          ntv=nx*(5 + 3*nt) + 2*nt +1
793          nitv=nx + nt + 1
794          err= ldisp + iepsx*nx + ntv + nitv/2 - lstk(bot)
795          if (err.gt.0) then
796             call error(17)
797             return
798          endif
799          if (iepsx.eq.1) then
800             lepsx=ldisp
801             ltv=lepsx + nx
802             do 119 i=0,nx-1
803                stk(lepsx+i)=zero
804  119        continue
805          else
806             ltv=ldisp
807          endif
808          litv=iadr(ltv+ntv)
809          lstk(top+1)=sadr(litv+nitv)
810          indopt=1
811          if (indtes.ne.0) indopt=indtes
812          if(isim.eq.1) then
813             call gcbd(indopt,foptim,nomsub,nx,stk(lx),
814      &           stk(lf),stk(lg),imp,io,zero,napm,itmax,epsf,epsg,
815      &           stk(lepsx),df0,stk(lbi),stk(lbs),nfac,
816      &           stk(ltv),ntv,istk(litv),nitv,
817      &           istk(lizs),sstk(lrzs),stk(ldzs))
818          else
819             call gcbd(indopt,boptim,nomsub,nx,stk(lx),stk(lf),
820      &           stk(lg),imp,io,zero,napm,itmax,epsf,epsg,
821      &           stk(lepsx),df0,stk(lbi),stk(lbs),nfac,
822      &           stk(ltv),ntv,istk(litv),nitv,
823      &           istk(lizs),sstk(lrzs),stk(ldzs))
824          endif
825          if(err.gt.0.or.err1.gt.0) return
826          go to 300
827       endif
828 c     
829 c     algorithme non implante
830       call error(136)
831       return
832 c     
833 c     laissons la pile aussi propre qu on aurait aime la trouver
834  300  top2=top3
835       top =top3 + lhs - 1
836       fun=0
837 c     
838       lhs1=lhs - ireci -irecd
839       if (lhs1.le.0) then
840          call error(41)
841          return
842       endif
843 c     
844 c     sauvegarde de f
845       il=iadr(lstk(top2))
846       istk(il)=1
847       istk(il+1)=1
848       istk(il+2)=1
849       istk(il+3)=0
850       l=sadr(il+4)
851       stk(l)=stk(lf)
852       lstk(top+1)=l+1
853       if(lhs.eq.1) go to 320
854 c     
855 c     sauvegarde de x
856       l=l+1
857       top2=top2 + 1
858       lstk(top2)=l
859       il=iadr(l)
860       if(itvx.eq.1) then
861          call icopy(4,istk(ilx),1,istk(il),1)
862          lx1=sadr(il+4)
863       else
864          call icopy(9+nx1*nx2,istk(ilx),1,istk(il),1)
865          lx1=sadr(il+9+nx1*nx2)
866       endif
867       ilx=il
868       lstk(top+1)=lx1+nx
869       call unsfdcopy(nx,stk(lx) ,1,stk(lx1),1)
870       if(lhs1.eq.2) go to 320
871 c     
872 c     sauvegarde de g
873       top2=top2 + 1
874       lstk(top2)=lstk(top+1)
875       il=iadr(lstk(top2))
876       if(itvx.eq.1) then
877          call icopy(4,istk(ilx),1,istk(il),1)
878          l=sadr(il+4)
879       else
880          call icopy(9+nx1*nx2,istk(ilx),1,istk(il),1)
881          l=sadr(il+9+nx1*nx2)
882       endif
883       call unsfdcopy(nx,stk(lg) ,1,stk(l),1)
884       lstk(top+1)=l+nx
885       if(lhs1.eq.3) goto 320
886 c     
887 c     sauvegarde de tv (tableau interne a l' optimiseur - pour hot start
888       if (lhs1.eq.4) then
889          istv=0
890          if(ialg.eq.1) istv=1
891          if (istv.eq.0) then
892 c     pas de hot start pour cet algorithme
893             call error(137)
894             return
895          endif
896          top2=top2 + 1
897          lstk(top2)=lstk(top+1)
898          il=iadr(lstk(top2))
899          istk(il)=1
900          istk(il+1)=1
901          istk(il+2)=ntv + nitv
902          istk(il+3)=0
903          l=sadr(il+4)
904 c     recopie eventuelle de dzs et izs pour ne pas les ecraser
905          if (indtv.eq.0.and.(ireci-1)*(irecd-1).eq.0) then
906             err=l+ntv+nitv+ireci*nizs+irecd*ndzs-lstk(bot)
907             if (err.gt.0) then
908                call error(17)
909                return
910             endif
911 c     if (l+ntv+nitv+1.ge.sadr(lizs)) then
912             ldzs2=lstk(bot)-ndzs
913             lizs2=iadr(ldzs2)-nizs
914             if (sadr(lizs2).le.ltv+ntv+nitv+1) then
915                call error(17)
916                return
917             endif
918             call unsfdcopy(ndzs,stk(ldzs),-1,stk(ldzs2),-1)
919             do 315 i=nizs-1,0,-1
920                istk(lizs2+i)=istk(lizs+i)
921  315        continue
922             ldzs=ldzs2
923             lizs=lizs2
924          endif
925       endif
926       call unsfdcopy(ntv,stk(ltv),1,stk(l),1)
927       if (nitv.gt.0) then
928          do 316 i=nitv-1,0,-1
929             istk(litv+2*i)=istk(litv+i)
930  316     continue
931          litv1=l+ntv
932          do 317 i=0,nitv-1
933             stk(litv1+i)=real(istk(litv+i))
934  317     continue
935       endif
936       lstk(top+1)=l + ntv + nitv
937 c     
938 c     sauvegarde de izs et dzs
939  320  if (lhs.eq.lhs1) go to 350
940       if (ireci.eq.1) then
941          top2=top2 + 1
942          lstk(top2)=lstk(top+1)
943          il=iadr(lstk(top2))
944          istk(il)=1
945          istk(il+1)=1
946          istk(il+2)=nizs
947          istk(il+3)=0
948          l=sadr(il+4)
949          do 325 i=0,nizs-1
950             stk(l+i)=real(istk(lizs+i))
951  325     continue
952          lstk(top+1)=l+nizs
953       endif
954       if (irecd.eq.1) then
955          top2=top2 + 1
956          lstk(top2)=lstk(top+1)
957          il=iadr(lstk(top2))
958          istk(il)=1
959          istk(il+1)=1
960          istk(il+2)=ndzs
961          istk(il+3)=0
962          l=sadr(il+4)
963          call unsfdcopy(ndzs,stk(ldzs) ,1,stk(l),1)
964          lstk(top+1)=l+ndzs
965       endif
966       go to 350
967 c     
968 c     commentaires finaux
969  350  continue
970       if (iarret.eq.1) return
971       if (indopt.gt.0) go to 360
972       if(indopt.eq.0) then
973          call error(131)
974          return
975       elseif(indopt.eq.-7) then
976          call error(134)
977          return
978       elseif (indopt.eq.-14) then
979          call error(133)
980          return
981       elseif(indopt.le.-10) then
982          call error(132)
983          return
984       endif
985
986  360  continue
987       if(imp.ne.0) then
988          if(indopt.eq.1) then
989             call writebufscioptim(buf,epsg)
990             call msgs(12,0)
991          elseif(indopt.eq.2) then
992             call writebufscioptim(buf,epsg)
993             call msgs(13,0)
994          elseif(indopt.eq.3)  then
995             call msgs(14,0)
996          elseif(indopt.eq.4)  then
997             call msgs(15,0)
998          elseif(indopt.eq.5)  then
999             call msgs(16,0)
1000          elseif(indopt.eq.6)  then
1001             call msgs(17,0)
1002          elseif(indopt.eq.7)  then
1003             call msgs(18,0)
1004          elseif(indopt.eq.8)  then
1005             call msgs(19,0)
1006          elseif(indopt.eq.9)  then
1007             call msgs(20,0)
1008          elseif(indopt.ge.10)  then
1009             call msgs(21,0)
1010          endif
1011       endif
1012       return
1013       end
1014 c     --------------------------
1015