← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/optim/n1fc1a.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/optim/n1fc1a.f
 
1
      subroutine n1fc1a(simul,prosca,n,mode,xn,fn,g,df0,eps0,dx,imp,
 
2
     &                  zero,io,ntot,iter,nsim,memax,s,gd,x,sa,gg,al,
 
3
     &                  aps,anc,poids,q,jc,ic,r,a,e,rr,xga,y,w1,w2,izs,
 
4
     &                  rzs,dzs)
 
5
C ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
6
c     Copyright INRIA
 
7
C ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
8
C minimisation d'une fonction hemiderivable par une methode de faisceau.
 
9
C la direction est obtenue par projection de l'origine
 
10
C sur un polyedre genere par un ensemble de gradients deja calcules
 
11
C et la recherche lineaire donne un pas de descente ou un pas nul.
 
12
C l'algorithme minimise f a eps0 pres (si convexite)
 
13
C ou eps0 est une tolerance donnee par l'utilisateur.
 
14
C
 
15
C        mode
 
16
C                <=0 mode=indic de simul
 
17
C                1 fin normale
 
18
C                2 appel incoherent
 
19
C                3 reduire l'echelle des x
 
20
C                4 max iterations
 
21
C                5 max simulations
 
22
C                6 impossible d'aller au dela de dx
 
23
C                7 fprf2 mis en echec
 
24
C                8 on commence a boucler
 
25
C      imp
 
26
C                <0 indic=1 toutes les -imp iterations
 
27
C                0 pas d'impressions
 
28
C                1 impressions initiales et finales
 
29
C                2 impressions a chaque convergence
 
30
C                3 une impression par iteration
 
31
C                4 informations n1fc1 et nlis2
 
32
C                >4 debug
 
33
C                         5 tolerances diverses
 
34
C                         6 poids
 
35
C                         >6 fprf2
 
36
C         --------------------------------------------------
 
37
C fait appel aux subroutines suivantes:
 
38
C --------subroutine fprf2 (calcul de la direction)
 
39
C --------subroutines fremf2 et ffinf1 (esclaves de fprf2)
 
40
C --------subroutine frdf1 (reduction du faisceau)
 
41
C --------subroutine nlis2 (recherche lineaire)
 
42
C --------subroutine simul (module de simulation)
 
43
C --------subroutine prosca (produit de dualite donnant le gradient)
 
44
C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
45
C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
46
      implicit double precision (a-h,o-z)
 
47
      external simul, prosca
 
48
      dimension xn(n), g(n), izs(*), dzs(*), x(n), gd(n), gg(n)
 
49
      dimension s(n), sa(n), jc(*), ic(*), r(*), a(*), e(*), rr(*),
 
50
     &          xga(*), y(*), w1(*), w2(*)
 
51
      dimension q(*), al(memax), aps(memax), anc(memax), poids(memax)
 
52
      real rzs(*)
 
53
      dimension i5(1), d3(1), d4(1)
 
54
C
 
55
C         initialisations
 
56
C
 
57
      itmax = iter
 
58
      iter = 0
 
59
      itimp = 0
 
60
      napmax = nsim
 
61
      nsim = 1
 
62
      logic = 1
 
63
      logic2 = 0
 
64
      tmax = 1.d20
 
65
      eps = df0
 
66
      epsm = eps
 
67
      df = df0
 
68
      mode = 1
 
69
      ntot = 0
 
70
      iflag = 0
 
71
C
 
72
C          initialisation du faisceau
 
73
C          calcul du diametre de l'epure et du test d'arret
 
74
C
 
75
      aps(1) = 0.d0
 
76
      anc(1) = 0.d0
 
77
      poids(1) = 0.d0
 
78
      nta = 0
 
79
      kgrad = 1
 
80
      memax1 = memax + 1
 
81
      do 50 i = 1,n
 
82
 50   q(i) = -g(i)
 
83
      call prosca(n,g,g,ps,izs,rzs,dzs)
 
84
      if (ps .gt. 0.d0) goto 60
 
85
      mode = 2
 
