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