/[MITgcm]/MITgcm_contrib/mlosch/optim_m1qn3/m1qn3_offline.F
ViewVC logotype

Contents of /MITgcm_contrib/mlosch/optim_m1qn3/m1qn3_offline.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.7 - (show annotations) (download)
Wed Mar 20 16:05:00 2019 UTC (14 months, 1 week ago) by mlosch
Branch: MAIN
CVS Tags: HEAD
Changes since 1.6: +2 -2 lines
add moderl to the list of stored variable

This fixes a problem with non-initialized data.
(thanks to Andrew MacRae: https://github.com/dorugeber/MITgcm/commit/2fd395f)

1 C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/m1qn3_offline.F,v 1.6 2012/06/04 13:07:40 mlosch Exp $
2 C $Name: $
3 C
4 C The original m1qn3 has been modified to work "offline", i.e.
5 C the simulator and the driver of m1qn3_offline are separate
6 C programs that are called alternatingly from a (shell-)script.
7 C This requires that the "state" of m1qn3 is saved before
8 C this program terminates. Communication with the routine
9 C writing and restoring the state of m1qn3 is achieved via
10 C 3 new common-blocks that are contained in 3 header files.
11 C
12 C Apr27, 2012, Martin.Losch@awi.de
13 C
14 subroutine m1qn3_offline
15 & (simul,prosca,ctonb,ctcab,n,x,f,g,dxmin,df1,
16 & epsg,normtype,impres,io,imode,omode,niter,nsim,
17 & iz,dz,ndz,reverse,indic,izs,rzs,dzs)
18 c
19 c-----------------------------------------------------------------------
20 c
21 c M1QN3, Version 3.3, October 2009
22 c
23 c M1qn3 has two running modes: the SID (Scalar Initial Scaling) mode
24 c and the DIS (Diagonal Initial Scaling) mode. Both do not require
25 c the same amount of storage, the same subroutines, ...
26 c In the description below, items that differ in the DIS mode with
27 c respect to the SIS mode are given in brakets.
28 c
29 c Use the following subroutines:
30 c M1QN3A
31 c DDD, DDDS
32 c NLIS0 + DCUBE (Dec 88)
33 c MUPDTS, DYSTBL.
34 c
35 c The following routines are proposed to the user in case the
36 c Euclidean scalar product is used:
37 c DUCLID, DTONBE, DTCABE.
38 c
39 c La sous-routine M1QN3 est une interface entre le programme
40 c appelant et la sous-routine M1QN3A, le minimiseur proprement dit.
41 c
42 c Le module PROSCA est sense realiser le produit scalaire de deux
43 c vecteurs de Rn; le module DTONB est sense realiser le changement
44 c de coordonnees correspondant au changement de bases: base
45 c euclidienne -> base orthonormale (pour le produit scalaire
46 c PROSCA); le module CTBAB fait la transformation inverse: base
47 c orthonormale -> base euclidienne.
48 c
49 c Iz is an integer working zone for M1QN3A, its dimension is 5.
50 c It is formed of 5 scalars that are set by the optimizer:
51 c - the dimension of the problem,
52 c - an identifier of the scaling mode,
53 c - the number of updates,
54 c - two pointers.
55 c
56 c Dz est la zone de travail pour M1QN3A, de dimension ndz.
57 c Elle est subdivisee en
58 c 3 [ou 4] vecteurs de dimension n: d,gg,[diag,]aux
59 c m scalaires: alpha
60 c m vecteurs de dimension n: ybar
61 c m vecteurs de dimension n: sbar
62 c
63 c m est alors le plus grand entier tel que
64 c m*(2*n+1)+3*n .le. ndz [m*(2*n+1)+4*n .le. ndz)]
65 c soit m := (ndz-3*n) / (2*n+1) [m := (ndz-4*n) / (2*n+1)].
66 c Il faut avoir m >= 1, donc ndz >= 5n+1 [ndz >= 6n+1].
67 c
68 c A chaque iteration la metrique est formee a partir d un multiple
69 c de l'identite [d'une matrice diagonale] D qui est mise a jour m
70 c fois par la formule de BFGS en utilisant les m couples {y,s} les
71 c plus recents.
72 c
73 c-----------------------------------------------------------------------
74 c
75 c Authors: Jean Charles Gilbert, Claude Lemarechal, INRIA.
76 c
77 c Copyright 2008, 2009, INRIA.
78 c
79 c M1QN3 is distributed under the terms of the GNU General Public
80 c License.
81 c
82 c-----------------------------------------------------------------------
83 c
84 c This program is free software: you can redistribute it and/or
85 c modify it under the terms of the GNU General Public License as
86 c published by the Free Software Foundation, either version 3 of
87 c the License, or (at your option) any later version.
88 c
89 c This program is distributed in the hope that it will be useful,
90 c but WITHOUT ANY WARRANTY; without even the implied warranty of
91 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
92 c General Public License for more details.
93 c
94 c You should have received a copy of the GNU General Public License
95 c along with this program. If not, see
96 c <http://www.gnu.org/licenses/>.
97 c
98 c-----------------------------------------------------------------------
99 c
100 implicit none
101 c
102 c arguments
103 c
104 character*3 normtype
105 integer n,impres,io,imode(3),omode,niter,nsim,iz(5),ndz,indic,
106 & reverse,izs(*)
107 real rzs(*)
108 double precision x(n),f,g(n),dxmin,df1,epsg,dz(ndz),dzs(*)
109 external simul,prosca,ctonb,ctcab
110 c
111 c --- local variables
112 c
113 c The variable 'reentry' is used to know where to jump when entering
114 c the subroutine in reverse commnunication (it is better than using
115 c 'reverse', which is known outside m1qn3 and could be changed by
116 c the user):
117 c = 0: no jump, start at the begining of the subroutine
118 c = 1: jump inside m1qn3a, where the simulator was called with
119 c indic=1
120 c = 2: jump inside m1qn3a/mlis3, where the simulator was called with
121 c indic=4
122 c The variable 'reentry' is set to 0 by the compiler and when m1qn3
123 c has completed its job (see below); since the value is saved, it
124 c has always the value 0 when entering the solver for the first
125 c time when solving a new problem, even if it is in the same
126 c program.
127 c
128 logical sscale
129 integer ntravu,m,mmemo
130 double precision gnorm,ps
131 #include "m1qn3_common.h"
132 CML logical inmemo,sscale
133 CML integer ntravu,id,igg,idiag,iaux,ialpha,iybar,isbar,m,mmemo,
134 CML & reentry
135 CML double precision gnorm,ps
136 CML save inmemo,sscale,ntravu,id,igg,idiag,iaux,ialpha,iybar,isbar,m,
137 CML & mmemo,reentry,gnorm,ps
138 c
139 c data
140 c
141 CML this needs to be done outside of m1qn3_offline
142 CML data reentry /0/
143 c
144 c --- function
145 c
146 double precision ddot,dnrmi
147 c
148 c --- stop if reverse < 0 (m1qn3 should not be called with reverse < 0)
149 c
150 if (reverse.lt.0) then
151 write (io,'(/a,a,i0,a/)')
152 & " >>> m1qn3 should not be called with a negative reverse",
153 & " (=", reverse, ")"
154 stop
155 endif
156 c
157 c --- possible jumps
158 c 9999: go directly in m1qn3a
159 c
160 if (reentry.gt.0) goto 9999
161 c
162 c---- license notice
163 c
164 if (impres.ge.5) then
165 write (io,934)
166 endif
167 934 format (1x,79("-")
168 & /1x,"M1QN3 Copyright (C) 2008, J. Ch. Gilbert, Cl. ",
169 & "Lemarechal."
170 & /1x,79("-")
171 & /1x,"This program comes with ABSOLUTELY NO WARRANTY. This is",
172 & " free software, and you"
173 & /1x,"are welcome to redistribute it under certain ",
174 & "conditions. See the file COPYING "
175 & /1x,"in the root directory of the M1QN3 distribution for ",
176 & "details."
177 & /1x,79("-"))
178 c
179 c---- impressions initiales et controle des arguments
180 c
181 if (impres.ge.1) then
182 write (io,900) n,dxmin,df1,epsg,normtype,niter,nsim,impres
183 if (reverse.gt.0) then
184 write (io,'(5x,a)') "reverse communication"
185 else
186 write (io,'(5x,a)') "direct communication"
187 endif
188 endif
189 900 format (/" M1QN3 (Version 3.3, October 2009): entry point"/
190 & 5x,"dimension of the problem (n):",i14/
191 & 5x,"absolute precision on x (dxmin):",9x,1pd9.2/
192 & 5x,"expected decrease for f (df1):",11x,1pd9.2/
193 & 5x,"relative precision on g (epsg):",10x,1pd9.2,
194 & " (",a3,"-norm)"/
195 & 5x,"maximal number of iterations (niter):",i6/
196 & 5x,"maximal number of simulations (nsim):",i6/
197 & 5x,"printing level (impres):",15x,i4)
198 c
199 c---- check the arguments
200 c
201 if (n.le.0) then
202 omode=2
203 if (reverse.gt.0) reverse = -1
204 if (impres.ge.1) write (io,'(/a)')
205 & " >>> m1qn3: n should be > 0"
206 return
207 endif
208 if (niter.le.0) then
209 omode=2
210 if (reverse.gt.0) reverse = -1
211 if (impres.ge.1) write (io,'(/a)')
212 & " >>> m1qn3: niter should be > 0"
213 return
214 endif
215 if (nsim.le.0) then
216 omode=2
217 if (reverse.gt.0) reverse = -1
218 if (impres.ge.1) write (io,'(/a)')
219 & " >>> m1qn3: nsim should be > 0"
220 return
221 endif
222 if (dxmin.le.0.d0) then
223 omode=2
224 if (reverse.gt.0) reverse = -1
225 if (impres.ge.1) write (io,'(/a)')
226 & " >>> m1qn3: dxmin should be > 0.d0"
227 return
228 endif
229 c if (epsg.le.0.d0 .or. epsg.gt.1.d0) then
230 if (epsg.le.0.d0) then
231 omode=2
232 if (reverse.gt.0) reverse = -1
233 if (impres.ge.1) write (io,'(/a)')
234 & " >>> m1qn3: epsg should be > 0.d0"
235 return
236 endif
237 if (epsg.ge.1.d0) then
238 omode=1
239 niter=0
240 nsim=0
241 epsg=1.d0
242 if (reverse.gt.0) reverse = -1
243 if (impres.ge.1) write (io,'(/a)')
244 & " >>> m1qn3: epsg is >= 1.d0, no need to make progress"
245 goto 1000
246 endif
247 if ((normtype.ne.'two') .and.
248 & (normtype.ne.'sup') .and.
249 & (normtype.ne.'dfn')) then
250 omode=2
251 if (reverse.gt.0) reverse = -1
252 write (io,'(/a,a,a/)') " >>> m1qn3: unknown norm type '",
253 & normtype, "'"
254 return
255 endif
256 if (impres.lt.0) then
257 omode=2
258 if (reverse.gt.0) reverse = -1
259 write (io,'(/a,i0/)')
260 & " >>> m1qn3: impres should be >= 0 and has the value ",
261 & impres
262 return
263 endif
264 c
265 c---- what method
266 c
267 if (imode(1).eq.0) then
268 if (impres.ge.1) write (io,920)
269 920 format (/" m1qn3: Diagonal Initial Scaling mode")
270 sscale=.false.
271 else
272 if (impres.ge.1) write (io,921)
273 921 format (/" m1qn3: Scalar Initial Scaling mode")
274 sscale=.true.
275 endif
276 c
277 if ((ndz.lt.5*n+1).or.((.not.sscale).and.(ndz.lt.6*n+1))) then
278 omode=2
279 if (reverse.gt.0) reverse = -1
280 if (impres.ge.1) write (io,922)
281 922 format (/" >>> m1qn3: not enough memory allocated")
282 return
283 endif
284 c
285 c---- Compute m
286 c
287 call mupdts (sscale,inmemo,n,m,ndz)
288 c
289 c --- Check the value of m (if (y,s) pairs in core, m will be >= 1)
290 c
291 if (m.lt.1) then
292 omode=2
293 if (reverse.gt.0) reverse = -1
294 if (impres.ge.1) write (io,930)
295 930 format (/" >>> m1qn3: m is set too small in mupdts")
296 return
297 endif
298 c
299 c --- mmemo = number of (y,s) pairs in core memory
300 c
301 mmemo=1
302 if (inmemo) mmemo=m
303 c
304 ntravu=2*(2+mmemo)*n+m
305 if (sscale) ntravu=ntravu-n
306 if (impres.ge.1) write (io,931) ndz,ntravu,m
307 931 format (/5x,"allocated memory (ndz) :",i9/
308 & 5x,"used memory : ",i9/
309 & 5x,"number of updates : ",i9)
310 if (ndz.lt.ntravu) then
311 omode=2
312 if (reverse.gt.0) reverse = -1
313 if (impres.ge.1) write (io,922)
314 return
315 endif
316 c
317 if (impres.ge.1) then
318 if (inmemo) then
319 write (io,932)
320 else
321 write (io,933)
322 endif
323 endif
324 932 format (5x,"(y,s) pairs are stored in core memory")
325 933 format (5x,"(y,s) pairs are stored by the user")
326 c
327 c---- cold start or warm restart ?
328 c check iz: iz(1)=n, iz(2)=(0 if DIS, 1 if SIS),
329 c iz(3)=m, iz(4)=jmin, iz(5)=jmax
330 c
331 if (imode(2).eq.0) then
332 if (impres.ge.1) write (io,940)
333 else
334 if (iz(1).ne.n .or. iz(2).ne.imode(1) .or. iz(3).ne.m .or.
335 & iz(4).lt.1 .or. iz(5).lt.0 .or. iz(4).gt.iz(3) .or.
336 & iz(5).gt.iz(3))
337 & then
338 omode=2
339 if (reverse.gt.0) reverse = -1
340 if (impres.ge.1) then
341 write (io,941)
342 if (iz(1).ne.n) write (io,942)
343 if (iz(2).ne.imode(1)) write (io,943)
344 if (iz(3).ne.m) write (io,944)
345 if (iz(4).lt.1 .or. iz(5).lt.0 .or. iz(4).gt.iz(3)
346 & .or. iz(5).gt.iz(3)) write (io,945)
347 endif
348 return
349 endif
350 if (impres.ge.1) write (io,946)
351 endif
352 940 format (/" m1qn3: cold start"/1x)
353 941 format (/" >>> m1qn3: inconsistent warm restart ")
354 942 format (" >>> m1qn3: (the number of variables has changed)")
355 943 format (" >>> m1qn3: (the scaling mode has changed)")
356 944 format (" >>> m1qn3: (the number of updates has changed)")
357 945 format (" >>> m1qn3: (wrong pointers)")
358 946 format (/" m1qn3: warm restart"/1x)
359 iz(1)=n
360 iz(2)=0
361 if (sscale) iz(2)=1
362 iz(3)=m
363 c
364 c---- split the working zone dz
365 c
366 idiag=1
367 iybar=idiag+n
368 if (sscale) iybar=1
369 isbar=iybar+n*mmemo
370 id=isbar+n*mmemo
371 igg=id+n
372 iaux=igg+n
373 ialpha=iaux+n
374 c
375 c---- call the optimization code
376 c
377 9999 continue
378 call m1qn3a (simul,prosca,ctonb,ctcab,n,x,f,g,dxmin,df1,epsg,
379 & normtype,impres,io,imode,omode,niter,nsim,inmemo,
380 & iz(3),iz(4),iz(5),dz(id),dz(igg),dz(idiag),dz(iaux),
381 & dz(ialpha),dz(iybar),dz(isbar),reverse,reentry,indic,
382 & izs,rzs,dzs)
383 if (reentry.gt.0) return
384 c
385 c---- impressions finales
386 c
387 1000 continue
388 if (impres.ge.1) write (io,960) omode,niter,nsim,epsg
389 960 format (/1x,79("-")/
390 & /" m1qn3: output mode is ",i2
391 & /5x,"number of iterations: ",i14
392 & /5x,"number of simulations: ",i13
393 & /5x,"realized relative precision on g: ",1pd9.2)
394 if (normtype.eq.'two') then
395 gnorm = sqrt(ddot(n,g,1,g,1))
396 elseif (normtype.eq.'sup') then
397 gnorm = dnrmi(n,g)
398 elseif (normtype.eq.'dfn') then
399 call prosca (n,g,g,ps,izs,rzs,dzs)
400 gnorm=dsqrt(ps)
401 endif
402
403 if (impres.ge.1) write (io,961) f,normtype,gnorm
404 961 format (5x,"f = ",1pd15.8
405 & /5x,a3,"-norm of g = ",1pd15.8)
406
407 return
408 end
409 c
410 c--------0---------0---------0---------0---------0---------0---------0--
411 c
412 subroutine m1qn3a (simul,prosca,ctonb,ctcab,n,x,f,g,dxmin,df1,
413 & epsg,normtype,impres,io,imode,omode,niter,nsim,
414 & inmemo,m,jmin,jmax,d,gg,diag,aux,alpha,ybar,
415 & sbar,reverse,reentry,indic,izs,rzs,dzs)
416 c----
417 c
418 c Code d optimisation proprement dit.
419 c
420 c----
421 c
422 implicit none
423 c
424 c arguments
425 c
426 logical inmemo
427 character*3 normtype
428 integer n,impres,io,imode(3),omode,niter,nsim,m,jmin,jmax,indic,
429 & reverse,reentry,izs(*)
430 real rzs(*)
431 double precision x(n),f,g(n),dxmin,df1,epsg,d(n),gg(n),diag(n),
432 & aux(n),alpha(m),ybar(n,*),sbar(n,*),dzs(*)
433 CML & aux(n),alpha(m),ybar(n,1),sbar(n,1),dzs(*)
434 external simul,prosca,ctonb,ctcab
435 c
436 c variables locales
437 c
438 logical skip_update
439 integer i
440 #include "m1qn3a_common.h"
441 CML logical sscale,cold,warm,skip_update
442 CML integer i,itmax,moderl,isim,jcour
443 CML double precision d1,t,tmin,tmax,gnorm,gnorms,eps1,ff,preco,precos,
444 CML & ys,den,dk,dk1,ps,ps2,hp0
445 CML save sscale,cold,warm,i,itmax,moderl,isim,jcour,d1,t,tmin,tmax,
446 CML & gnorm,gnorms,eps1,ff,preco,precos,ys,den,dk,dk1,ps,ps2,hp0
447 c
448 c parametres
449 c
450 double precision rm1,rm2
451 parameter (rm1=0.0001d+0,rm2=0.99d+0)
452 double precision pi
453 parameter (pi=3.1415927d+0)
454 double precision rmin
455 c
456 c function
457 c
458 double precision ddot,dnrmi
459 c
460 c --- possible jumps
461 c 9998: call of the simulator in m1qn3a with indic = 1
462 c 9999: call of the simulator in mlis3 with indic = 4
463 c
464 if (reentry.eq.1) goto 9998
465 if (reentry.eq.2) goto 9999
466 c
467 c---- initialisation
468 c
469 rmin=1.d-20
470 c
471 sscale=.true.
472 if (imode(1).eq.0) sscale=.false.
473 c
474 warm=.false.
475 if (imode(2).eq.1) warm=.true.
476 cold=.not.warm
477 c
478 skip_update = .false.
479 c
480 itmax=niter
481 niter=0
482 isim=1
483 eps1=1.d+0
484 c
485 call prosca (n,g,g,ps,izs,rzs,dzs)
486 gnorm = dsqrt(ps)
487 if (normtype.eq.'two') then
488 gnorms = sqrt(ddot(n,g,1,g,1))
489 elseif (normtype.eq.'sup') then
490 gnorms = dnrmi(n,g)
491 elseif (normtype.eq.'dfn') then
492 gnorms = gnorm
493 endif
494 if (impres.ge.1) write (io,900) f,normtype,gnorms
495 900 format (5x,"f = ",1pd15.8
496 & /5x,a3,"-norm of g = ",1pd15.8)
497 if (gnorms.lt.rmin) then
498 omode=2
499 if (impres.ge.1) write (io,901)
500 goto 1000
501 endif
502 901 format (/" >>> m1qn3a: initial gradient is too small")
503 c
504 c --- initialisation pour dd
505 c
506 if (cold) then
507 jmin=1
508 jmax=0
509 endif
510 jcour=1
511 if (inmemo) jcour=jmax
512 c
513 c --- mise a l'echelle de la premiere direction de descente
514 c
515 if (cold) then
516 c
517 c --- use Fletcher's scaling and initialize diag to 1.
518 c
519 precos=2.d+0*df1/gnorm**2
520 do 10 i=1,n
521 d(i)=-g(i)*precos
522 diag(i)=1.d+0
523 10 continue
524 if (impres.ge.5) write(io,902) precos
525 902 format (/" m1qn3a: descent direction -g: precon = ",d10.3)
526 else
527 c
528 c --- use the matrix stored in [diag and] the (y,s) pairs
529 c
530 if (sscale) then
531 call prosca (n,ybar(1,jcour),ybar(1,jcour),ps,izs,rzs,dzs)
532 precos=1.d+0/ps
533 endif
534 do 11 i=1,n
535 d(i)=-g(i)
536 11 continue
537 if (inmemo) then
538 call dd (prosca,ctonb,ctcab,n,sscale,m,d,aux,jmin,jmax,
539 & precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
540 else
541 call dds (prosca,ctonb,ctcab,n,sscale,m,d,aux,jmin,jmax,
542 & precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
543 endif
544 endif
545 c
546 if (impres.eq.3) write(io,903)
547 if (impres.eq.4) write(io,903)
548 903 format (/1x,79("-"))
549 904 format (1x)
550 c
551 c --- initialisation pour mlis3
552 c
553 tmax=1.d+20
554 call prosca (n,d,g,hp0,izs,rzs,dzs)
555 if (hp0.ge.0.d+0) then
556 omode=7
557 if (impres.ge.1) write (io,905) niter,hp0
558 goto 1000
559 endif
560 905 format (/" >>> m1qn3 (iteration ",i2,"): "
561 & /5x," the search direction d is not a ",
562 & "descent direction: (g,d) = ",d12.5)
563 c
564 c --- compute the angle (-g,d)
565 c
566 if (warm.and.impres.ge.5) then
567 call prosca (n,g,g,ps,izs,rzs,dzs)
568 ps=dsqrt(ps)
569 call prosca (n,d,d,ps2,izs,rzs,dzs)
570 ps2=dsqrt(ps2)
571 ps=hp0/ps/ps2
572 ps=dmin1(-ps,1.d+0)
573 ps=dacos(ps)
574 d1=ps*180.d+0/pi
575 write (io,906) sngl(d1)
576 endif
577 906 format (/" m1qn3: descent direction d: ",
578 & "angle(-g,d) = ",f5.1," degrees")
579 c
580 c---- Debut de l iteration. on cherche x(k+1) de la forme x(k) + t*d,
581 c avec t > 0. On connait d.
582 c
583 c Debut de la boucle: etiquette 100,
584 c Sortie de la boucle: goto 1000.
585 c
586 100 niter=niter+1
587 if (impres.ge.5) write(io,903)
588 if (impres.ge.4) write(io,904)
589 if (impres.ge.4) write (io,910) niter,isim,f,hp0
590 910 format (" m1qn3: iter ",i0,", simul ",i0,
591 & ", f=",1pd15.8,", h'(0)=",d12.5)
592 c
593 c --- free simulation if desired
594 c
595 if (imode(3).gt.0) then
596 if (mod(niter-1,imode(3)).eq.0) then
597 indic=1
598 if (reverse.gt.0) then
599 reentry = 1
600 return
601 else
602 call simul(indic,n,x,f,g,izs,rzs,dzs)
603 endif
604 endif
605 endif
606 9998 continue
607 c
608 c --- recherche lineaire et nouveau point x(k+1)
609 c
610 do 101 i=1,n
611 gg(i)=g(i)
612 101 continue
613 ff=f
614 if (impres.ge.5) write (io,911)
615 911 format (/" m1qn3: line search")
616 c
617 c --- calcul de tmin
618 c
619 tmin=0.d+0
620 do 200 i=1,n
621 tmin=dmax1(tmin,dabs(d(i)))
622 200 continue
623 tmin=dxmin/tmin
624 t=1.d+0
625 d1=hp0
626 c
627 9999 continue
628 call mlis3 (n,simul,prosca,x,f,d1,t,tmin,tmax,d,g,rm2,rm1,impres,
629 & io,moderl,isim,nsim,aux,reverse,reentry,indic,izs,rzs,
630 & dzs)
631 if (reentry.gt.0) return
632 c
633 c --- mlis3 renvoie les nouvelles valeurs de x, f et g
634 c
635 if (moderl.ne.0) then
636 if (moderl.lt.0) then
637 c
638 c --- calcul impossible
639 c t, g: ou les calculs sont impossibles
640 c x, f: ceux du t_gauche (donc f <= ff)
641 c
642 omode=moderl
643 goto 1000
644 elseif (moderl.eq.1) then
645 c
646 c --- descente bloquee sur tmax
647 c
648 skip_update = .true.
649 c omode=3
650 c if (impres.ge.1) write(io,912) niter
651 c 912 format (/" >>> m1qn3 (iteration ",i0,
652 c & "): line search blocked on tmax: "/
653 c & " >>> possible reasons: bad scaling,",
654 c & " unbounded problem")
655 elseif (moderl.eq.4) then
656 c
657 c --- nsim atteint
658 c x, f: ceux du t_gauche (donc f <= ff)
659 c
660 omode=5
661 goto 1000
662 elseif (moderl.eq.5) then
663 c
664 c --- arret demande par l utilisateur (indic = 0)
665 c x, f: ceux en sortie du simulateur
666 c
667 omode=0
668 goto 1000
669 elseif (moderl.eq.6) then
670 c
671 c --- arret sur dxmin ou appel incoherent
672 c x, f: ceux du t_gauche (donc f <= ff)
673 c
674 omode=6
675 goto 1000
676 endif
677 else
678 skip_update = .false.
679 endif
680 c
681 c NOTE: stopping tests are now done after having updated the matrix, so
682 c that update information can be stored in case of a later warm restart
683 c
684 c --- mise a jour de la matrice
685 c
686 if (skip_update) then
687 if (impres.ge.5) write(io,'(/a)')
688 & " m1qn3: matrix update is skipped"
689 elseif (m.gt.0) then
690 c
691 c --- mise a jour des pointeurs
692 c
693 jmax=jmax+1
694 if (jmax.gt.m) jmax=jmax-m
695 if ((cold.and.niter.gt.m).or.(warm.and.jmin.eq.jmax)) then
696 jmin=jmin+1
697 if (jmin.gt.m) jmin=jmin-m
698 endif
699 if (inmemo) jcour=jmax
700 c
701 c --- y, s et (y,s)
702 c
703 do 400 i=1,n
704 sbar(i,jcour)=t*d(i)
705 ybar(i,jcour)=g(i)-gg(i)
706 400 continue
707 if (impres.ge.5) then
708 call prosca (n,sbar(1,jcour),sbar(1,jcour),ps,izs,rzs,dzs)
709 dk1=dsqrt(ps)
710 if (niter.gt.1) write (io,930) dk1/dk
711 930 format (/" m1qn3: convergence rate, s(k)/s(k-1) = ",
712 & 1pd12.5)
713 dk=dk1
714 endif
715 call prosca (n,ybar(1,jcour),sbar(1,jcour),ys,izs,rzs,dzs)
716 if (ys.le.0.d+0) then
717 omode=7
718 if (impres.ge.1) write (io,931) niter,ys
719 931 format (/" >>> m1qn3 (iteration ",i2,
720 & "): the scalar product (y,s) = ",d12.5
721 & /27x,"is not positive")
722 goto 1000
723 endif
724 c
725 c --- ybar et sbar
726 c
727 d1=dsqrt(1.d+0/ys)
728 do 410 i=1,n
729 sbar(i,jcour)=d1*sbar(i,jcour)
730 ybar(i,jcour)=d1*ybar(i,jcour)
731 410 continue
732 if (.not.inmemo) call ystbl (.true.,ybar,sbar,n,jmax)
733 c
734 c --- compute the scalar or diagonal preconditioner
735 c
736 if (impres.ge.5) write(io,932)
737 932 format (/" m1qn3: matrix update:")
738 c
739 c --- Here is the Oren-Spedicato factor, for scalar scaling
740 c
741 if (sscale) then
742 call prosca (n,ybar(1,jcour),ybar(1,jcour),ps,izs,rzs,dzs)
743 precos=1.d+0/ps
744 c
745 if (impres.ge.5) write (io,933) precos
746 933 format (5x,"Oren-Spedicato factor = ",d10.3)
747 c
748 c --- Scale the diagonal to Rayleigh s ellipsoid.
749 c Initially (niter.eq.1) and for a cold start, this is
750 c equivalent to an Oren-Spedicato scaling of the
751 c identity matrix.
752 c
753 else
754 call ctonb (n,ybar(1,jcour),aux,izs,rzs,dzs)
755 ps=0.d0
756 do 420 i=1,n
757 ps=ps+diag(i)*aux(i)*aux(i)
758 420 continue
759 d1=1.d0/ps
760 if (impres.ge.5) then
761 write (io,934) d1
762 934 format(5x,"fitting the ellipsoid: factor = ",1pd10.3)
763 endif
764 do 421 i=1,n
765 diag(i)=diag(i)*d1
766 421 continue
767 c
768 c --- update the diagonal
769 c (gg is used as an auxiliary vector)
770 c
771 call ctonb (n,sbar(1,jcour),gg,izs,rzs,dzs)
772 ps=0.d0
773 do 430 i=1,n
774 ps=ps+gg(i)*gg(i)/diag(i)
775 430 continue
776 den=ps
777 do 431 i=1,n
778 diag(i)=1.d0/
779 & (1.d0/diag(i)+aux(i)**2-(gg(i)/diag(i))**2/den)
780 if (diag(i).le.0.d0) then
781 if (impres.ge.5) write (io,935) i,diag(i),rmin
782 diag(i)=rmin
783 endif
784 431 continue
785 935 format (/" >>> m1qn3-WARNING: diagonal element ",i8,
786 & " is negative (",d10.3,"), reset to ",d10.3)
787 c
788 if (impres.ge.5) then
789 ps=0.d0
790 do 440 i=1,n
791 ps=ps+diag(i)
792 440 continue
793 ps=ps/n
794 preco=ps
795 c
796 ps2=0.d0
797 do 441 i=1,n
798 ps2=ps2+(diag(i)-ps)**2
799 441 continue
800 ps2=dsqrt(ps2/n)
801 write (io,936) preco,ps2
802 936 format (5x,"updated diagonal: average value = ",
803 & 1pd10.3,", sqrt(variance) = ",d10.3)
804 endif
805 endif
806 endif
807 c
808 c --- printings
809 c
810 c
811 c --- tests d arret
812 c
813 call prosca(n,g,g,ps,izs,rzs,dzs)
814 if (normtype.eq.'two') then
815 gnorm = sqrt(ddot(n,g,1,g,1))
816 elseif (normtype.eq.'sup') then
817 gnorm = dnrmi(n,g)
818 elseif (normtype.eq.'dfn') then
819 gnorm = dsqrt(ps)
820 endif
821 eps1 = gnorm/gnorms
822 c
823 if (impres.eq.3) then
824 if (mod(niter-1,50).eq.0) write(io,'(/a,a)')
825 & " iter simul stepsize f |g|",
826 & " |g|/|g0|"
827 write(io,
828 & '(1x,i5,2x,i5,2x,1pd8.2,2x,d21.14,2x,d11.5,2x,d10.4)')
829 & niter, isim, t, f, gnorm, eps1
830 endif
831 if (impres.ge.5) write (io,940) eps1
832 940 format (/" m1qn3: stopping criterion on g: ",1pd12.5)
833 if (eps1.lt.epsg) then
834 omode=1
835 goto 1000
836 endif
837 if (niter.eq.itmax) then
838 omode=4
839 if (impres.ge.1) write (io,941) niter
840 941 format (/" >>> m1qn3 (iteration ",i0,
841 & "): maximal number of iterations")
842 goto 1000
843 endif
844 if (isim.gt.nsim) then
845 omode=5
846 if (impres.ge.1) write (io,942) niter,isim
847 942 format (/" >>> m1qn3 (iteration ",i3,"): ",i6,
848 & " simulations (maximal number reached)")
849 goto 1000
850 endif
851 c
852 c --- calcul de la nouvelle direction de descente d = - H.g
853 c
854 if (m.eq.0) then
855 preco=2.d0*(ff-f)/ps
856 do 500 i=1,n
857 d(i)=-g(i)*preco
858 500 continue
859 else
860 do 510 i=1,n
861 d(i)=-g(i)
862 510 continue
863 if (inmemo) then
864 call dd (prosca,ctonb,ctcab,n,sscale,m,d,aux,jmin,jmax,
865 & precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
866 else
867 call dds (prosca,ctonb,ctcab,n,sscale,m,d,aux,jmin,jmax,
868 & precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
869 endif
870 endif
871 c
872 c --- test: la direction d est-elle de descente ?
873 c hp0 sera utilise par mlis3
874 c
875 call prosca (n,d,g,hp0,izs,rzs,dzs)
876 if (hp0.ge.0.d+0) then
877 omode=7
878 if (impres.ge.1) write (io,905) niter,hp0
879 goto 1000
880 endif
881 if (impres.ge.5) then
882 call prosca (n,g,g,ps,izs,rzs,dzs)
883 ps=dsqrt(ps)
884 call prosca (n,d,d,ps2,izs,rzs,dzs)
885 ps2=dsqrt(ps2)
886 ps=hp0/ps/ps2
887 ps=dmin1(-ps,1.d+0)
888 ps=dacos(ps)
889 d1=ps
890 d1=d1*180.d0/pi
891 write (io,906) sngl(d1)
892 endif
893 c
894 c---- on poursuit les iterations
895 c
896 goto 100
897 c
898 c --- n1qn3 has finished for ever
899 c
900 1000 continue
901 if (reverse.ne.0) reverse = -1
902 reentry = 0
903 nsim=isim
904 epsg=eps1
905 return
906 end
907 c
908 c--------0---------0---------0---------0---------0---------0---------0--
909 c
910 subroutine dd (prosca,ctonb,ctcab,n,sscale,nm,depl,aux,jmin,jmax,
911 & precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
912 c----
913 c
914 c calcule le produit H.g ou
915 c . H est une matrice construite par la formule de bfgs inverse
916 c a nm memoires a partir de la matrice diagonale diag
917 c dans un espace hilbertien dont le produit scalaire
918 c est donne par prosca
919 c (cf. J. Nocedal, Math. of Comp. 35/151 (1980) 773-782)
920 c . g est un vecteur de dimension n (en general le gradient)
921 c
922 c la matrice diag apparait donc comme un preconditionneur diagonal
923 c
924 c depl = g (en entree), = H g (en sortie)
925 c
926 c la matrice H est memorisee par les vecteurs des tableaux
927 c ybar, sbar et les pointeurs jmin, jmax
928 c
929 c alpha(nm) est une zone de travail
930 c
931 c izs(*),rzs(*),dzs(*) sont des zones de travail pour prosca
932 c
933 c----
934 c
935 implicit none
936 c
937 c arguments
938 c
939 logical sscale
940 integer n,nm,jmin,jmax,izs(*)
941 real rzs(*)
942 double precision depl(n),precos,diag(n),alpha(nm),ybar(n,1),
943 & sbar(n,1),aux(n),dzs(*)
944 external prosca,ctonb,ctcab
945 c
946 c variables locales
947 c
948 integer jfin,i,j,jp
949 double precision r,ps
950 c
951 jfin=jmax
952 if (jfin.lt.jmin) jfin=jmax+nm
953 c
954 c phase de descente
955 c
956 do 100 j=jfin,jmin,-1
957 jp=j
958 if (jp.gt.nm) jp=jp-nm
959 call prosca (n,depl,sbar(1,jp),ps,izs,rzs,dzs)
960 r=ps
961 alpha(jp)=r
962 do 20 i=1,n
963 depl(i)=depl(i)-r*ybar(i,jp)
964 20 continue
965 100 continue
966 c
967 c preconditionnement
968 c
969 if (sscale) then
970 do 150 i=1,n
971 depl(i)=depl(i)*precos
972 150 continue
973 else
974 call ctonb (n,depl,aux,izs,rzs,dzs)
975 do 151 i=1,n
976 aux(i)=aux(i)*diag(i)
977 151 continue
978 call ctcab (n,aux,depl,izs,rzs,dzs)
979 endif
980 c
981 c remontee
982 c
983 do 200 j=jmin,jfin
984 jp=j
985 if (jp.gt.nm) jp=jp-nm
986 call prosca (n,depl,ybar(1,jp),ps,izs,rzs,dzs)
987 r=alpha(jp)-ps
988 do 120 i=1,n
989 depl(i)=depl(i)+r*sbar(i,jp)
990 120 continue
991 200 continue
992 return
993 end
994 c
995 c--------0---------0---------0---------0---------0---------0---------0--
996 c
997 subroutine dds (prosca,ctonb,ctcab,n,sscale,nm,depl,aux,jmin,
998 & jmax,precos,diag,alpha,ybar,sbar,izs,rzs,dzs)
999 c----
1000 c
1001 c This subroutine has the same role as dd (computation of the
1002 c product H.g). It supposes however that the (y,s) pairs are not
1003 c stored in core memory, but on a devise chosen by the user.
1004 c The access to this devise is performed via the subroutine ystbl.
1005 c
1006 c----
1007 c
1008 implicit none
1009 c
1010 c arguments
1011 c
1012 logical sscale
1013 integer n,nm,jmin,jmax,izs(*)
1014 real rzs(*)
1015 double precision depl(n),precos,diag(n),alpha(nm),ybar(n),sbar(n),
1016 & aux(n),dzs(*)
1017 external prosca,ctonb,ctcab
1018 c
1019 c variables locales
1020 c
1021 integer jfin,i,j,jp
1022 double precision r,ps
1023 c
1024 jfin=jmax
1025 if (jfin.lt.jmin) jfin=jmax+nm
1026 c
1027 c phase de descente
1028 c
1029 do 100 j=jfin,jmin,-1
1030 jp=j
1031 if (jp.gt.nm) jp=jp-nm
1032 call ystbl (.false.,ybar,sbar,n,jp)
1033 call prosca (n,depl,sbar,ps,izs,rzs,dzs)
1034 r=ps
1035 alpha(jp)=r
1036 do 20 i=1,n
1037 depl(i)=depl(i)-r*ybar(i)
1038 20 continue
1039 100 continue
1040 c
1041 c preconditionnement
1042 c
1043 if (sscale) then
1044 do 150 i=1,n
1045 depl(i)=depl(i)*precos
1046 150 continue
1047 else
1048 call ctonb (n,depl,aux,izs,rzs,dzs)
1049 do 151 i=1,n
1050 aux(i)=aux(i)*diag(i)
1051 151 continue
1052 call ctcab (n,aux,depl,izs,rzs,dzs)
1053 endif
1054 c
1055 c remontee
1056 c
1057 do 200 j=jmin,jfin
1058 jp=j
1059 if (jp.gt.nm) jp=jp-nm
1060 call ystbl (.false.,ybar,sbar,n,jp)
1061 call prosca (n,depl,ybar(1),ps,izs,rzs,dzs)
1062 r=alpha(jp)-ps
1063 do 120 i=1,n
1064 depl(i)=depl(i)+r*sbar(i)
1065 120 continue
1066 200 continue
1067 return
1068 end
1069 c
1070 c--------0---------0---------0---------0---------0---------0---------0--
1071 c
1072 subroutine mlis3 (n,simul,prosca,x,f,fpn,t,tmin,tmax,d,g,
1073 1 amd,amf,imp,io,logic,nap,napmax,xn,
1074 1 reverse,reentry,indic,izs,rzs,dzs)
1075 c
1076 c ----
1077 c
1078 c mlis3 + minuscules + commentaires
1079 c + version amelioree (XII 88): interpolation cubique systematique
1080 c et anti-overflows
1081 c + declaration variables (II/89, JCG).
1082 c + barr is also progressively decreased (12/93, CL & JChG).
1083 c barmul is set to 5.
1084 c
1085 c ----------------------------------------------------------------
1086 c
1087 c en sortie logic =
1088 c
1089 c 0 descente serieuse
1090 c 1 descente bloquee sur tmax
1091 c 4 nap > napmax
1092 c 5 arret demande par le simulateur
1093 c 6 fonction et gradient pas d accord
1094 c < 0 contrainte implicite active
1095 c
1096 c indic on entry (in reverse communication)
1097 c
1098 c < 0: the simulator cannot compute f and g
1099 c = 0: the simulator wants to stop
1100 c > 0: the simulator has done its work
1101 c
1102 c indic on return (in reverse communication)
1103 c
1104 c = 4: the simulator is asked to compute f and g
1105 c
1106 c reverse
1107 c
1108 c = 0: direct communication
1109 c = 1: reverse communication
1110 c
1111 c reentry on entry
1112 c
1113 c = 2: reverse communication, return from a simulation, skip to
1114 c the place where the simulator was called (9999)
1115 c
1116 c reentry return
1117 c
1118 c = 0: reverse communication, the linesearch has finished its job
1119 c = 2: reverse communication, a simulation is required
1120 c
1121 c ----
1122 c
1123 implicit none
1124 c
1125 c --- arguments
1126 c
1127 integer n,imp,io,logic,nap,napmax,indic,reverse,reentry,izs(*)
1128 real rzs(*)
1129 double precision x(n),f,fpn,t,tmin,tmax,d(n),g(n),amd,amf,xn(n),
1130 & dzs(*)
1131 external simul,prosca
1132 c
1133 c --- variables locales
1134 c
1135 #include "mlis3_common.h"
1136 integer i
1137 double precision taa,ps
1138 CML save t_increased,i,indica,indicd,tesf,tesd,tg,fg,fpg,td,ta,fa,
1139 CML & fpa,d2,fn,fp,ffn,fd,fpd,z,test,barmin,barmul,barmax,barr,
1140 CML & gauche,droite,taa,ps
1141 c
1142 1000 format (/4x," mlis3 ",4x,"fpn=",1pd10.3," d2=",d9.2,
1143 1 " tmin=",d9.2," tmax=",d9.2)
1144 1001 format (/4x," mlis3",3x,"stop on tmin",5x,
1145 1 "stepsizes",11x,"functions",8x,"derivatives")
1146 1002 format (4x," mlis3",37x,1pd10.3,2d11.3)
1147 1003 format (4x," mlis3",1pd14.3,2d11.3)
1148 1004 format (4x," mlis3",37x,1pd10.3," indic=",i3)
1149 1005 format (4x," mlis3",14x,1pd18.8,2x,d21.14,2x,d11.4)
1150 1006 format (4x," mlis3",14x,1pd18.8," indic=",i3)
1151 1007 format (/4x," mlis3",10x,"tmin forced to tmax")
1152 1008 format (/4x," mlis3",10x,"inconsistent call")
1153 c
1154 c --- possible jump
1155 c
1156 if (reentry.eq.2) goto 9999
1157 c
1158 if (n.gt.0 .and. fpn.lt.0.d0 .and. t.gt.0.d0
1159 1 .and. tmax.gt.0.d0 .and. amf.gt.0.d0
1160 1 .and. amd.gt.amf .and. amd.lt.1.d0) go to 5
1161 logic=6
1162 go to 999
1163 5 tesf=amf*fpn
1164 tesd=amd*fpn
1165 barmin=0.01d0
1166 barmul=5.d0
1167 barmax=0.3d0
1168 barr=barmin
1169 td=0.d0
1170 tg=0.d0
1171 fn=f
1172 fg=fn
1173 fpg=fpn
1174 ta=0.d0
1175 fa=fn
1176 fpa=fpn
1177 call prosca (n,d,d,ps,izs,rzs,dzs)
1178 d2=ps
1179 c
1180 c elimination d un t initial ridiculement petit
1181 c
1182 c<
1183 c if (t.gt.tmin) go to 20
1184 c t=tmin
1185 c if (t.le.tmax) go to 20
1186 c if (imp.gt.0) write (io,1007)
1187 c tmin=tmax
1188 c 20 if (fn+t*fpn.lt.fn+0.9d0*t*fpn) go to 30
1189 c t=2.d0*t
1190 c go to 20
1191 c changed into
1192 if (t.lt.tmin) then
1193 t=tmin
1194 if (imp.ge.4) write (io,'(a)')
1195 & ' mlis3: initial step-size increased to tmin'
1196 if (t.gt.tmax) then
1197 if (imp.gt.0) write (io,1007)
1198 tmin=tmax
1199 endif
1200 endif
1201 t_increased = .false.
1202 do while (fn+t*fpn.ge.fn+0.9d0*t*fpn)
1203 t_increased = .true.
1204 t = 2.d0*t
1205 enddo
1206 if (t_increased .and. (imp.ge.4)) write (io,'(a,1pd10.3)')
1207 & ' mlis3: initial step-size increased to ',t
1208 c>
1209 30 indica=1
1210 logic=0
1211 if (t.gt.tmax) then
1212 t=tmax
1213 logic=1
1214 endif
1215 if (imp.ge.4) write (io,1000) fpn,d2,tmin,tmax
1216 c
1217 c --- nouveau x
1218 c initialize xn to the current iterate
1219 c use x as the trial iterate
1220 c
1221 do i=1,n
1222 xn(i)=x(i)
1223 x(i)=xn(i)+t*d(i)
1224 enddo
1225 c
1226 c --- boucle
1227 c
1228 100 nap=nap+1
1229 if(nap.gt.napmax) then
1230 logic=4
1231 fn=fg
1232 do i=1,n
1233 xn(i)=xn(i)+tg*d(i)
1234 enddo
1235 go to 999
1236 endif
1237 c
1238 c --- appel simulateur
1239 c
1240 indic=4
1241 if (reverse.gt.0) then
1242 reentry = 2
1243 return
1244 else
1245 call simul(indic,n,x,f,g,izs,rzs,dzs)
1246 endif
1247 9999 continue
1248 if (indic.eq.0) then
1249 c
1250 c --- arret demande par l utilisateur
1251 c
1252 logic=5
1253 fn=f
1254 do 170 i=1,n
1255 xn(i)=x(i)
1256 170 continue
1257 go to 999
1258 endif
1259 if(indic.lt.0) then
1260 c
1261 c --- les calculs n ont pas pu etre effectues par le simulateur
1262 c
1263 td=t
1264 indicd=indic
1265 logic=0
1266 if (imp.ge.4) write (io,1004) t,indic
1267 t=tg+0.1d0*(td-tg)
1268 go to 905
1269 endif
1270 c
1271 c --- les tests elementaires sont faits, on y va
1272 c
1273 call prosca (n,d,g,ps,izs,rzs,dzs)
1274 fp=ps
1275 c
1276 c --- premier test de Wolfe
1277 c
1278 ffn=f-fn
1279 if(ffn.gt.t*tesf) then
1280 td=t
1281 fd=f
1282 fpd=fp
1283 indicd=indic
1284 logic=0
1285 if(imp.ge.4) write (io,1002) t,ffn,fp
1286 go to 500
1287 endif
1288 c
1289 c --- test 1 ok, donc deuxieme test de Wolfe
1290 c
1291 if(imp.ge.4) write (io,1003) t,ffn,fp
1292 if(fp.gt.tesd) then
1293 logic=0
1294 go to 320
1295 endif
1296 if (logic.eq.0) go to 350
1297 c
1298 c --- test 2 ok, donc pas serieux, on sort
1299 c
1300 320 fn=f
1301 do 330 i=1,n
1302 xn(i)=x(i)
1303 330 continue
1304 go to 999
1305 c
1306 c
1307 c
1308 350 tg=t
1309 fg=f
1310 fpg=fp
1311 if(td.ne.0.d0) go to 500
1312 c
1313 c extrapolation
1314 c
1315 taa=t
1316 gauche=(1.d0+barmin)*t
1317 droite=10.d0*t
1318 call ecube (t,f,fp,ta,fa,fpa,gauche,droite)
1319 ta=taa
1320 if(t.lt.tmax) go to 900
1321 logic=1
1322 t=tmax
1323 go to 900
1324 c
1325 c interpolation
1326 c
1327 500 if(indica.le.0) then
1328 ta=t
1329 t=0.9d0*tg+0.1d0*td
1330 go to 900
1331 endif
1332 test=barr*(td-tg)
1333 gauche=tg+test
1334 droite=td-test
1335 taa=t
1336 call ecube (t,f,fp,ta,fa,fpa,gauche,droite)
1337 ta=taa
1338 if (t.gt.gauche .and. t.lt.droite) then
1339 barr=dmax1(barmin,barr/barmul)
1340 c barr=barmin
1341 else
1342 barr=dmin1(barmul*barr,barmax)
1343 endif
1344 c
1345 c --- fin de boucle
1346 c - t peut etre bloque sur tmax
1347 c (venant de l extrapolation avec logic=1)
1348 c
1349 900 fa=f
1350 fpa=fp
1351 905 indica=indic
1352 c
1353 c --- faut-il continuer ?
1354 c
1355 if (td.eq.0.d0) go to 950
1356 if (td-tg.lt.tmin) go to 920
1357 c
1358 c --- limite de precision machine (arret de secours) ?
1359 c
1360 do 910 i=1,n
1361 z=xn(i)+t*d(i)
1362 if (z.ne.xn(i).and.z.ne.x(i)) go to 950
1363 910 continue
1364 if (imp.gt.3) write (io,'(5x,a)') "mlis3 no change in x"
1365 c
1366 c --- arret sur dxmin ou de secours
1367 c
1368 920 logic=6
1369 c
1370 c si indicd<0, derniers calculs non faits par simul
1371 c
1372 if (indicd.lt.0) logic=indicd
1373 c
1374 c si tg=0, xn = xn_depart,
1375 c sinon on prend xn=x_gauche qui fait decroitre f
1376 c
1377 if (tg.eq.0.d0) go to 940
1378 fn=fg
1379 do 930 i=1,n
1380 930 xn(i)=xn(i)+tg*d(i)
1381 940 if (imp.le.3) go to 999
1382 write (io,1001)
1383 write (io,1005) tg,fg,fpg
1384 if (logic.eq.6) write (io,1005) td,fd,fpd
1385 if (logic.eq.7) write (io,1006) td,indicd
1386 go to 999
1387 c
1388 c recopiage de x et boucle
1389 c
1390 950 do 960 i=1,n
1391 960 x(i)=xn(i)+t*d(i)
1392 go to 100
1393 c --- linesearch finished, no skip at next entry in mlis3
1394 999 if (reverse.ne.0) reentry = 0
1395 return
1396 end
1397 c
1398 c
1399 c--------0---------0---------0---------0---------0---------0---------0--
1400 c
1401 subroutine ecube (t,f,fp,ta,fa,fpa,tlower,tupper)
1402 c
1403 implicit none
1404 c
1405 c --- arguments
1406 c
1407 double precision sign,den,anum,t,f,fp,ta,fa,fpa,tlower,tupper
1408 c
1409 c --- variables locales
1410 c
1411 double precision z1,b,discri
1412 c
1413 c Using f and fp at t and ta, computes new t by cubic formula
1414 c safeguarded inside [tlower,tupper].
1415 c
1416 z1=fp+fpa-3.d0*(fa-f)/(ta-t)
1417 b=z1+fp
1418 c
1419 c first compute the discriminant (without overflow)
1420 c
1421 if (dabs(z1).le.1.d0) then
1422 discri=z1*z1-fp*fpa
1423 else
1424 discri=fp/z1
1425 discri=discri*fpa
1426 discri=z1-discri
1427 if (z1.ge.0.d0 .and. discri.ge.0.d0) then
1428 discri=dsqrt(z1)*dsqrt(discri)
1429 go to 120
1430 endif
1431 if (z1.le.0.d0 .and. discri.le.0.d0) then
1432 discri=dsqrt(-z1)*dsqrt(-discri)
1433 go to 120
1434 endif
1435 discri=-1.d0
1436 endif
1437 if (discri.lt.0.d0) then
1438 if (fp.lt.0.d0) t=tupper
1439 if (fp.ge.0.d0) t=tlower
1440 go to 900
1441 endif
1442 c
1443 c discriminant nonnegative, compute solution (without overflow)
1444 c
1445 discri=dsqrt(discri)
1446 120 if (t-ta.lt.0.d0) discri=-discri
1447 sign=(t-ta)/dabs(t-ta)
1448 if (b*sign.gt.0.d+0) then
1449 t=t+fp*(ta-t)/(b+discri)
1450 else
1451 den=z1+b+fpa
1452 anum=b-discri
1453 if (dabs((t-ta)*anum).lt.(tupper-tlower)*dabs(den)) then
1454 t=t+anum*(ta-t)/den
1455 else
1456 t=tupper
1457 endif
1458 endif
1459 900 t=dmax1(t,tlower)
1460 t=dmin1(t,tupper)
1461 return
1462 end
1463 c
1464 c--------0---------0---------0---------0---------0---------0---------0--
1465 c
1466 subroutine mupdts (sscale,inmemo,n,m,nrz)
1467 c
1468 implicit none
1469 c
1470 c arguments
1471 c
1472 logical sscale,inmemo
1473 integer n,m,nrz
1474 c----
1475 c
1476 c On entry:
1477 c sscale: .true. if scalar initial scaling,
1478 c .false. if diagonal initial scaling
1479 c n: number of variables
1480 c
1481 c This routine has to return:
1482 c m: the number of updates to form the approximate Hessien H,
1483 c inmemo: .true., if the vectors y and s used to form H are stored
1484 c in core memory,
1485 c .false. otherwise (storage of y and s on disk, for
1486 c instance).
1487 c When inmemo=.false., the routine "ystbl", which stores and
1488 c restores (y,s) pairs, has to be rewritten.
1489 c
1490 c----
1491 c
1492 if (sscale) then
1493 m=(nrz-3*n)/(2*n+1)
1494 else
1495 m=(nrz-4*n)/(2*n+1)
1496 endif
1497 inmemo=.true.
1498 return
1499 end
1500 c
1501 c--------0---------0---------0---------0---------0---------0---------0--
1502 c
1503 subroutine ystbl (store,ybar,sbar,n,j)
1504 c----
1505 c
1506 c This subroutine should store (if store = .true.) or restore
1507 c (if store = .false.) a pair (ybar,sbar) at or from position
1508 c j in memory. Be sure to have 1 <= j <= m, where m in the number
1509 c of updates specified by subroutine mupdts.
1510 c
1511 c The subroutine is used only when the (y,s) pairs are not
1512 c stored in core memory in the arrays ybar(.,.) and sbar(.,.).
1513 c In this case, the subroutine has to be written by the user.
1514 c
1515 c----
1516 c
1517 implicit none
1518 c
1519 c arguments
1520 c
1521 logical store
1522 integer n,j
1523 double precision ybar(n),sbar(n)
1524 c
1525 return
1526 end
1527 c
1528 c--------0---------0---------0---------0---------0---------0---------0--
1529 c
1530 subroutine ctonbe (n,u,v,izs,rzs,dzs)
1531 c
1532 implicit none
1533 integer n,izs(*)
1534 real rzs(*)
1535 double precision u(1),v(1),dzs(*)
1536 c
1537 integer i
1538 c
1539 do i=1,n
1540 v(i)=u(i)
1541 enddo
1542 return
1543 end
1544 c
1545 c--------0---------0---------0---------0---------0---------0---------0--
1546 c
1547 subroutine ctcabe (n,u,v,izs,rzs,dzs)
1548 c
1549 implicit none
1550 integer n,izs(*)
1551 real rzs(*)
1552 double precision u(1),v(1),dzs(*)
1553 c
1554 integer i
1555 c
1556 do i=1,n
1557 v(i)=u(i)
1558 enddo
1559 return
1560 end
1561 c
1562 c--------0---------0---------0---------0---------0---------0---------0--
1563 c
1564 subroutine euclid (n,x,y,ps,izs,rzs,dzs)
1565 c
1566 implicit none
1567 integer n,izs(*)
1568 real rzs(*)
1569 double precision x(n),y(n),ps,dzs(*)
1570 c
1571 integer i
1572 c
1573 ps=0.d0
1574 do i=1,n
1575 ps=ps+x(i)*y(i)
1576 enddo
1577 return
1578 end
1579 c
1580 c--------0---------0---------0---------0---------0---------0---------0--
1581 c
1582 subroutine simul_rc (indic,n,x,f,g,izs,rzs,dzs)
1583 c
1584 implicit none
1585 integer indic,n,izs(*)
1586 real rzs(*)
1587 double precision x(n),f,g(n),dzs(*)
1588 c
1589 return
1590 end
1591 c
1592 c--------0---------0---------0---------0---------0---------0---------0--
1593 c
1594 double precision function dnrmi (n,v)
1595 c
1596 integer n
1597 double precision v(n)
1598 c
1599 c----
1600 c
1601 c Commputes the infinty-norm of the vector v(n)
1602 c
1603 c----
1604 c
1605 c --- local variables
1606 c
1607 integer i
1608 double precision norm
1609 c
1610 c --- compute
1611 c
1612 norm = 0.d0
1613 if (n.gt.0) then
1614 do i=1,n
1615 norm = max(norm,abs(v(i)))
1616 enddo
1617 endif
1618 dnrmi = norm
1619 c
1620 return
1621 end

  ViewVC Help
Powered by ViewVC 1.1.22