86
      if (imp .ne. 0) call n1fc1o(io,3,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
87
      goto 900
 
88
 60   diam2 = 100. * df0 * df0 / ps
 
89
      eta2 = 1.d-2 * eps0 * eps0 / diam2
 
90
      ap = zero * df0 / diam2
 
91
      if (imp .gt. 2) call n1fc1o(io,4,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
92
C
 
93
C              boucle
 
94
C
 
95
 100  iter = iter + 1
 
96
      itimp = itimp + 1
 
97
      if (iter .lt. itmax) goto 110
 
98
      if (imp .gt. 0) call n1fc1o(io,5,iter,i2,i3,i4,i5,d1,d2,d3,d4)
 
99
      mode = 4
 
100
      goto 900
 
101
 110  ntot = ntot + 1
 
102
      if (logic .eq. 3) ro = ro * dsqrt(s2)
 
103
      if (itimp .ne. -imp) goto 200
 
104
      itimp = 0
 
105
      indic = 1
 
106
      call simul(indic,n,xn,f,g,izs,rzs,dzs)
 
107
C
 
108
C         calcul de la direction
 
109
C
 
110
 200  eps = dmin1(eps,epsm)
 
111
      eps = dmax1(eps,eps0)
 
112
      call fremf2(prosca,iflag,n,ntot,nta,memax1,q,poids,e,a,r,izs,rzs,
 
113
     &            dzs)
 
114
      call fprf2(iflag,ntot,nv,io,zero,s2,eps,al,imp,u,eta2,memax1,jc,
 
115
     &           ic,r,a,e,rr,xga,y,w1,w2)
 
116
C
 
117
C         fin anormale de fprf2
 
118
C
 
119
      if (iflag .eq. 0) goto 250
 
120
      if (imp .gt. 0) call n1fc1o(io,6,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
121
      mode = 7
 
122
      goto 900
 
123
 250  nta = ntot
 
124
      call ffinf1(n,nv,jc,xga,q,s)
 
125
      u = dmax1(u,0.d0)
 
126
      s2 = dmax1(s2,0.d0)
 
127
C
 
128
C          tests d'arret (nb. si nr g est interieur a g
 
129
C                                alors nr g est "nul")
 
130
C
 
131
      if (nv .lt. n+2) goto 260
 
132
      eta2 = dmax1(eta2,10.d0*s2)
 
133
      if (imp .ge. 2) call n1fc1o(io,7,i1,i2,i3,i4,i5,eta2,d2,d3,d4)
 
134
 260  if (s2 .gt. eta2) goto 300
 
135
C
 
136
C         calcul de la precision
 
137
      z = 0.d0
 
138
      do 270 k = 1,nv
 
139
        j = jc(k) - 1
 
140
        if (j .gt. 0) z = z + xga(k)*poids(j)
 
141
 270  continue
 
142
      epsm = dmin1(eps,z)
 
143
      if (imp.ge.2) call n1fc1o(io,8,iter,nsim,i3,i4,i5,fn,epsm,s2,d4)
 
144
      if (epsm .gt. eps0) goto 280
 
145
      mode = 1
 
146
      if (imp .gt. 0) call n1fc1o(io,9,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
147
      goto 900
 
148
C
 
149
C         diminution de epsilon
 
150
 280  epsm = dmax1(0.1d0*epsm,eps0)
 
151
      eps = epsm
 
152
      if (logic .eq. 3) tol = 0.01d0 * eps
 
153
      iflag = 2
 
154
      goto 200
 
155
C
 
156
C                 suite des iterations
 
157
C                    impressions
 
158
C
 
159
 300  if (imp .gt. 3) call n1fc1o(io,10,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
160
      if (imp .gt. 2) call n1fc1o(io,11,iter,nsim,nv,i4,i5,fn,eps,s2,u)
 
161
      if (imp .ge. 6) call n1fc1o(io,12,ntot,i2,i3,i4,i5,d1,d2,d3,poids)
 
162
C                test de non-pivotage
 
163
      if (logic .ne. 3) goto 350
 
164
      z = 0.d0
 
165
      do 310 i = 1,n
 
166
        z1 = s(i) - sa(i)
 
167
 310  z = z + z1*z1
 
168
      if (z .gt. 10.d0*zero*zero*s2) goto 350
 
169
      if (imp .gt. 0) call n1fc1o(io,13,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
170
      mode = 8
 
171
      goto 900
 
172
C
 
173
C                recherche lineaire
 
174
C
 
175
 350  iflag = 3
 
176
      s3 = s2 + u*eps
 
177
      if (logic .eq. 3) goto 365
 
178
      ro = 2. * df / s3
 
179
      tol = 0.01d0 * eps
 
180
      goto 370
 
181
 365  ro = ro / dsqrt(s2)
 
182
      tol = dmax1(0.6d0*tol,0.01d0*eps0)
 
183
 370  fa = fn
 
184
      alfa = 0.2d0
 
185
      beta = 0.1d0
 
186
      fpn = -s3
 
187
      if (memax .eq. 1) tol = 0.d0
 
188
C                 calcul de la resolution minimale, fonction de dx
 
189
      tmin = 0.d0
 
190
      do 372 i = 1,n
 
191
 372  tmin = dmax1(tmin,dabs(s(i)/dx))
 
192
      tmin = 1.d0 / tmin
 
193
      if (iter .eq. 1) roa = ro
 
194
      call nlis2(simul,prosca,n,xn,fn,fpn,ro,tmin,tmax,s,s2,g,gd,alfa,
 
195
     &           beta,imp,io,logic,nsim,napmax,x,tol,ap,tps,tnc,gg,izs,
 
196
     &           rzs,dzs)
 
197
      if (logic.eq.0 .or. logic.eq.2 .or. logic.eq.3) goto 380
 
198
C                 sortie par anomalie dans nlis2
 
199
      if (imp .le. 0) goto 375
 
200
      if (logic.eq.6 .or. logic.lt.0)
 
201
     &  call n1fc1o(io,14,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
202
      if (logic .eq. 4) call n1fc1o(io,15,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
203
      if (logic .eq. 5) call n1fc1o(io,16,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
204
      if (logic .eq. 1) call n1fc1o(io,17,i1,i2,i3,i4,i5,d1,d2,d3,d4)
 
205
 375  if (logic .eq. 1) mode = 3
 
206
      if (logic .eq. 4) mode = 5
 
207
      if (logic .eq. 5) mode = 0
 
208
      if (logic .eq. 6) mode = 6
 
209
      if (logic .lt. 0) mode = logic
 
210
      goto 900
 
211
 380  if (logic .ne. 3) goto 385
 
212
      do 382 i = 1,n
 
213
 382  sa(i) = s(i)
 
214
 385  if (iter .gt. 1) goto 390
 
215
C
 
216
C              1ere iteration, ajustement de ap, diam et eta
 
217
      if (logic .eq. 0) tps = (fn-fa) - ro*fpn
 
218
      ap = zero * zero * dabs(tps) / (s2*ro*ro)
 
219
      ajust = ro / roa
 
220
      if (logic .ne. 3) diam2 = diam2 * ajust * ajust
 
221
      if (logic .ne. 3) eta2 = eta2 / (ajust*ajust)
 
222
      if (imp .ge. 2) call n1fc1o(io,18,i1,i2,i3,i4,i5,diam2,eta2,ap,d4)
 
223
 390  mm = memax - 1
 
224
      if (logic .eq. 2) mm = memax - 2
 
225
      if (ntot .le. mm) goto 400
 
226
C
 
227
C      reduction du faisceau pour entrer le nouvel element
 
228
C
 
229
      call frdf1(prosca,n,ntot,mm,kgrad,al,q,s,poids,aps,anc,memax1,r,e,
 
230
     &           ic,izs,rzs,dzs)
 
231
      iflag = 1
 
232
      nta = ntot
 
233
      if (imp .ge. 2)
 
234
     &  call n1fc1o(io,19,iter,nsim,ntot,i4,i5,fn,d2,d3,d4)
 
235
C
 
236
 400  if (imp .ge. 5) call n1fc1o(io,20,logic,i2,i3,i4,i5,ro,tps,tnc,d4)
 
237
      if (logic .eq. 3) goto 500
 
238
C
 
239
C                 iteration de descente
 
240
C
 
241
      iflag = min0(iflag,2)
 
242
      df = fa - fn
 
243
      if (ntot .eq. 0) goto 500
 
244
C
 
245
C               actualisation des poids
 
246
C
 
247
      s3n = ro * dsqrt(s2)
 
248
      do 430 k = 1,ntot
 
249
        nk = (k-1)*n + 1
 
250
        call prosca(n,q(nk),s,ps,izs,rzs,dzs)
 
251
        z1 = dabs(aps(k)+(-df+ro*ps))
 
252
        z2 = anc(k) + s3n
 
253
        poids(k) = dmax1(z1,ap*z2*z2)
 
254
        aps(k) = z1
 
255
        anc(k) = z2
 
256
 430  continue
 
257
C
 
258
C                actualisation de eps
 
259
C
 
260
      eps = ro * s3
 
261
      kgrad = ntot + 1
 
262
C
 
263
C       nouvel element du faisceau (pour les trois types de pas)
 
264
C
 
265
 500  nt1 = ntot + 1
 
266
      if (logic .eq. 3) goto 510
 
267
      aps(nt1) = 0.d0
 
268
      anc(nt1) = 0.d0
 
269
      poids(nt1) = 0.d0
 
270
      goto 520
 
271
 510  aps(nt1) = tps
 
272
      anc(nt1) = dsqrt(tnc)
 
273
      poids(nt1) = dmax1(tps,ap*tnc)
 
274
 520  nk = ntot * n
 
275
      do 530 i = 1,n
 
276
        nki = nk + i
 
277
 530  q(nki) = -g(i)
 
278
C
 
279
C      traitement pour logic=2 (on ajoute encore un gradient)
 
280
      if (logic .ne. 2) goto 550
 
281
      ntot = ntot + 1
 
282
      logic = 3
 
283
      logic2 = 1
 
284
      do 540 i = 1,n
 
285
 540  g(i) = gd(i)
 
286
      goto 390
 
287
 550  logic = logic - logic2
 
288
      logic2 = 0
 
289
      goto 100
 
290
C
 
291
C                epilogue
 
292
C
 
293
 900  if (iter .le. 1) goto 990
 
294
      do 910 i = 1,n
 
295
 910  g(i) = -s(i)
 
296
 990  return
 
297
      end