 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:06:29 00004 00031 00033 MODULE linesrch 00034 IMPLICIT NONE 00035 00036 INTEGER, PARAMETER :: msfd=20 00037 INTEGER :: nsfd !< number of function calls 00038 INTEGER :: idgl !< index of smallest negative slope 00039 INTEGER :: idgr !< index of smallest positive slope 00040 INTEGER :: idgm !< index of minimal slope 00041 INTEGER :: minf=1 !< min. number of function calls 00042 INTEGER :: maxf=5 !< max. number of function calls 00043 INTEGER :: lsinfo !< (status) information 00044 DOUBLE PRECISION, DIMENSION(4,msfd) :: sfd !< abscissa; function value; slope; predicted zero 00045 DOUBLE PRECISION :: stmx=0.9 !< maximum slope ratio 00046 DOUBLE PRECISION :: gtol !< slope ratio 00047 00048 END MODULE linesrch 00049 00067 00068 SUBROUTINE ptline(n,x,f,g,s,step, info) ! - 2 arguments 00069 USE linesrch 00070 00071 IMPLICIT NONE 00072 INTEGER, INTENT(IN) :: n 00073 DOUBLE PRECISION, INTENT(IN OUT) :: x(n) 00074 DOUBLE PRECISION, INTENT(IN OUT) :: f 00075 DOUBLE PRECISION, INTENT(IN OUT) :: g(n) 00076 DOUBLE PRECISION, INTENT(IN OUT) :: s(n) 00077 DOUBLE PRECISION, INTENT(OUT) :: step 00078 INTEGER, INTENT(OUT) :: info 00079 00080 INTEGER:: i1 00081 INTEGER :: i2 00082 INTEGER :: i ! internal 00083 INTEGER :: im ! internal 00084 DOUBLE PRECISION :: alpha ! internal 00085 DOUBLE PRECISION :: dginit ! internal 00086 DOUBLE PRECISION :: dg ! internal 00087 DOUBLE PRECISION :: fsaved ! internal 00088 DOUBLE PRECISION :: tot ! internal 00089 DOUBLE PRECISION :: fp1 ! internal 00090 DOUBLE PRECISION :: fp2 ! internal 00091 SAVE 00092 00093 ! initialization --------------------------------------------------- 00094 00095 info=0 ! reset INFO flag 00096 dg=0.0D0 00097 DO i=1,n ! 00098 dg=dg-g(i)*s(i) ! DG = scalar product: grad x search 00099 END DO! 00100 00101 IF(nsfd == 0) THEN ! initial call 00102 dginit=dg ! DG = initial directional gradient 00103 IF(dginit >= 0.0D0) GO TO 100 ! error: step not decreasing 00104 step=1.0D0 ! initial step factor is one 00105 alpha=step ! get initial step factor 00106 tot=0.0D0 ! reset total step 00107 idgl=1 ! index of smallest negative slope 00108 idgr=0 ! index of smallest positive slope 00109 fsaved=f ! initial Function value 00110 nsfd=1 ! starting point of iteration 00111 sfd(1,1)=0.0 ! abscissa 00112 sfd(2,1)=0.0 ! reference function value 00113 sfd(3,1)=dginit ! slope 00114 sfd(4,1)=0.0 ! predicted zero 00115 im=1 ! optimum 00116 ELSE ! subsequent call 00117 nsfd=nsfd+1 00118 sfd(1,nsfd)=tot ! abscissa 00119 sfd(2,nsfd)=f-fsaved ! function value difference to reference 00120 sfd(3,nsfd)=dg ! slope 00121 sfd(4,nsfd)=0.0 ! predicted zero (see below) 00122 IF(dg < sfd(3,im)) THEN 00123 im=nsfd 00124 END IF 00125 00126 ! define interval indices IDGL and IDGR 00127 IF(dg <= 0.0D0) THEN 00128 IF(dg >= sfd(3,idgl)) idgl=nsfd 00129 END IF 00130 IF(dg >= 0.0D0) THEN ! limit to the right 00131 IF(idgr == 0) idgr=nsfd 00132 IF(dg <= sfd(3,idgr)) idgr=nsfd 00133 END IF 00134 00135 IF(idgr == 0) THEN 00136 i1=nsfd-1 00137 i2=nsfd 00138 ELSE 00139 i1=idgl 00140 i2=idgr 00141 END IF 00142 fp1=sfd(3,i1) 00143 fp2=sfd(3,i2) ! interpolation 00144 sfd(4,nsfd)=(sfd(1,i1)*fp2-sfd(1,i2)*fp1)/(fp2-fp1) 00145 00146 ! convergence tests 00147 IF(nsfd >= minf.AND.ABS(dg) <= ABS(dginit)*gtol) THEN 00148 ! normal convergence return with INFO=1 ---------------------- 00149 alpha=tot+alpha ! total ALPHA is returned 00150 step =alpha 00151 idgm=idgl 00152 IF(idgr /= 0) THEN 00153 IF(sfd(3,idgr)+sfd(3,idgl) < 0.0D0) idgm=idgr 00154 END IF 00155 GO TO 101 00156 END IF 00157 IF(nsfd >= maxf) GO TO 102 ! max number of function calls 00158 alpha=MIN(sfd(4,nsfd),stmx)-tot ! new step from previous 00159 IF(ABS(alpha) < 1.0D-3.AND.sfd(4,nsfd) > stmx) GO TO 103 00160 IF(ABS(alpha) < 1.0D-3) GO TO 104 00161 END IF 00162 00163 ! prepare next function call --------------------------------------- 00164 00165 DO i=1,n 00166 x(i)=x(i)+alpha*s(i) ! step by ALPHA -> new X 00167 END DO 00168 tot=tot+alpha ! 00169 step=tot 00170 info=-1 ! recalculate function and gradient 00171 lsinfo=info 00172 RETURN 00173 00174 ! error exits ------------------------------------------------------ 00175 104 info=info+1 ! 4: step small 00176 103 info=info+1 ! 3: maximum reached 00177 102 info=info+1 ! 2: too many function calls 00178 101 info=info+1 ! 1: normal convergence 00179 lsinfo=info 00180 im=1 00181 DO i=1,nsfd 00182 IF(ABS(sfd(3,i)) < ABS(sfd(3,im))) im=i 00183 END DO 00184 alpha=sfd(1,im)-sfd(1,nsfd) 00185 IF(im == nsfd) RETURN ! already at minimum 00186 DO i=1,n 00187 x(i)=x(i)+alpha*s(i) ! step by ALPHA to slope minimum 00188 END DO 00189 f=sfd(2,im)+fsaved ! F at minimum 00190 step=sfd(1,im) ! total step at convergence 00191 IF(im /= 1) RETURN ! improvement 00192 info=5 ! no improvement 00193 100 step=0.0D0 ! 0: initial slope not negative 00194 lsinfo=info 00195 RETURN 00196 END SUBROUTINE ptline 00197 00210 00211 SUBROUTINE ptldef(gtole,stmax,minfe,maxfe) 00212 USE linesrch 00213 00214 IMPLICIT NONE 00215 INTEGER, INTENT(IN) :: minfe 00216 INTEGER, INTENT(IN) :: maxfe 00217 REAL, INTENT(IN) :: gtole 00218 REAL, INTENT(IN) :: stmax 00219 00220 gtol=MAX(1.0E-4,MIN(gtole,0.9E0)) ! slope ratio 00221 IF(gtole == 0.0) gtol=0.9D0 ! default slope ratio 00222 stmx=stmax ! maximum total step 00223 IF(stmx == 0.0D0) stmx=10.0D0 ! default limit 00224 minf=MAX(1,MIN(minfe,msfd-2)) ! minimum number of evaluations 00225 maxf=MAX(2,MIN(maxfe,msfd-1)) ! maximum number of evaluations 00226 IF(maxfe == 0) maxf=5 ! default max number of values 00227 nsfd=0 ! reset 00228 END SUBROUTINE ptldef 00229 00236 00237 SUBROUTINE ptlopt(nf,m,slopes,steps) 00238 USE linesrch 00239 IMPLICIT NONE 00240 00241 INTEGER, INTENT(OUT) :: nf 00242 INTEGER, INTENT(OUT) :: m 00243 REAL, DIMENSION(3), INTENT(OUT) :: slopes 00244 REAL, DIMENSION(3), INTENT(OUT) :: steps 00245 INTEGER :: i 00246 00247 ! ... 00248 nf=nsfd 00249 IF(nsfd == 0) THEN ! no values 00250 m=0 00251 DO i=1,3 00252 slopes(i)=0.0 00253 steps(i) =0.0 00254 END DO 00255 ELSE ! values exist 00256 m=1 00257 DO i=1,nsfd 00258 IF(ABS(sfd(3,i)) < ABS(sfd(3,m))) m=i 00259 END DO 00260 slopes(1)=REAL(sfd(3,1)) 00261 slopes(2)=REAL(sfd(3,nsfd)) 00262 slopes(3)=REAL(sfd(3,m)) 00263 steps(1) =REAL(sfd(1,1)) 00264 steps(2) =REAL(sfd(1,nsfd)) 00265 steps(3) =REAL(sfd(1,m)) 00266 END IF 00267 END SUBROUTINE ptlopt 00268 00272 00273 SUBROUTINE ptlprt(lunp) 00274 USE linesrch 00275 00276 IMPLICIT NONE 00277 INTEGER :: i 00278 INTEGER :: j 00279 INTEGER :: im 00280 INTEGER :: lun 00281 INTEGER, INTENT(IN) :: lunp 00282 REAL :: ratio 00283 CHARACTER (LEN=2) :: tlr 00284 ! ... 00285 lun=lunp 00286 IF(lun == 0) lun=6 00287 IF(nsfd <= 0) RETURN 00288 WRITE(lun,*) ' ' 00289 WRITE(lun,*) 'PTLINE: line-search method based on slopes', & 00290 ' with sufficient slope-decrease' 00291 WRITE(lun,*) 'PTLDEF: slope ratio limit=',gtol 00292 WRITE(lun,*) 'PTLDEF: maximum step =',stmx 00293 WRITE(lun,*) 'PTLDEF:',minf,' <= nr of calls <=',maxf 00294 WRITE(lun,101) 00295 im=1 00296 DO i=1,nsfd 00297 IF(ABS(sfd(3,i)) < ABS(sfd(3,im))) im=i 00298 END DO 00299 DO i=1,nsfd 00300 tlr=' ' 00301 IF(i == im) tlr='**' 00302 IF(i == idgl) tlr(1:1)='L' 00303 IF(i == idgr) tlr(2:2)='R' 00304 IF(i == 1) THEN 00305 WRITE(lun,102) i-1, sfd(1,i),tlr,(sfd(j,i),j=2,4) 00306 ELSE 00307 ratio=REAL(ABS(sfd(3,i)/sfd(3,1))) 00308 WRITE(lun,103) i-1, sfd(1,i),tlr,(sfd(j,i),j=2,4),ratio 00309 END IF 00310 00311 END DO 00312 IF(lsinfo == 0) WRITE(lun,*) & 00313 'PTLINE: INFO=0 input error (e.g. gradient not negative)' 00314 IF(lsinfo == 1) WRITE(lun,*) 'PTLINE: INFO=1 convergence reached' 00315 IF(lsinfo == 2) WRITE(lun,*) 'PTLINE: INFO=2 too many function calls' 00316 IF(lsinfo == 3) WRITE(lun,*) 'PTLINE: INFO=3 maximum step reached' 00317 IF(lsinfo == 4) WRITE(lun,*) 'PTLINE: INFO=4 step too small (< 0.001)' 00318 WRITE(lun,*) ' ' 00319 00320 101 FORMAT(' i x F(x) F''(X)', & 00321 ' minimum F''(X)') 00322 102 FORMAT(i3,f12.6,1X,a2,g15.6,g14.6,f12.6,' ratio') 00323 103 FORMAT(i3,f12.6,1X,a2,g15.6,g14.6,f12.6,f10.3) 00324 00325 END SUBROUTINE ptlprt 00326
1.8.0
