1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2c Copyright (C) 1988 - INRIA - Jean-Charles GILBERT
3c Copyright (C) 1988 - INRIA - Claude LEMARECHAL
4c
5c Copyright (C) 2012 - 2016 - Scilab Enterprises
6c
7c This file is hereby licensed under the terms of the GNU GPL v2.0,
8c pursuant to article 5.3.4 of the CeCILL v.2.1.
9c This file was originally licensed under the terms of the CeCILL v2.1,
10c and continues to be available under such terms.
11c For more information, see the COPYING file which you should have received
12c along with this program.
13c
14      subroutine n1qn2 (simul,prosca,n,x,f,g,dxmin,df1,epsg,iprint,io,
15     /                  mode,niter,nsim,dz,ndz,izs,rzs,dzs)
16c!But
17c     Minimisation sans contrainte par un algorithme
18c     de quasi-Newton a memoire limitee.
19c!Origine
20c     Version 1.0 de n1qn2 (Modulopt, INRIA), septembre 1988.
21c!Commentaires
22c     Ce code est en principe destine aux problemes de grande taille,
23c     n grand, mais convient egalement pour n quelconque. La methode
24c     utilisee est du type quasi-Newton (BFGS) a encombrement variable,
25c     ce qui permet d'utiliser au maximum la memoire declaree dispo-
26c     nible. On estine que plus la memoire utilisee est importante, plus
27c     rapide sera la decroissance du critere f.
28c!Sous-routines appelees
29c     n1qn2:    routine-chapeau qui structure la memoire declaree dispo-
30c               nible et appelle n1qn2a,
31c     n1qn2a:   optimiseur proprement dit,
32c     strang:   routine de calcul de la direction de descente,
33c     nlis0:    routine de recherche lineaire.
34c
35c     De son cote, l'utilisateur doit fournir:
36c     1) une routine qui appelle le module d'optimisation n1qn2,
37c     2) une routine de simulation, appelee simul par n1qn2, qui
38c        calcule la valeur de f et de son gradient en un point donne,
39c     3) une routine, appelee prosca par n1qn2, qui realise le produit
40c        scalaire de deux vecteurs, ce produit scalaire doit etre
41c        celui utilise pour calculer le gradient de f dans simul.
42c!Liste d'appel
43c     subroutine n1qn2 (simul,prosca,n,x,f,g,dxmin,df1,epsg,iprint,io,
44c    /                  mode,niter,nsim,dz,ndz,izs,rzs,dzs)
45c
46c     Dans la description des arguments qui suit, (e) signifie que
47c     l'argument doit etre initialise avant l'appel de n1qn2, (s)
48c     signifie que l'argument est une variable n'ayant de signification
49c     qu'en sortie et (es) = (e)+(s). Les arguments du type (s) et (es)
50c     sont en general modifies par n1qn2 et ne peuvent donc pas etre
51c     des constantes.
52c
53c     simul:     Nom d'appel de la sous-routine de simulation qui
54c                qui calcule la valeur de f et de son gradient g
55c                a l'itere courant. Ce module doit se presenter comme
56c                suit:
57c                   subroutine simul (indic,n,x,f,g,izs,rzs,dzs).
58c                Le nom de la sous-routine doit etre declare external
59c                dans le module appelant n1qn2. Les arguments n, x, f,
60c                g, izs, rzs et dzs ont la meme signification que ci-
61c                dessous. N1qn2 appelle simul soit avec indic=1, dans
62c                ce cas le simulateur fera ce qu'il veut mais ne chan-
63c                gera pas la valeur des arguments, ou avec indic=4, dans
64c                cas le simulateur calculera a la fois f et g.
65c     prosca:    Nom d'appel de la sous-routine effectuant le produit
66c                scalaire de deux vecteurs u et v. Ce module doit se
67c                presente sous la forme:
68c                   subroutine prosca (n,u,v,ps,izs,rzs,dzs).
69c                Le nom de la sous-routine doit etre declare external
70c                dans le module appelant n1qn2. Les argument n, izs,
71c                rzs et dzs ont la meme signification que ci-dessous.
72c                Les arguments u, v et ps sont des vecteurs de dimension
73c                n du type double precision. Ps donne le produit
74c                scalaire de u et v.
75c     n(e):      Scalaire du type integer. Donne la dimension n de la
76c                variable x.
77c     x(es):     Vecteur de dimension n du type double precision. En
78c                entree il faut fournir la valeur du point initial, en
79c                sortie, c'est le point final calcule par n1qn2.
80c     f(es):     Scalaire du type double precision. En entree, c'est la
81c                valeur de f en x (initial), valeur que l'on obtiendra
82c                en appelant le simulateur simul avant d'appeler n1qn2.
83c                En sortie, c'est la valeur de f en x (final).
84c     g(es):     Vecteur de dimension n du type double precision.
85c                En entree, il faut fournir la valeur du gradient de f
86c                en x (initial), valeur que l'on obtiendra en appelant
87c                le simulateur simul avant d'appeler n1qn2. En sortie en
88c                mode 1, c'est la valeur du gradient de f en x (final).
89c     dxmin(e):  Scalaire du type double precision, strictement positif.
90c                Cet argument definit la resolution sur x en norme
91c                l-infini: deux points dont la distance en norme l-
92c                infini est superieure a dxmin seront consideres comme
93c                non distinguables par la routine de recherche lineaire.
94c     df1(e):    Scalaire du type double precision, strictement positif.
95c                Cet argument donne une estimation de la diminution
96c                escomptee pour f lors de la premiere iteration.
97c     epsg(es):  Scalaire du type double precision, strictement positif
98c                et strictement inferieur a 1. Sa valeur en entree,
99c                determine le test d'arret qui porte sur la norme
100c                (prosca) du gradient. Le minimiseur considere que la
101c                convergence est atteinte en x(k) et s'arrete en mode 1
102c                si E(k) := |g(k)|/|g(1)| < epsg, ou g(1) et g(k) sont
103c                les gradients au point d'entree et a l'iteration k,
104c                respectivement. En sortie, epsg = E(k).
105c     iprint(e): Scalaire du type integer qui controle les sorties.
106c                <0:  Rien n'est imprime et n1qn2 appelle le simulateur
107c                     avec indic=1 toutes les (-iprint) iterations.
108c                =0:  Rien n'est imprime.
109c                >=1: Impressions initiales et finales, messages
110c                     d'erreurs.
111c                >=3: Une ligne d'impression par iteration donnant
112c                     l'ordre k de l'iteration courante menant de x(k)
113c                     a x(k+1), le nombre d'appels au simulateur avant
114c                     cette iteration, la valeur du critere f et sa
115c                     derivee directionnelle suivant d(k).
116c                >=4: Impression de nlis0.
117c                >=5: Impressions supplentaires en fin d'iteration k:
118c                     le test d'arret E(k+1), prosca(y(k),s(k)) qui
119c                     doit etre positif, le facteur de Oren-Spedicato
120c                     et l'angle de la direction de descente d(k) avec
121c                     -g(k).
122c     io(e):     Scalaire du type integer qui sera pris comme numero
123c                de canal de sortie pour les impressions controlees
124c                par iprint.
125c     mode(s):   Scalaire du type integer donnant le mode de sortie de
126c                n1qn2.
127c                <0: Impossibilite de poursuivre la recherche lineaire
128c                    car le simulateur a repondu avec indic<0. Mode
129c                    renvoie cette valeur de indic.
130c                =0: Arret demande par le simulateur qui a repondu avec
131c                    indic=0.
132c                =1: Fin normale avec test sur le gradient satisfait.
133c                =2: Arguments d'entree mal initialises. Il peut s'agir
134c                    de n<=0, niter<=0, nsim<=0, dxmin<=0, epsg<=0,
135c                    epsg>1 ou de nrz<5n+1 (pas assez de memoire
136c                    allouee).
137c                =3: La recherche lineaire a ete bloquee sur tmax=10^20
138c                    (mode tres peu probable).
139c                =4: Nombre maximal d'iterations atteint.
140c                =5: Nombre maximal de simulations atteint.
141c                =6: Arret sur dxmin lors de la recherche lineaire. Ce
142c                    mode de sortie peut avoir des origines tres
143c                    diverses. Si le nombre d'essais de pas lors de la
144c                    derniere recherche lineaire est faible, cela peut
145c                    signifier que dxmin a ete pris trop grand. Il peut
146c                    aussi s'agir d'erreurs ou d'imprecision dans le
147c                    calcul du gradient. Dans ce cas, la direction de
148c                    recherche d(k) peut ne plus etre une direction de
149c                    descente de f en x(k), etant donne que n1qn2
150c                    s'autorise des directions d(k) pouvant faire avec
151c                    -g(k) un angle proche de 90 degres. On veillera
152c                    donc a calculer le gradient avec precision.
153c                =7: Soit (g,d) soit (y,s) ne sont pas positifs (mode
154c                    de sortie tres improbable).
155c     niter(es): Scalaire du type integer, strictement positif. En
156c                entree, c'est le nombre maximal d'iterations admis.
157c                En sortie, c'est le nombre d'iterations effectuees.
158c     nsim(es):  Scalaire du type integer, strictement positif. En
159c                entree, c'est le nombre maximal de simulations admis.
160c                En sortie, c'est le nombre de simulations effectuees.
161c     rz(s):     Vecteur de dimension nrz du type double precision.
162c                C'est l'adresse d'une zone de travail pour n1qn2.
163c     nrz(e):    Scalaire du type integer, strictement positif, donnant
164c                la dimension de la zone de travail rz. En plus des
165c                vecteurs x et g donnes en arguments, n1qn2 a besoin
166c                d'une zone de travail composee d'au moins trois
167c                vecteurs de dimension n et chaque mise a jour demandee
168c                necessite 1 scalaire et 2 vecteurs de dimension n
169c                supplementaires. Donc si m est le nombre de mises a
170c                jour desire pour la construction de la metrique
171c                locale, il faudra prendre
172c                   nrz >= 3*n + m*(2*n+1).
173c                En fait, le nombre m est determine par n1qn2 qui prend:
174c                   m = partie entiere par defaut de ((nrz-3*n)/(2*n+1))
175c                Ce nombre doit etre >= 1. Il faut donc nrz >= 5*n +1,
176c                sinon n1qn2 s'arrete en mode 2.
177c     izs, rzs, dzs:  Adresses de zones-memoire respectivement du type
178c                integer, real et double precision. Elles sont reservees
179c                a l'utilisateur. N1qn2 ne les utilise pas et les
180c                transmet a simul et prosca.
181c!
182c-----------------------------------------------------------------------
183c
184c     arguments
185c
186      integer n,iprint,io,mode,niter,nsim,ndz,izs(*)
187      real rzs(*)
188      double precision x(*),f,g(*),dxmin,df1,epsg,dz(*)
189      double precision dzs(*)
190      external simul,prosca
191c
192c     variables locales
193c
194      integer m,ndzu,l1memo,id,igg,iaux,ialpha,iybar,isbar
195      double precision r1,r2
196      double precision ps
197      character bufstr*(4096)
198c
199c---- impressions initiales et controle des arguments
200c
201      if (iprint.ge.1) then
202         write (bufstr,900)
203         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
204
205         write (bufstr,9001) n
206         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
207
208         write (bufstr,9002) dxmin
209         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
210
211         write (bufstr,9003) df1
212         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
213
214         write (bufstr,9004) epsg
215         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
216
217         write (bufstr,9005) niter
218         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
219
220         write (bufstr,9006) nsim
221         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
222
223         write (bufstr,9006) iprint
224         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
225      endif
226900   format (' n1qn2: point d''entree')
2279001  format (5x,'dimension du probleme (n)              :',i6)
2289002  format (5x,'precision absolue en x (dxmin)         :',d9.2)
2299003  format (5x,'decroissance attendue pour f (df1)     :',d9.2)
2309004  format (5x,'precision relative en g (epsg)         :',d9.2)
2319005  format (5x,'nombre maximal d''iterations (niter)    :',i6)
2329006  format (5x,'nombre maximal d''appels a simul (nsim) :',i6)
2339007  format (5x,'niveau d''impression (iprint)           :',i4)
234      if (n.le.0.or.niter.le.0.or.nsim.le.0.or.dxmin.le.0.0d+0
235     /    .or.epsg.le.0.0d+0.or.epsg.gt.1.0d+0) then
236          mode=2
237          if (iprint.ge.1) then
238            write (bufstr,901)
239            call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
240            endif
241901       format (' >>> n1qn2 : appel incoherent')
242          goto 904
243      endif
244      if (ndz.lt.5*n+1) then
245          mode=2
246          if (iprint.ge.1) then
247            write (bufstr,902)
248            call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
249            endif
250902       format (' >>> n1qn2: memoire allouee insuffisante')
251          goto 904
252      endif
253c
254c---- calcul de m et des pointeurs subdivisant la zone de travail dz
255c
256      ndzu=ndz-3*n
257      l1memo=2*n+1
258      m=ndzu/l1memo
259      ndzu=m*l1memo+3*n
260      if (iprint.ge.1) then
261        write (bufstr,903) ndz
262        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
263        write (bufstr,9031) ndzu
264        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
265        write (bufstr,9032) m
266        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
267      endif
268903   format (5x,'memoire allouee (ndz)  :',i7)
2699031  format (5x,'memoire utilisee       :',i7)
2709032  format (5x,'nombre de mises a jour :',i6)
271      id=1
272      igg=id+n
273      iaux=igg+n
274      ialpha=iaux+n
275      iybar=ialpha+m
276      isbar=iybar+n*m
277c
278c---- appel du code d'optimisation
279c
280      call n1qn2a (simul,prosca,n,x,f,g,dxmin,df1,epsg,
281     /             iprint,io,mode,niter,nsim,m,
282     /             dz(id),dz(igg),dz(iaux),
283     /             dz(ialpha),dz(iybar),dz(isbar),izs,rzs,dzs)
284c
285c---- impressions finales
286c
287904   continue
288      if (iprint.ge.1) then
289        write (bufstr,905)
290        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
291        write (bufstr,9051) mode
292        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
293        write (bufstr,9052) niter
294        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
295        write (bufstr,9053) nsim
296        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
297        write (bufstr,9054) epsg
298        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
299      endif
300905   format (1x,79('-'))
3019051  format (1x,'n1qn2 : sortie en mode ',i2)
3029052  format (5x,'nombre d''iterations              : ',i4)
3039053  format (5x,'nombre d''appels a simul          : ',i6)
3049054  format (5x,'precision relative atteinte sur g: ',d9.2)
305      call prosca (n,x,x,ps,izs,rzs,dzs)
306      r1=sqrt(ps)
307      call prosca (n,g,g,ps,izs,rzs,dzs)
308      r2=sqrt(ps)
309      if (iprint.ge.1) then
310        write (bufstr,906) r1
311        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
312
313        write (bufstr,9061) f
314        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
315
316        write (bufstr,9062) r2
317        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
318
319        endif
320906   format (5x,'norme de x = ',d15.8)
3219061  format (5x,'f          = ',d15.8)
3229062  format (5x,'norme de g = ',d15.8)
323c
324      return
325      end
326