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

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

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


Revision 1.2 - (hide annotations) (download)
Fri Apr 27 12:00:50 2012 UTC (11 years, 11 months ago) by mlosch
Branch: MAIN
Changes since 1.1: +14 -2 lines
- store more of the mlis3 variables (probably too many now, but better
save than sorry). This does change the results during longer line
searches
- add more information and header information

1 mlosch 1.2 C $Header: $
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 program that are called alternatingly from a (shell-)script.
7     C This requires that the "state" of m1qn3 is saved at 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 mlosch 1.1 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,1),sbar(n,1),dzs(*)
433     external simul,prosca,ctonb,ctcab
434     c
435     c variables locales
436     c
437     logical skip_update
438     integer i,moderl
439     double precision preco,precos,ys,den,dk1,ps,ps2,hp0
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     logical t_increased
1137     integer i,indica,indicd
1138 mlosch 1.2 double precision taa,ps
1139 mlosch 1.1 CML save t_increased,i,indica,indicd,tesf,tesd,tg,fg,fpg,td,ta,fa,
1140     CML & fpa,d2,fn,fp,ffn,fd,fpd,z,test,barmin,barmul,barmax,barr,
1141     CML & gauche,droite,taa,ps
1142     c
1143     1000 format (/4x," mlis3 ",4x,"fpn=",1pd10.3," d2=",d9.2,
1144     1 " tmin=",d9.2," tmax=",d9.2)
1145     1001 format (/4x," mlis3",3x,"stop on tmin",5x,
1146     1 "stepsizes",11x,"functions",8x,"derivatives")
1147     1002 format (4x," mlis3",37x,1pd10.3,2d11.3)
1148     1003 format (4x," mlis3",1pd14.3,2d11.3)
1149     1004 format (4x," mlis3",37x,1pd10.3," indic=",i3)
1150     1005 format (4x," mlis3",14x,1pd18.8,2x,d21.14,2x,d11.4)
1151     1006 format (4x," mlis3",14x,1pd18.8," indic=",i3)
1152     1007 format (/4x," mlis3",10x,"tmin forced to tmax")
1153     1008 format (/4x," mlis3",10x,"inconsistent call")
1154     c
1155     c --- possible jump
1156     c
1157     if (reentry.eq.2) goto 9999
1158     c
1159     if (n.gt.0 .and. fpn.lt.0.d0 .and. t.gt.0.d0
1160     1 .and. tmax.gt.0.d0 .and. amf.gt.0.d0
1161     1 .and. amd.gt.amf .and. amd.lt.1.d0) go to 5
1162     logic=6
1163     go to 999
1164     5 tesf=amf*fpn
1165     tesd=amd*fpn
1166     barmin=0.01d0
1167     barmul=5.d0
1168     barmax=0.3d0
1169     barr=barmin
1170     td=0.d0
1171     tg=0.d0
1172     fn=f
1173     fg=fn
1174     fpg=fpn
1175     ta=0.d0
1176     fa=fn
1177     fpa=fpn
1178     call prosca (n,d,d,ps,izs,rzs,dzs)
1179     d2=ps
1180     c
1181     c elimination d un t initial ridiculement petit
1182     c
1183     c<
1184     c if (t.gt.tmin) go to 20
1185     c t=tmin
1186     c if (t.le.tmax) go to 20
1187     c if (imp.gt.0) write (io,1007)
1188     c tmin=tmax
1189     c 20 if (fn+t*fpn.lt.fn+0.9d0*t*fpn) go to 30
1190     c t=2.d0*t
1191     c go to 20
1192     c changed into
1193     if (t.lt.tmin) then
1194     t=tmin
1195     if (imp.ge.4) write (io,'(a)')
1196     & ' mlis3: initial step-size increased to tmin'
1197     if (t.gt.tmax) then
1198     if (imp.gt.0) write (io,1007)
1199     tmin=tmax
1200     endif
1201     endif
1202     t_increased = .false.
1203     do while (fn+t*fpn.ge.fn+0.9d0*t*fpn)
1204     t_increased = .true.
1205     t = 2.d0*t
1206     enddo
1207     if (t_increased .and. (imp.ge.4)) write (io,'(a,1pd10.3)')
1208     & ' mlis3: initial step-size increased to ',t
1209     c>
1210     30 indica=1
1211     logic=0
1212     if (t.gt.tmax) then
1213     t=tmax
1214     logic=1
1215     endif
1216     if (imp.ge.4) write (io,1000) fpn,d2,tmin,tmax
1217     c
1218     c --- nouveau x
1219     c initialize xn to the current iterate
1220     c use x as the trial iterate
1221     c
1222     do i=1,n
1223     xn(i)=x(i)
1224     x(i)=xn(i)+t*d(i)
1225     enddo
1226     c
1227     c --- boucle
1228     c
1229     100 nap=nap+1
1230     if(nap.gt.napmax) then
1231     logic=4
1232     fn=fg
1233     do i=1,n
1234     xn(i)=xn(i)+tg*d(i)
1235     enddo
1236     go to 999
1237     endif
1238     c
1239     c --- appel simulateur
1240     c
1241     indic=4
1242     if (reverse.gt.0) then
1243     reentry = 2
1244     return
1245     else
1246     call simul(indic,n,x,f,g,izs,rzs,dzs)
1247     endif
1248     9999 continue
1249     if (indic.eq.0) then
1250     c
1251     c --- arret demande par l utilisateur
1252     c
1253     logic=5
1254     fn=f
1255     do 170 i=1,n
1256     xn(i)=x(i)
1257     170 continue
1258     go to 999
1259     endif
1260     if(indic.lt.0) then
1261     c
1262     c --- les calculs n ont pas pu etre effectues par le simulateur
1263     c
1264     td=t
1265     indicd=indic
1266     logic=0
1267     if (imp.ge.4) write (io,1004) t,indic
1268     t=tg+0.1d0*(td-tg)
1269     go to 905
1270     endif
1271     c
1272     c --- les tests elementaires sont faits, on y va
1273     c
1274     call prosca (n,d,g,ps,izs,rzs,dzs)
1275     fp=ps
1276     c
1277     c --- premier test de Wolfe
1278     c
1279     ffn=f-fn
1280     if(ffn.gt.t*tesf) then
1281     td=t
1282     fd=f
1283     fpd=fp
1284     indicd=indic
1285     logic=0
1286     if(imp.ge.4) write (io,1002) t,ffn,fp
1287     go to 500
1288     endif
1289     c
1290     c --- test 1 ok, donc deuxieme test de Wolfe
1291     c
1292     if(imp.ge.4) write (io,1003) t,ffn,fp
1293     if(fp.gt.tesd) then
1294     logic=0
1295     go to 320
1296     endif
1297     if (logic.eq.0) go to 350
1298     c
1299     c --- test 2 ok, donc pas serieux, on sort
1300     c
1301     320 fn=f
1302     do 330 i=1,n
1303     xn(i)=x(i)
1304     330 continue
1305     go to 999
1306     c
1307     c
1308     c
1309     350 tg=t
1310     fg=f
1311     fpg=fp
1312     if(td.ne.0.d0) go to 500
1313     c
1314     c extrapolation
1315     c
1316     taa=t
1317     gauche=(1.d0+barmin)*t
1318     droite=10.d0*t
1319     call ecube (t,f,fp,ta,fa,fpa,gauche,droite)
1320     ta=taa
1321     if(t.lt.tmax) go to 900
1322     logic=1
1323     t=tmax
1324     go to 900
1325     c
1326     c interpolation
1327     c
1328     500 if(indica.le.0) then
1329     ta=t
1330     t=0.9d0*tg+0.1d0*td
1331     go to 900
1332     endif
1333     test=barr*(td-tg)
1334     gauche=tg+test
1335     droite=td-test
1336     taa=t
1337     call ecube (t,f,fp,ta,fa,fpa,gauche,droite)
1338     ta=taa
1339     if (t.gt.gauche .and. t.lt.droite) then
1340     barr=dmax1(barmin,barr/barmul)
1341     c barr=barmin
1342     else
1343     barr=dmin1(barmul*barr,barmax)
1344     endif
1345     c
1346     c --- fin de boucle
1347     c - t peut etre bloque sur tmax
1348     c (venant de l extrapolation avec logic=1)
1349     c
1350     900 fa=f
1351     fpa=fp
1352     905 indica=indic
1353     c
1354     c --- faut-il continuer ?
1355     c
1356     if (td.eq.0.d0) go to 950
1357     if (td-tg.lt.tmin) go to 920
1358     c
1359     c --- limite de precision machine (arret de secours) ?
1360     c
1361     do 910 i=1,n
1362     z=xn(i)+t*d(i)
1363     if (z.ne.xn(i).and.z.ne.x(i)) go to 950
1364     910 continue
1365     if (imp.gt.3) write (io,'(5x,a)') "mlis3 no change in x"
1366     c
1367     c --- arret sur dxmin ou de secours
1368     c
1369     920 logic=6
1370     c
1371     c si indicd<0, derniers calculs non faits par simul
1372     c
1373     if (indicd.lt.0) logic=indicd
1374     c
1375     c si tg=0, xn = xn_depart,
1376     c sinon on prend xn=x_gauche qui fait decroitre f
1377     c
1378     if (tg.eq.0.d0) go to 940
1379     fn=fg
1380     do 930 i=1,n
1381     930 xn(i)=xn(i)+tg*d(i)
1382     940 if (imp.le.3) go to 999
1383     write (io,1001)
1384     write (io,1005) tg,fg,fpg
1385     if (logic.eq.6) write (io,1005) td,fd,fpd
1386     if (logic.eq.7) write (io,1006) td,indicd
1387     go to 999
1388     c
1389     c recopiage de x et boucle
1390     c
1391     950 do 960 i=1,n
1392     960 x(i)=xn(i)+t*d(i)
1393     go to 100
1394     c --- linesearch finished, no skip at next entry in mlis3
1395     999 if (reverse.ne.0) reentry = 0
1396     return
1397     end
1398     c
1399     c
1400     c--------0---------0---------0---------0---------0---------0---------0--
1401     c
1402     subroutine ecube (t,f,fp,ta,fa,fpa,tlower,tupper)
1403     c
1404     implicit none
1405     c
1406     c --- arguments
1407     c
1408     double precision sign,den,anum,t,f,fp,ta,fa,fpa,tlower,tupper
1409     c
1410     c --- variables locales
1411     c
1412     double precision z1,b,discri
1413     c
1414     c Using f and fp at t and ta, computes new t by cubic formula
1415     c safeguarded inside [tlower,tupper].
1416     c
1417     z1=fp+fpa-3.d0*(fa-f)/(ta-t)
1418     b=z1+fp
1419     c
1420     c first compute the discriminant (without overflow)
1421     c
1422     if (dabs(z1).le.1.d0) then
1423     discri=z1*z1-fp*fpa
1424     else
1425     discri=fp/z1
1426     discri=discri*fpa
1427     discri=z1-discri
1428     if (z1.ge.0.d0 .and. discri.ge.0.d0) then
1429     discri=dsqrt(z1)*dsqrt(discri)
1430     go to 120
1431     endif
1432     if (z1.le.0.d0 .and. discri.le.0.d0) then
1433     discri=dsqrt(-z1)*dsqrt(-discri)
1434     go to 120
1435     endif
1436     discri=-1.d0
1437     endif
1438     if (discri.lt.0.d0) then
1439     if (fp.lt.0.d0) t=tupper
1440     if (fp.ge.0.d0) t=tlower
1441     go to 900
1442     endif
1443     c
1444     c discriminant nonnegative, compute solution (without overflow)
1445     c
1446     discri=dsqrt(discri)
1447     120 if (t-ta.lt.0.d0) discri=-discri
1448     sign=(t-ta)/dabs(t-ta)
1449     if (b*sign.gt.0.d+0) then
1450     t=t+fp*(ta-t)/(b+discri)
1451     else
1452     den=z1+b+fpa
1453     anum=b-discri
1454     if (dabs((t-ta)*anum).lt.(tupper-tlower)*dabs(den)) then
1455     t=t+anum*(ta-t)/den
1456     else
1457     t=tupper
1458     endif
1459     endif
1460     900 t=dmax1(t,tlower)
1461     t=dmin1(t,tupper)
1462     return
1463     end
1464     c
1465     c--------0---------0---------0---------0---------0---------0---------0--
1466     c
1467     subroutine mupdts (sscale,inmemo,n,m,nrz)
1468     c
1469     implicit none
1470     c
1471     c arguments
1472     c
1473     logical sscale,inmemo
1474     integer n,m,nrz
1475     c----
1476     c
1477     c On entry:
1478     c sscale: .true. if scalar initial scaling,
1479     c .false. if diagonal initial scaling
1480     c n: number of variables
1481     c
1482     c This routine has to return:
1483     c m: the number of updates to form the approximate Hessien H,
1484     c inmemo: .true., if the vectors y and s used to form H are stored
1485     c in core memory,
1486     c .false. otherwise (storage of y and s on disk, for
1487     c instance).
1488     c When inmemo=.false., the routine "ystbl", which stores and
1489     c restores (y,s) pairs, has to be rewritten.
1490     c
1491     c----
1492     c
1493     if (sscale) then
1494     m=(nrz-3*n)/(2*n+1)
1495     else
1496     m=(nrz-4*n)/(2*n+1)
1497     endif
1498     inmemo=.true.
1499     return
1500     end
1501     c
1502     c--------0---------0---------0---------0---------0---------0---------0--
1503     c
1504     subroutine ystbl (store,ybar,sbar,n,j)
1505     c----
1506     c
1507     c This subroutine should store (if store = .true.) or restore
1508     c (if store = .false.) a pair (ybar,sbar) at or from position
1509     c j in memory. Be sure to have 1 <= j <= m, where m in the number
1510     c of updates specified by subroutine mupdts.
1511     c
1512     c The subroutine is used only when the (y,s) pairs are not
1513     c stored in core memory in the arrays ybar(.,.) and sbar(.,.).
1514     c In this case, the subroutine has to be written by the user.
1515     c
1516     c----
1517     c
1518     implicit none
1519     c
1520     c arguments
1521     c
1522     logical store
1523     integer n,j
1524     double precision ybar(n),sbar(n)
1525     c
1526     return
1527     end
1528     c
1529     c--------0---------0---------0---------0---------0---------0---------0--
1530     c
1531     subroutine ctonbe (n,u,v,izs,rzs,dzs)
1532     c
1533     implicit none
1534     integer n,izs(*)
1535     real rzs(*)
1536     double precision u(1),v(1),dzs(*)
1537     c
1538     integer i
1539     c
1540     do i=1,n
1541     v(i)=u(i)
1542     enddo
1543     return
1544     end
1545     c
1546     c--------0---------0---------0---------0---------0---------0---------0--
1547     c
1548     subroutine ctcabe (n,u,v,izs,rzs,dzs)
1549     c
1550     implicit none
1551     integer n,izs(*)
1552     real rzs(*)
1553     double precision u(1),v(1),dzs(*)
1554     c
1555     integer i
1556     c
1557     do i=1,n
1558     v(i)=u(i)
1559     enddo
1560     return
1561     end
1562     c
1563     c--------0---------0---------0---------0---------0---------0---------0--
1564     c
1565     subroutine euclid (n,x,y,ps,izs,rzs,dzs)
1566     c
1567     implicit none
1568     integer n,izs(*)
1569     real rzs(*)
1570     double precision x(n),y(n),ps,dzs(*)
1571     c
1572     integer i
1573     c
1574     ps=0.d0
1575     do i=1,n
1576     ps=ps+x(i)*y(i)
1577     enddo
1578     return
1579     end
1580     c
1581     c--------0---------0---------0---------0---------0---------0---------0--
1582     c
1583     subroutine simul_rc (indic,n,x,f,g,izs,rzs,dzs)
1584     c
1585     implicit none
1586     integer indic,n,izs(*)
1587     real rzs(*)
1588     double precision x(n),f,g(n),dzs(*)
1589     c
1590     return
1591     end
1592     c
1593     c--------0---------0---------0---------0---------0---------0---------0--
1594     c
1595     double precision function dnrmi (n,v)
1596     c
1597     integer n
1598     double precision v(n)
1599     c
1600     c----
1601     c
1602     c Commputes the infinty-norm of the vector v(n)
1603     c
1604     c----
1605     c
1606     c --- local variables
1607     c
1608     integer i
1609     double precision norm
1610     c
1611     c --- compute
1612     c
1613     norm = 0.d0
1614     if (n.gt.0) then
1615     do i=1,n
1616     norm = max(norm,abs(v(i)))
1617     enddo
1618     endif
1619     dnrmi = norm
1620     c
1621     return
1622     end

  ViewVC Help
Powered by ViewVC 1.1.22