1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
subroutine qnbd(indqn,simul,n,x,f,g,imp,io,zero,
& napmax,itmax,epsf,epsg,epsx,df0,binf,bsup,nfac,
& trav,ntrav,itrav,nitrav,izs,rzs,dzs)
c!but
c code de minimisation d une fonction reguliere sous contraintes
c de bornes , aux normes modulopt
c!origine
c f. bonnans , inria , 1986
c!methode
c principe de l algorithme : quasi-newton + projection
c details dans le rapport inria n. 242,1983
c cette version permet de tester plusieurs variantes de l algorithme
c en modifiant certains parametres internes (cf comment dans le code)
c taille memoire necessaire de l ordre de n**2/2
c pour les problemes de grande taille le code gcbd est mieux adapte
c
c!sous programmes appeles
c zqnbd optimiseur effectif
c proj projection
c calmaj mise a jour du hessien
c ajour mise a jour des facteurs de choleski
c rlbd,satur recherche lineaire avec bornes
c
c!liste d'appel
c
c subroutine qnbd(indqn,simul,n,x,f,g,imp,io,zero,
c & napmax,itmax,epsf,epsg,epsx,df0,binf,bsup,nfac,
c & trav,ntrav,itrav,nitrav,izs,rzs,dzs)
c
c indqn indicateur de qnbd es
c en entree =1 standard
c =2 dh et indic initialises au debut de trav et itrav
c ifac,f,g initialises
c en sortie
c si < 0 incapacite de calculer un point meilleur que le point initial
c si = 0 arret demande par l utilisateur
c si > 0 on fournit un point meilleur que le point de depart
c < -10 parametres d entree non convenables
c = -6 arret lors du calcul de la direction de descente et iter=1
c = -5 arret lors du calcul de l approximation du hessien iter=1
c = -3 anomalie de simul : indic negatif en un point ou
c f et g ont ete precedemment calcules
c = -2 echec de la recherche lineaire a la premiere iteration
c = -1 f non definie au point initial
c = 1 arret sur epsg
c = 2 epsf
c = 3 epsx
c = 4 napmax
c = 5 itmax
c = 6 pente dans la direction opposee au gradient trop petite
c = 7 arret lors du calcul de la direction de descente
c = 8 arret lors du calcul de l approximation du hessien
c = 10 arret par echec de la recherche lineaire , cause non precisee
c = 11 idem avec indsim < 0
c = 12 un pas trop petit proche d un pas trop grand
c ceci peut resulter d une erreur dans le gradient
c = 13 trop grand nombre d appels dans une recherche lineaire
c simul voir les normes modulopt
c n dim de x e
c binf,bsup bornes inf,sup,de dim n e
c x variables a optimiser (controle) es
c f valeur du critere s
c g gradient de f s
c zero proche zero machine e
c napmax nombre maximum d appels de simul e
c itmax nombre maximum d iterations de descente e
c itrav vect travail dim nitrav=2n , se decompose en indic et izig
c nfac nombre de variables factorisees (e si indqn=2) s
c imp facteur d impression e
c varie de 0 (pas d impressions) a 3 (nombreuses impressions)
c io numero du fichier de resultats e
c epsx vect dim n precision sur x e
c epsf critere arret sur f e
c epsg arret si sup a norm2(g+)/n e
c trav vect travail dim ntrav
c il faut ntrav > n(n+1)/2 +6n
c df0>0 decroissance f prevue (prendre 1. par defaut) e
c izs,rzs,dzs : cf normes modulopt es
c!
c indications sur les variables internes a qnbd et zqnbd
c izig sert a la memorisation des contraintes (actif si izag>1)
c si i ne change pas d ens on enleve 1 a izig (positif)
c sinon on ajoute izag
c factorisation seulement si izig est nul
c dh estimation hessien dim n(n+1)/2 rangee en troismorceaux es
c indic(i) nouvel indice de l indice i
c indic vect dim n ordre de rangement des indices es
c pas necessaire de l initialiser si indqn=1
c
c parametres de la recherche lineaire
c amd,amf param. du test de wolfe . (.7,.1)
c napm nombre max d appels dans la rl (=15)
c
implicit double precision (a-h,o-z)
real rzs(*)
double precision dzs(*)
dimension binf(n),bsup(n),x(n),g(n),epsx(n)
dimension trav(ntrav),itrav(nitrav),izs(*)
external simul
c
if(imp.ge.1)write(io,1010)
1010 format(' *********** qnbd ****************')
c
c
c parametres caracteristiques de l algorithme
c si les parametres sont nuls l algorithme est celui du rr 242
c ig=1 test sur grad(i) pour relach var
c in=1 limite le nombre de factorisations par iter a n/10
c irel=1 test sur decroissance grad pour rel a iter courante
c epsrel taux de decroissance permettant le relachement (cf irit)
c iact blocage variables dans ib (gestion contraintes actives)
c ieps1 =1 eps1 egal a zero
c note eps1 correspond a eps(xk)
ig=0
in=0
irel=1
epsrel=.50d+0
izag=0
iact=1
ieps1=0
c
c decoupage du vecteur trav
n1=(n*(n+1))/2 +1
n2=n1+n
n3=n2+n
n4=n3+n
n5=n4+n-1
if(ntrav.lt.n5) then
if(imp.gt.0)write(io,110)ntrav,n5
110 format(' qnbd : ntrav=',i8,' devrait valoir ',i8)
indqn=-11
return
endif
ni1=n+1
if(nitrav.lt.2*n) then
ni2=2*n
if(imp.gt.0)write(io,111)nitrav,ni2
111 format(' qnbd : nitrav=',i8,'devrait valoir',i8)
indqn=-12
return
endif
call zqnbd(indqn,simul,trav(1),n,binf,bsup,x,f,g,zero,napmax,
&itmax,itrav,itrav(ni1),nfac,imp,io,epsx,epsf,epsg,trav(n1),
&trav(n2),trav(n3),trav(n4),df0,ig,in,irel,izag,iact,
&epsrel,ieps1,izs,rzs,dzs)
return
end
|
