1 subroutine dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax, 2 + isave,dsave) 3 character*(*) task 4 integer isave(2) 5 double precision f, g, stp, ftol, gtol, xtol, stpmin, stpmax 6 double precision dsave(13) 7c ********** 8c 9c Subroutine dcsrch 10c 11c This subroutine finds a step that satisfies a sufficient 12c decrease condition and a curvature condition. 13c 14c Each call of the subroutine updates an interval with 15c endpoints stx and sty. The interval is initially chosen 16c so that it contains a minimizer of the modified function 17c 18c psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). 19c 20c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the 21c interval is chosen so that it contains a minimizer of f. 22c 23c The algorithm is designed to find a step that satisfies 24c the sufficient decrease condition 25c 26c f(stp) <= f(0) + ftol*stp*f'(0), 27c 28c and the curvature condition 29c 30c abs(f'(stp)) <= gtol*abs(f'(0)). 31c 32c If ftol is less than gtol and if, for example, the function 33c is bounded below, then there is always a step which satisfies 34c both conditions. 35c 36c If no step can be found that satisfies both conditions, then 37c the algorithm stops with a warning. In this case stp only 38c satisfies the sufficient decrease condition. 39c 40c A typical invocation of dcsrch has the following outline: 41c 42c Evaluate the function at stp = 0.0d0; store in f. 43c Evaluate the gradient at stp = 0.0d0; store in g. 44c Choose a starting step stp. 45c 46c task = 'START' 47c 10 continue 48c call dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax, 49c + isave,dsave) 50c if (task .eq. 'FG') then 51c Evaluate the function and the gradient at stp 52c go to 10 53c end if 54c 55c NOTE: The user must not alter work arrays between calls. 56c 57c The subroutine statement is 58c 59c subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, 60c task,isave,dsave) 61c where 62c 63c stp is a double precision variable. 64c On entry stp is the current estimate of a satisfactory 65c step. On initial entry, a positive initial estimate 66c must be provided. 67c On exit stp is the current estimate of a satisfactory step 68c if task = 'FG'. If task = 'CONV' then stp satisfies 69c the sufficient decrease and curvature condition. 70c 71c f is a double precision variable. 72c On initial entry f is the value of the function at 0. 73c On subsequent entries f is the value of the 74c function at stp. 75c On exit f is the value of the function at stp. 76c 77c g is a double precision variable. 78c On initial entry g is the derivative of the function at 0. 79c On subsequent entries g is the derivative of the 80c function at stp. 81c On exit g is the derivative of the function at stp. 82c 83c ftol is a double precision variable. 84c On entry ftol specifies a nonnegative tolerance for the 85c sufficient decrease condition. 86c On exit ftol is unchanged. 87c 88c gtol is a double precision variable. 89c On entry gtol specifies a nonnegative tolerance for the 90c curvature condition. 91c On exit gtol is unchanged. 92c 93c xtol is a double precision variable. 94c On entry xtol specifies a nonnegative relative tolerance 95c for an acceptable step. The subroutine exits with a 96c warning if the relative difference between sty and stx 97c is less than xtol. 98c On exit xtol is unchanged. 99c 100c task is a character variable of length at least 60. 101c On initial entry task must be set to 'START'. 102c On exit task indicates the required action: 103c 104c If task(1:2) = 'FG' then evaluate the function and 105c derivative at stp and call dcsrch again. 106c 107c If task(1:4) = 'CONV' then the search is successful. 108c 109c If task(1:4) = 'WARN' then the subroutine is not able 110c to satisfy the convergence conditions. The exit value of 111c stp contains the best point found during the search. 112c 113c If task(1:5) = 'ERROR' then there is an error in the 114c input arguments. 115c 116c On exit with convergence, a warning or an error, the 117c variable task contains additional information. 118c 119c stpmin is a double precision variable. 120c On entry stpmin is a nonnegative lower bound for the step. 121c On exit stpmin is unchanged. 122c 123c stpmax is a double precision variable. 124c On entry stpmax is a nonnegative upper bound for the step. 125c On exit stpmax is unchanged. 126c 127c isave is an integer work array of dimension 2. 128c 129c dsave is a double precision work array of dimension 13. 130c 131c Subprograms called 132c 133c MINPACK-2 ... dcstep 134c 135c MINPACK-1 Project. June 1983. 136c Argonne National Laboratory. 137c Jorge J. More' and David J. Thuente. 138c 139c MINPACK-2 Project. November 1993. 140c Argonne National Laboratory and University of Minnesota. 141c Brett M. Averick, Richard G. Carter, and Jorge J. More'. 142c 143c ********** 144 double precision zero, p5, p66 145 parameter (zero=0.0d0,p5=0.5d0,p66=0.66d0) 146 double precision xtrapl, xtrapu 147 parameter (xtrapl=1.1d0,xtrapu=4.0d0) 148 149 logical brackt 150 integer stage 151 double precision finit, ftest, fm, fx, fxm, fy, fym, ginit, gtest, 152 + gm, gx, gxm, gy, gym, stx, sty, stmin, stmax, 153 + width, width1 154 155 external dcstep 156 157c Initialization block. 158 159 if (task(1:5) .eq. 'START') then 160 161c Check the input arguments for errors. 162 163 if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN' 164 if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX' 165 if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO' 166 if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO' 167 if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO' 168 if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO' 169 if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO' 170 if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN' 171 172c Exit if there are errors on input. 173 174 if (task(1:5) .eq. 'ERROR') return 175 176c Initialize local variables. 177 178 brackt = .false. 179 stage = 1 180 finit = f 181 ginit = g 182 gtest = ftol*ginit 183 width = stpmax - stpmin 184 width1 = width/p5 185 186c The variables stx, fx, gx contain the values of the step, 187c function, and derivative at the best step. 188c The variables sty, fy, gy contain the value of the step, 189c function, and derivative at sty. 190c The variables stp, f, g contain the values of the step, 191c function, and derivative at stp. 192 193 stx = zero 194 fx = finit 195 gx = ginit 196 sty = zero 197 fy = finit 198 gy = ginit 199 stmin = zero 200 stmax = stp + xtrapu*stp 201 task = 'FG' 202 203 go to 10 204 205 else 206 207c Restore local variables. 208 209 if (isave(1) .eq. 1) then 210 brackt = .true. 211 else 212 brackt = .false. 213 end if 214 stage = isave(2) 215 ginit = dsave(1) 216 gtest = dsave(2) 217 gx = dsave(3) 218 gy = dsave(4) 219 finit = dsave(5) 220 fx = dsave(6) 221 fy = dsave(7) 222 stx = dsave(8) 223 sty = dsave(9) 224 stmin = dsave(10) 225 stmax = dsave(11) 226 width = dsave(12) 227 width1 = dsave(13) 228 229 end if 230 231c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the 232c algorithm enters the second stage. 233 234 ftest = finit + stp*gtest 235 if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) stage = 2 236 237c Test for warnings. 238 239 if (brackt .and. (stp .le. stmin .or. stp .ge. stmax)) 240 + task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS' 241 if (brackt .and. stmax-stmin .le. xtol*stmax) 242 + task = 'WARNING: XTOL TEST SATISFIED' 243 if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest) 244 + task = 'WARNING: STP = STPMAX' 245 if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest)) 246 + task = 'WARNING: STP = STPMIN' 247 248c Test for convergence. 249 250 if (f .le. ftest .and. abs(g) .le. gtol*(-ginit)) 251 + task = 'CONVERGENCE' 252 253c Test for termination. 254 255 if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') go to 10 256 257c A modified function is used to predict the step during the 258c first stage if a lower function value has been obtained but 259c the decrease is not sufficient. 260 261 if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then 262 263c Define the modified function and derivative values. 264 265 fm = f - stp*gtest 266 fxm = fx - stx*gtest 267 fym = fy - sty*gtest 268 gm = g - gtest 269 gxm = gx - gtest 270 gym = gy - gtest 271 272c Call dcstep to update stx, sty, and to compute the new step. 273 274 call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,brackt,stmin, 275 + stmax) 276 277c Reset the function and derivative values for f. 278 279 fx = fxm + stx*gtest 280 fy = fym + sty*gtest 281 gx = gxm + gtest 282 gy = gym + gtest 283 284 else 285 286c Call dcstep to update stx, sty, and to compute the new step. 287 288 call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,brackt,stmin,stmax) 289 290 end if 291 292c Decide if a bisection step is needed. 293 294 if (brackt) then 295 if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty-stx) 296 width1 = width 297 width = abs(sty-stx) 298 end if 299 300c Set the minimum and maximum steps allowed for stp. 301 302 if (brackt) then 303 stmin = min(stx,sty) 304 stmax = max(stx,sty) 305 else 306 stmin = stp + xtrapl*(stp-stx) 307 stmax = stp + xtrapu*(stp-stx) 308 end if 309 310c Force the step to be within the bounds stpmax and stpmin. 311 312 stp = max(stp,stpmin) 313 stp = min(stp,stpmax) 314 315c If further progress is not possible, let stp be the best 316c point obtained during the search. 317 318 if (brackt .and. (stp .le. stmin .or. stp .ge. stmax) .or. 319 + (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx 320 321c Obtain another function and derivative. 322 323 task = 'FG' 324 325 10 continue 326 327c Save local variables. 328 329 if (brackt) then 330 isave(1) = 1 331 else 332 isave(1) = 0 333 end if 334 isave(2) = stage 335 dsave(1) = ginit 336 dsave(2) = gtest 337 dsave(3) = gx 338 dsave(4) = gy 339 dsave(5) = finit 340 dsave(6) = fx 341 dsave(7) = fy 342 dsave(8) = stx 343 dsave(9) = sty 344 dsave(10) = stmin 345 dsave(11) = stmax 346 dsave(12) = width 347 dsave(13) = width1 348 349 end 350