1 |
heimbach |
1.1 |
|
2 |
|
|
subroutine lsline( |
3 |
|
|
& simul |
4 |
|
|
& , nn, ifail, lphprint |
5 |
|
|
& , ifunc, nfunc |
6 |
|
|
& , ff, dotdg |
7 |
|
|
& , sflag, tflag, tact, epsx |
8 |
|
|
& , dd, gg, xx, xdiff |
9 |
|
|
& ) |
10 |
|
|
|
11 |
|
|
c ================================================================== |
12 |
|
|
c SUBROUTINE lsline |
13 |
|
|
c ================================================================== |
14 |
|
|
c |
15 |
|
|
c o line search algorithm for determining control vector update; |
16 |
|
|
c After computing updated control vector from given gradient, |
17 |
|
|
c a forward and adjoint model run are performed (simul.F) |
18 |
|
|
c using the updated control vector. |
19 |
|
|
c Tests are then applied to see whether solution has improved. |
20 |
|
|
c |
21 |
|
|
c o Reference: J.C. Gilbert & C. Lemarechal |
22 |
|
|
c Some numerical experiments with variable-storage |
23 |
|
|
c quasi-Newton algorithms |
24 |
|
|
c Mathematical Programming 45 (1989), pp. 407-435 |
25 |
|
|
c |
26 |
|
|
c o started: ??? not reproducible |
27 |
|
|
c |
28 |
|
|
c o changed: Benny Cheng, NASA/JPL |
29 |
|
|
c |
30 |
|
|
c ================================================================== |
31 |
|
|
c SUBROUTINE lsline |
32 |
|
|
c ================================================================== |
33 |
|
|
|
34 |
|
|
#include "blas1.h" |
35 |
|
|
|
36 |
|
|
implicit none |
37 |
|
|
include 'mpif.h' |
38 |
|
|
c---------------------------------- |
39 |
|
|
c declare arguments |
40 |
|
|
c---------------------------------- |
41 |
|
|
integer nmax |
42 |
|
|
parameter( nmax = MAX_INDEPEND ) |
43 |
|
|
integer nn, ifail, ifunc, nfunc |
44 |
|
|
double precision ff, dotdg, tact, epsx |
45 |
|
|
double precision xx(nn), gg(nn), xdiff(nn),dd(nn) |
46 |
|
|
logical tflag,sflag,isforward |
47 |
|
|
logical lphprint |
48 |
|
|
|
49 |
|
|
external simul |
50 |
|
|
|
51 |
|
|
integer pidlen,myindx(2) |
52 |
|
|
integer status(MPI_STATUS_SIZE),ierr |
53 |
|
|
common /mpi_parm/myid,nprocs,mystart,myend |
54 |
|
|
integer myid,nprocs,mystart,myend |
55 |
|
|
c---------------------------------- |
56 |
|
|
c declare local variables |
57 |
|
|
c---------------------------------- |
58 |
|
|
|
59 |
|
|
double precision xpara1, xpara2 |
60 |
|
|
parameter( xpara1 = 0.001D0, xpara2=0.9D0 ) |
61 |
|
|
|
62 |
|
|
double precision factor |
63 |
|
|
parameter( factor = 0.2 ) |
64 |
|
|
|
65 |
|
|
double precision barmax |
66 |
|
|
parameter( barmax = 0.3 ) |
67 |
|
|
double precision barmul |
68 |
|
|
parameter( barmul = 5.0 ) |
69 |
|
|
double precision barmin |
70 |
|
|
parameter( barmin = 0.01 ) |
71 |
|
|
double precision tmax |
72 |
|
|
parameter( tmax = 5. ) |
73 |
|
|
|
74 |
|
|
integer i |
75 |
|
|
|
76 |
|
|
double precision tg, fg, td, ta |
77 |
|
|
double precision fa, fpa, fp |
78 |
|
|
double precision fnew, fdiff,ftemp |
79 |
|
|
double precision z, test, barr |
80 |
|
|
double precision left, right, told |
81 |
|
|
double precision t,a,b |
82 |
|
|
double precision xtemp(nn),gtemp(nn),gold(nn) |
83 |
|
|
external DDOT,DNRM2 |
84 |
|
|
double precision DDOT, DNRM2 |
85 |
|
|
|
86 |
|
|
double precision temp,gnorm0 |
87 |
|
|
real*8,allocatable:: vv(:) |
88 |
|
|
integer, allocatable:: displs(:) |
89 |
|
|
integer, allocatable:: counts(:) |
90 |
|
|
|
91 |
|
|
c---------------------------------- |
92 |
|
|
c check parameters |
93 |
|
|
c---------------------------------- |
94 |
|
|
|
95 |
|
|
if ( (nn.le.0) |
96 |
|
|
& .or. (xpara1.le.0.0) .or. (xpara1.ge.0.5) |
97 |
|
|
& .or. (xpara2.le.xpara1) .or. (xpara2.ge.1.0) ) then |
98 |
|
|
ifail = 9 |
99 |
|
|
go to 999 |
100 |
|
|
endif |
101 |
|
|
|
102 |
|
|
c---------------------------------- |
103 |
|
|
c initialization |
104 |
|
|
c---------------------------------- |
105 |
|
|
|
106 |
|
|
c======================================================================= |
107 |
|
|
c begin of simulation iter. |
108 |
|
|
c======================================================================= |
109 |
|
|
|
110 |
|
|
if (lphprint .and. myid .eq. 0) |
111 |
|
|
& print *, 'lsopt: simul.' |
112 |
|
|
|
113 |
|
|
c------------------------------------ |
114 |
|
|
c compute cost function and gradient |
115 |
|
|
c------------------------------------ |
116 |
|
|
|
117 |
|
|
gold = gg |
118 |
|
|
temp = DNRM2(nn,gg,1) |
119 |
|
|
call MPI_ALLREDUCE(temp*temp,gnorm0,1,MPI_DOUBLE_PRECISION, |
120 |
|
|
& MPI_SUM,MPI_COMM_WORLD,ierr) |
121 |
|
|
gnorm0 = dsqrt(gnorm0) |
122 |
|
|
|
123 |
|
|
if (sflag ) then |
124 |
|
|
|
125 |
|
|
!get estimated cost from hessian update control, this is the trial step |
126 |
|
|
isforward = .true. |
127 |
|
|
call simul( 0, isforward,nn, xdiff, fnew, gg ) |
128 |
|
|
|
129 |
|
|
xtemp = xdiff |
130 |
|
|
! new direction vector dd |
131 |
|
|
|
132 |
|
|
temp = DDOT( nn, xdiff - xx, 1, gold, 1 ) |
133 |
|
|
call MPI_ALLREDUCE(temp,b,1,MPI_DOUBLE_PRECISION,MPI_SUM, |
134 |
|
|
& MPI_COMM_WORLD,ierr) |
135 |
|
|
dotdg = b |
136 |
|
|
|
137 |
|
|
c----------------------------------------- |
138 |
|
|
c quadratic line search |
139 |
|
|
c----------------------------------------- |
140 |
|
|
|
141 |
|
|
a = fnew - ff - b |
142 |
|
|
t = -b/(2.*a) |
143 |
|
|
t = dmax1(t,0.D0) |
144 |
|
|
|
145 |
|
|
if (t <= 1.d-2 ) then |
146 |
|
|
tact = 1. |
147 |
|
|
xdiff = xx - 2.*(.01*ff)*gold*tact/(gnorm0*gnorm0) |
148 |
|
|
tflag = .false. |
149 |
|
|
|
150 |
|
|
if (lphprint .and. myid .eq. 0) |
151 |
|
|
& print *, 'lsopt: using steepest descent, t value ',t,'< .01' |
152 |
|
|
|
153 |
|
|
else |
154 |
|
|
|
155 |
|
|
xdiff = (1. - t)*xx + t*xdiff |
156 |
|
|
tact = t |
157 |
|
|
tflag = .true. |
158 |
|
|
endif |
159 |
|
|
|
160 |
|
|
if (lphprint .and. myid .eq. 0) |
161 |
|
|
& print *, 'lsopt: trial stepsize =',t |
162 |
|
|
|
163 |
|
|
! no hessian update , get new adjoint run |
164 |
|
|
goto 999 |
165 |
|
|
|
166 |
|
|
else |
167 |
|
|
|
168 |
|
|
!no line search |
169 |
|
|
!read new ctrl vector and new cost gradient and new cost |
170 |
|
|
isforward = .false. |
171 |
|
|
call simul( 0, isforward, nn, xdiff, fnew, gg ) |
172 |
|
|
|
173 |
|
|
if (tflag) then |
174 |
|
|
|
175 |
|
|
!read in trial step adjoint output |
176 |
|
|
call simul( ifunc+1, isforward, nn, xtemp, ftemp, gtemp) |
177 |
|
|
|
178 |
|
|
dd = xtemp - xx |
179 |
|
|
c----------------------------------------- |
180 |
|
|
c apply 1st Wolfe test |
181 |
|
|
c----------------------------------------- |
182 |
|
|
fdiff = fnew - ff |
183 |
|
|
|
184 |
|
|
temp = DDOT( nn, dd, 1, gold, 1 ) |
185 |
|
|
call MPI_ALLREDUCE(temp,dotdg,1,MPI_DOUBLE_PRECISION,MPI_SUM, |
186 |
|
|
& MPI_COMM_WORLD,ierr) |
187 |
|
|
|
188 |
|
|
if (fdiff .gt. tact*xpara1*dotdg) then |
189 |
|
|
ifail = 7 |
190 |
|
|
if (lphprint .and. myid .eq. 0) |
191 |
|
|
& print *, 'lsopt: fails wolfe test 1' |
192 |
|
|
|
193 |
|
|
endif |
194 |
|
|
|
195 |
|
|
c----------------------------------------- |
196 |
|
|
c 1st Wolfe test 1 ok, apply 2nd Wolf test |
197 |
|
|
c----------------------------------------- |
198 |
|
|
|
199 |
|
|
temp = DDOT( nn, dd, 1, gg, 1 ) |
200 |
|
|
call MPI_ALLREDUCE(temp,fp,1,MPI_DOUBLE_PRECISION,MPI_SUM, |
201 |
|
|
& MPI_COMM_WORLD,ierr) |
202 |
|
|
|
203 |
|
|
if (fp .lt. xpara2*dotdg) then |
204 |
|
|
ifail = 7 |
205 |
|
|
if (lphprint .and. myid .eq. 0) |
206 |
|
|
& print *, 'lsopt: fails wolfe test 2' |
207 |
|
|
endif |
208 |
|
|
|
209 |
|
|
c----------------------------------------- |
210 |
|
|
c 2nd Wolfe test 2 ok, donc pas serieux, on sort |
211 |
|
|
c----------------------------------------- |
212 |
|
|
|
213 |
|
|
if (ifail .eq. 7 ) then |
214 |
|
|
|
215 |
|
|
if (ifunc .lt. nfunc) then |
216 |
|
|
t = dmin1(tact,tmax) |
217 |
|
|
if (t < 5.D0 ) t = .5*t |
218 |
|
|
xdiff = (1. - t)*xx + t*xtemp |
219 |
|
|
tact = t |
220 |
|
|
ifunc = ifunc + 1 |
221 |
|
|
goto 999 |
222 |
|
|
endif |
223 |
|
|
|
224 |
|
|
endif !end ifail = 7 |
225 |
|
|
|
226 |
|
|
else |
227 |
|
|
|
228 |
|
|
if ( fnew > .999D0*ff ) then |
229 |
|
|
|
230 |
|
|
if (ifunc .lt. nfunc) then |
231 |
|
|
tact = tact*.5 |
232 |
|
|
xdiff = xx - 2.*(.01*ff)*gold*tact/(gnorm0*gnorm0) |
233 |
|
|
ifunc = ifunc + 1 |
234 |
|
|
goto 999 |
235 |
|
|
endif |
236 |
|
|
|
237 |
|
|
endif |
238 |
|
|
|
239 |
|
|
endif ! end of tflag |
240 |
|
|
|
241 |
|
|
dd = xdiff - xx |
242 |
|
|
ff = fnew |
243 |
|
|
xx =xdiff |
244 |
|
|
ifunc = 0 |
245 |
|
|
tact = 1. |
246 |
|
|
endif ! end of sflag |
247 |
|
|
c----------------------------------------- |
248 |
|
|
c end of routine |
249 |
|
|
c----------------------------------------- |
250 |
|
|
999 continue |
251 |
|
|
|
252 |
|
|
return |
253 |
|
|
|
254 |
|
|
end |