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